FIR and spectrum tools for VC2005 Express

All replies

• 

  

  0

  Sign in to vote

  Program for converting h(t) coefficients of the FIR  to AFR  data

  // irtoafr.cpp : Defines the entry point for the console application.
  //

  #include "stdafx.h"

  #include <cstdio> // or  cstdio.h
  #include <cstdlib> //or  cstdlib.h
  #include <cmath>   //or  cmath.h

  #define M_PI 3.1415926535897932384626433832795
  #define PI M_PI // 3.14159265
  #define  NMAX 512
  FILE *fp, *fp1;
  int w,t,N;
  double h[NMAX];
  double h_im[NMAX];
  double Re_H[NMAX];
  double Im_H[NMAX];
  double Re_x[NMAX];
  double Im_x[NMAX];
  double abs_H[NMAX];
  double arg_h[NMAX];
  double f[NMAX];
  double angle;
  double   K, fs,T,Tmax;
  char str;
  char key;
  int k,n;

  
  /*
  float x[N];
  float y[N];

  double my_FIR(double sample_data, int Nf)
  {
    double result = 0;
    for ( int i = Nf - 2 ; i >= 0 ; i-- )
    {
      x[i + 1] = x[i];
      y[i + 1] = y[i];
    }
    x[0] = (double)sample_data;
    for (int k = 0; k < N; k++)
    {
      result = result + x[k]*h[k];
    }
    y[0] = result;
    return ((double)result);
  }

  */

  
  float heavi( double t)
  {
  if (t<0 ) return 0;
  if (t>=0) return 1;
  }

  float heavi1( int n)
  {
  if (n<0 ) return 0;
  if (n>=0) return 1;
  }

  
  float dirac( int n)
  {
  if (n!=0 ) return 0;
  if (n=0) return 1;
  }

  double getarg(float Re_H,float Im_H)
  {
  return atan2(Im_H,Re_H) * 180/ PI;
  }

  
  int _tmain(int argc, _TCHAR* argv[])
  {
  //double tau=0.02;
  //printf(" Input tau , s ");
  //scanf("%le",&tau);
  printf ("\n For input  h(t) from file input '1' ");
  printf ("\n For manual  input  h(t) input   '0' ");
  scanf("%s",&key);
  if (key=='0' ) { goto label_02;}
  if (key=='1' ) { goto label_01;}

  label_01:
    printf("\n Open htkoefs.txt for reading h[n] coefficients ");
  if ((fp=fopen("htkoefs.txt","r"))==NULL)  printf ("\n File opening error");
  fscanf(fp,"%le",&fs);
  fscanf(fp,"%d",&N);
  Tmax=N*1/fs;
  for (n=0 ; n<N;n++)
  {
   fscanf (fp,"%d",&k);
   fscanf (fp,"%le",&h[n]);
  }
  fclose(fp);
  goto label_03; 
  label_02:
   
   
  printf("\n Input fs ,Hz ");
  scanf("%le",&fs);
  printf(" Input  number of points ");
  scanf("%d",&N);
  Tmax=N*1/fs;
  printf("\n Tmax =%lf s \n", Tmax); 
   
  for (n=0 ; n<N;n++)
  {
   T=n*1/(N*fs); 
   printf("\nT=%le s; n =%d ; input h(n): ",T,n);
   scanf ("%le",&h[n]);

  /*
  h[0]=-0.0485469;
  h[1]=+0.0307880;
  h[2]=+0.0678165;
  h[3]=-0.0079250;
  h[4]=-0.0980449;
  h[5]=-0.0384873;
  h[6]=+0.1898828;
  h[7]=+0.4045168;
  h[8]=+0.4045168;
  h[9]=+0.1898828;
  h[10]=-0.0384873;
  h[11]=-0.0980449;
  h[12]=-0.0079250;
  h[13]=+0.0678165;
  h[14]=+0.0307880;
  h[15]=-0.0485469;
  */
  //tau=0.1*Tmax;
  //h[t]=1-exp(-T/tau);
  }

  label_03:
  printf("\n Fs=%lf Hz",fs);
  printf("\n N=%d  points",N);
  for (n=0 ; n<N;n++)
  {
   T=n*1/(N*fs); 
   printf("\n t=%le s; h[t]=%le   ",T,h[n]);
  }

  printf("\n Input 0 for next ");
  scanf("%s",&key);
  printf("\n OK. Parsing ... ");
   
  // H(jw)=integr(0;inf;h(tau)*exp(-j*w*tau) ; d tau)
  //
  /*
  DFT
  X[k]=sum(n=0;N-1; x[n]*exp(-2*pi*j*k*n/N) ), k=0,...,N-1
  IDFT
  x[n]=(1/N)*sum(k=0;N-1;X[k]*exp(2*pi*k*n/N) ), n=0,...,N-1
  */

  /*********************************************************/

  label1:
   printf("\n Open afrcoefs.txt for reading AFR coefficients ");
  if ((fp1=fopen("afrcoefs.txt","a"))==NULL)  printf ("\n File opening error");
  fprintf(fp1,"\n   %le  %d ",fs,N);
  /*
  fprintf(fp1,"\n Re(H(jw))  ; Im(H(jw)) ");
  for (k=0 ; k<N; k++)
  {
   Re_H[k]=0;
   Im_H[k]=0;
    for(n=0;n<N;n++)
    {    
     angle=-2*PI*k*n/N;
   //  Re_H[k] = Re_H[k]+ h[n]*cos(angle) - h_im[n]*sin(angle);
   //  Im_H[k] = Im_H[k]+ h[n]*sin(angle) + h_im[n]*cos(angle);
       Re_H[k] = Re_H[k]+ h[n]*cos(angle) ;//- h_im[n]*sin(angle);
       Im_H[k] = Im_H[k]+ h[n]*sin(angle) ;//+ h_im[n]*cos(angle);
    }

  // 1/n scaled DFT(h(t)*1(t))
   //Re_H[k]=Re_H[k]/N ;
   //Im_H[k]=Im_H[k]/N ;
     abs_H[k]=pow(((Re_H[k]*Re_H[k])+(Im_H[k]*Im_H[k])), 0.5);
     arg_h[k]=getarg(Re_H[k],Im_H[k]);
     f[k]=k*fs/N; 
       printf("\n n=%d Re(H(jw))= %lf ; Im(H(jw))= %lf ; f=%lf Hz",k, Re_H[k] ,Im_H[k], f[k]) ;
       fprintf(fp1,"\n n=%d Re(H(jw))= %lf ; Im(H(jw))= %lf ; f=%lf Hz",k, Re_H[k] ,Im_H[k], f[k]) ;
  }

  printf("\n Input 0 for continue ");
  scanf("%s",&key);
  */
   //fprintf(fp1,"\n |H(jw)| , arg( H(jw) ) ");
   
  for (k=0 ; k<N; k++)
  {
   Re_H[k]=0;
   Im_H[k]=0;
    for(n=0;n<N;n++)
    {   
     angle=-2*PI*k*n/N;
   //  Re_H[k] = Re_H[k]+ h[n]*cos(angle) - h_im[n]*sin(angle);
   //  Im_H[k] = Im_H[k]+ h[n]*sin(angle) + h_im[n]*cos(angle);
       Re_H[k] = Re_H[k]+ h[n]*cos(angle) ;//- h_im[n]*sin(angle);
       Im_H[k] = Im_H[k]+ h[n]*sin(angle) ;//+ h_im[n]*cos(angle);
    }

   // 1/n scaled DFT(h(t)*1(t))
   //Re_H[k]=Re_H[k]/N ;
   //Im_H[k]=Im_H[k]/N ;
   
     abs_H[k]=pow(((Re_H[k]*Re_H[k])+(Im_H[k]*Im_H[k])), 0.5);
     arg_h[k]=getarg(Re_H[k],Im_H[k]);
     f[k]=k*fs/N;
     printf("\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
      //fprintf(fp1,"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
     fprintf(fp1,"\n  %d   %le    %le    %le  ",k, abs_H[k] ,arg_h[k], f[k]) ;
  }
  fclose(fp1);
  printf("\n Input 0 for exit ");
  scanf("%s",&key);

  return 0;
  }

  //cidft
  /*
  for (n=0 ; n<N; n++)
  {
   Re_x[n]=0;
   Im_x[n]=0;
    for(k=0;k<N;n++)
    {  
  //  angle = 2 * PI * k * N / numidft
  //  idft_xnre(N + 1) = idft_xnre(N + 1) + idft_Xkre(k + 1) * Cos(angle) - idft_Xkim(k + 1) * Sin(angle)
  //  idft_xnim(N + 1) = idft_xnim(N + 1) + idft_Xkre(k + 1) * Sin(angle) + idft_Xkim(k + 1) * Cos(angle) 
     angle=2*PI*k*n/N;
     Re_H[k]=abs_H[k]*cos(angle);
     Im_H[k]=abs_H[k]*sin(angle);
      Re_x[n] = Re_x[n]+ Re_H[k]*cos(angle)- Im_H[k]*sin(angle);
      Im_x[n] = Im_x[n]+ Re_H[k]*sin(angle)+ Im_H[k]*cos(angle);
    }
   Re_x[n]=Re_x[n]/N ;
   Im_x[n]=Im_x[n]/N ;
   // n scaled CIDFT( )
   //Re_x[n]=Re_x[n]*N;  
   //Im_x[n]=Im_x[n]*N;  
   //  abs_H[k]=pow(((Re_H[k]*Re_H[k])+(Im_H[k]*Im_H[k])), 0.5);
   //  arg_h[k]=getarg(Re_H[k],Im_H[k]);
   //  T[k]=k*fs/N;
     printf("\n n=%d  Re(h(n))=%d   ",n, Re_x[n] ,Im_x[n]  ) ;

  }

  */

  // stdafx.cpp : source file that includes just the standard includes
  // irtoafr.pch will be the pre-compiled header
  // stdafx.obj will contain the pre-compiled type information

  #include "stdafx.h"

  // TODO: reference any additional headers you need in STDAFX.H
  // and not in this file

  // stdafx.h : include file for standard system include files,
  // or project specific include files that are used frequently, but
  // are changed infrequently
  //

  #pragma once

  
  #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
  #include <stdio.h>
  #include <tchar.h>

  // TODO: reference additional headers your program requires here

  Tuesday, January 23, 2018 10:01 PM
• 

  

  0

  Sign in to vote

  Implementation of the FIR

  // fir.cpp : Defines the entry point for the console application.
  //

  #include "stdafx.h"

  #include <cstdio> // or  cstdio.h
  #include <cstdlib> //or  cstdlib.h
  #include <cmath>   //or  cmath.h

   
  #define  NMAX 512
  FILE *fp, *fp1, *fp2;
  int w,t,N,N1;
  double h[NMAX];
  double xin[NMAX];
  double yout[NMAX];
  double f[NMAX];
   
  double   K, fs,T,Tmax;
  char str;
  char key;
  int k,n,i;

  
  double my_FIR(double sample_data, int Nf,  double *h)
  {
  static double x[NMAX];
  static double y[NMAX];

    double result = 0;
    for ( int i = Nf - 2 ; i >= 0 ; i-- )
    {
      x[i + 1] = x[i];
      y[i + 1] = y[i];
    }
    x[0] = (double)sample_data;
    for (int k = 0; k < Nf; k++)
    {
      result = result + x[k]*h[k];
    }
    y[0] = result;
    return ((double)result);
  }

   

  
  float heavi( double t)
  {
  if (t<0 ) return 0;
  if (t>=0) return 1;
  }

  float heavi1( int n)
  {
  if (n<0 ) return 0;
  if (n>=0) return 1;
  }

  
  float dirac( int n)
  {
  if (n!=0 ) return 0;
  if (n=0) return 1;
  }

   

  
  int _tmain(int argc, _TCHAR* argv[])
  {
   
  printf ("\n For input  h(t) from file input '1' ");
  printf ("\n For manual  input  h(t) input   '0' ");
  scanf("%s",&key);
  if (key=='0' ) { goto label_02;}
  if (key=='1' ) { goto label_01;}

  label_01:

   printf("\n Open htkoefs.txt for reading h[n] coefficients ");
  if ((fp=fopen("htkoefs.txt","r"))==NULL)  printf ("\n File opening error");
  fscanf(fp,"%le",&fs);
  fscanf(fp,"%d",&N);

  for (n=0 ; n<N;n++)
  {
   fscanf (fp,"%d",&k);
   fscanf (fp,"%le",&h[n]);
  }

  
  fclose(fp);
  goto label_03;
   
  label_02:
   
   
  printf("\n Input fs ,Hz " );
  scanf("%le",&fs);
  printf(" Input order ");
  scanf("%d",&N);
    
   
  for (n=0 ; n<N;n++)
  {
   
   printf("\n   n =%d ; input h[n]: " ,n);
   scanf ("%le",&h[n]);

  }

  label_03:
  printf("\n Fs coefs=%lf Hz ",fs);
  printf("\n order=  %d     ",N);

  for (n=0 ; n<N;n++)
  {
   printf("\n  n =%d ; h[n]=%le   ",n,h[n]);
  }

  printf("\n Input 0 for continue ");
  scanf("%s",&key);
  printf("\n OK.   ");

  
  printf ("\n For input  x(t) from file input '1' ");
  printf ("\n For manual  input  x(t) input   '0' ");
  scanf("%s",&key);
  if (key=='1' ) { goto label_05;}
  if (key=='0' ) { goto label_06;}

  label_05:

   printf("\n Open input.txt for reading x[n]   ");
  if ((fp1=fopen("input.txt","r"))==NULL)  printf ("\n File opening error ");
  fscanf(fp1,"%le",&fs);
  fscanf(fp1,"%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   fscanf (fp1,"%d",&k);
   fscanf (fp1,"%le",&xin[n]);
  }

  
  fclose(fp1);

  goto label_07;

  label_06:

  printf("\n Input fs,Hz , ( %lf ) ",fs) ;
  scanf( "%le",&fs);
  printf("\n Input number of samples of xin[n]  ") ;
  scanf( "%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   T=n*1/(N1*fs); 
   printf ("\n n=%d ; t=%le s; xin[n]= ",n,T);
   scanf ( "%le",&xin[n]);
  }

   label_07:
  printf( "\n fs=%le Hz ",fs);
  printf( "\n Nsamples= %d ",N1);

  for (n=0 ; n<N1;n++)
  {
   T=n*1/(N1*fs); 
   printf("\n n=%d x[n]=%le,  t=%le s",n,xin[n],T);
  }

  printf("\n Input 0 for continue ");
  scanf("%s",&key);
  printf("\n OK.   ");

  /*********************************************************/

  label1:
    printf("\n Open output.txt for writing sequence   ");

  if ((fp2=fopen("output.txt","a"))==NULL)  printf ("\n File opening error");
  fprintf(fp2,"\n   %le  %d \n",fs,N);

  for (i=0 ;  i<N1;i++)
  {
   T=i*1/(N1*fs);
  yout[i]= my_FIR(xin[i], N,h);
   printf( "\n n=%d t=%le s; yout[n]= %le ",i,T,yout[i] );
  fprintf(fp2,"\n %d  %le ",i,yout[i] );
  }

  
  fclose(fp2);
  printf("\n Input 0 for exit ");
  scanf("%s",&key);

  return 0;
  }

   

  // stdafx.cpp : source file that includes just the standard includes
  // fir.pch will be the pre-compiled header
  // stdafx.obj will contain the pre-compiled type information

  #include "stdafx.h"

  // TODO: reference any additional headers you need in STDAFX.H
  // and not in this file

  // stdafx.h : include file for standard system include files,
  // or project specific include files that are used frequently, but
  // are changed infrequently
  //

  #pragma once

  
  #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
  #include <stdio.h>
  #include <tchar.h>

  // TODO: reference additional headers your program requires here

  input.txt
     4.410000e+004   16

   0  1.000000e+001  
   1  7.071068e+000  
   2  6.123032e-016  
   3  -7.071068e+000  
   4  -1.000000e+001  
   5  -7.071068e+000  
   6  -1.836910e-015  
   7  7.071068e+000  
   8  1.000000e+001  
   9  7.071068e+000  
   10  3.061516e-015  
   11  -7.071068e+000  
   12  -1.000000e+001  
   13  -7.071068e+000  
   14  -4.286122e-015  
   15  7.071068e+000 

  htkoefs.txt

    4.410000e+004   16
    0   0 
    1   1 
    2   0  
    3   0  
    4   0  
    5   0  
    6   0 
    7   0  
    8   0  
    9   0  
    10  0  
    11  0  
    12  0  
    13  0  
    14  0  
    15  0 

  output.txt

  
     4.410000e+004  16

   0  0.000000e+000
   1  1.000000e+001
   2  7.071068e+000
   3  6.123032e-016
   4  -7.071068e+000
   5  -1.000000e+001
   6  -7.071068e+000
   7  -1.836910e-015
   8  7.071068e+000
   9  1.000000e+001
   10  7.071068e+000
   11  3.061516e-015
   12  -7.071068e+000
   13  -1.000000e+001
   14  -7.071068e+000
   15  -4.286122e-015

  
  

  • Edited by USERPC01 Tuesday, January 23, 2018 10:08 PM fixed bugs

  Tuesday, January 23, 2018 10:07 PM
• 

  

  0

  Sign in to vote

  Program for spectrum calculating, using CDFT,CIDFT, Mathlab's algorithm for Hilbert transform

  You can use complex.h, vector.h for data arrays and pointered arrays for arrays.

  How to use radix more then 20000 in  #define NMAX  , how to fix  memory mapping error in arrays if  NMAX>20000...32768?

  // getspec.cpp : Defines the entry point for the console application.
  //

  #include "stdafx.h"

  #include <cstdio>
  #include <cstdlib>
  #include <cmath>
  #include <limits>
  #include <iostream>
  using namespace std;

  #define NMAX 2048

  #define M_PI 3.1415926535897932384626433832795

  typedef struct {

  double re[NMAX];
  double im[NMAX];
  } CData;

   typedef struct  {
       double data[NMAX];
   } Array;

  
  CData  cdft(CData xn  ,int N )
   {
     
   
   CData Xk ;
   int k,n  ;
   double angle;

   for (k=0;k<N;k++)
    {
    Xk.re[k]=0;
    Xk.im[k]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*M_PI*k*n/N;  
      Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
      Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
      }

   }

   return Xk;
   }

   CData  cidft(CData idft_Xk, int N  )
   {

   int k,n ;
   double angle;
   CData  idft_xn;
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    idft_xn.re[n]=0;
    idft_xn.im[n]=0;
    for(k=0;k<N;k++)
     {
   
     angle=2*M_PI*k*n/N;
     idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
     idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
     }

   idft_xn.im[n]=idft_xn.im[n]/N;
   idft_xn.re[n]=idft_xn.re[n]/N;

   }

  return idft_xn ;
   }

  CData ScaleDFT(CData xk,int  num)
   {

   int i;
   for (i=0; i<num ; i++)
     {
      xk.re[i]=xk.re[i]/num;
      xk.im[i]=xk.im[i]/num;
     }
   return xk;
   }

  CData ScaleIQidft(CData Xk, int num)
   {
   int i;
     for (i=0;  i<num ; i++)
    {
     Xk.re[i]=Xk.re[i]*num;
     Xk.im[i]=Xk.im[i]*num;
    }

  return Xk;

  }

   CData hilbert(double *xr , int N)
   {
    int  k,n,i, no2 ;
    double angle;
    double h[NMAX];
    double dft__xn_re[NMAX];
    double dft__xn_im[NMAX];
     double dft__Xk_re[NMAX];
     double dft__Xk_im[NMAX];
       double idft__XK_re[NMAX];
       double idft__XK_im[NMAX];
     CData idft__XN;

   
    if (N==0)
    {
    printf("\n Error");
    idft__XN.re[0]=0;
    idft__XN.im[0]=0;
    return idft__XN;
    }
    //xnreal=real(xr);
    no2=(int) floor ((double)N/2);
   
    for (i=0;i<N;i++)
    {
     dft__xn_re[i]=xr[i];
     dft__xn_im[i]=0;
    }

   //fft
   // dft__Xk=dft(dft__xn,n);

  for (k=0;k<N;k++)
    {
    dft__Xk_re[k]=0;
    dft__Xk_im[k]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*M_PI*k*n/N;
      dft__Xk_re[k]+= dft__xn_re[n]*cos(angle)- dft__xn_im[n]*sin(angle);
      dft__Xk_im[k]+= dft__xn_re[n]*sin(angle)+ dft__xn_im[n]*cos(angle);
      }

   }

   //X={Xk.re,Xk.im};

   for (i=0;i<=N;i++)
    {
    h[i]=0;
    }

  
    if((N>0)&&(( 2*floor((double)N/2) )==(N))) // n is even
    {

    h[0]=1;
     h[N/2]=1;
     for(k=1;k<N/2;k++)
     {
     h[k]=2;
     }
    }
    else //n is odd
    {
      if(N>0)
      {
       h[0]=1;
       for(k=1;k<(N+1)/2;k++)
       {
       h[k]=2;
       }
      }
    }

   for(i=0;i<N;i++)
    {
    idft__XK_re[i]= dft__Xk_re[i]*h[i];
    idft__XK_im[i]= dft__Xk_im[i]*h[i];
    }

  
   //idft__XN= idft(idft__XK, n);
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    idft__XN.re[n]=0;
    idft__XN.im[n]=0;
    for(k=0;k<N;k++)
     {
     //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
     angle=2*M_PI*k*n/N;
     idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
     idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
     }
   idft__XN.im[n]=idft__XN.im[n]/N;
   idft__XN.re[n]=idft__XN.re[n]/N;

   }

  
  return idft__XN;
   }

  CData getanalyticsignal(double *xn, int num)
   {
   CData x_n,Hilb;
   int i;
     //x=x+jy,y=0
     //hilb(x)=a+jb
     //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
     //I=x(t), Q=j*Hilb(xre(t))

  Hilb=hilbert(xn,num);
   
   for (i=0; i<num; i++)
   {
    x_n.re[i]=xn[i];
    x_n.im[i]=Hilb.im[i];
   }
   
   return x_n;
   }

  
   
  double Arg(double re, double im)
  {
   
  return atan2 (im,re)  ; 
  }

  CData  get_mod_fi(CData x, int n, int mode   )
   {
    int i;
    CData res;
    double fi;
    for (i=0;i<n;i++)
    {
    res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
    //double treshh=1e9;
    // x.im[i]=x.im[i]*treshh;
    // x.im[i]=floor(x.im[i])/treshh;
    // x.re[i]=x.re[i]*treshh;
    // x.re[i]=floor(x.re[i])/treshh;
    res.im[i]=Arg(x.re[i],x.im[i]);
    if (mode==1)  { res.im[i]=res.im[i]*180/M_PI;  }
            
    }
    return  res;
   }

  CData GetSpectrum(double *xn, int num)
  {
   
  CData xn_iq, Xf_IQ, Xf_modfi;
  xn_iq=getanalyticsignal( xn,  num);
   
  Xf_IQ=cdft(xn_iq  ,num );
  Xf_IQ=ScaleDFT(Xf_IQ,num);
  Xf_modfi=get_mod_fi(Xf_IQ, num,1   );
   
  return Xf_modfi;
  }

  Array getenvelope(double *x, int num)
   {
   int i;
   CData  y;
   Array env;
   y=hilbert(x,num);
   for (i=0;i<num;i++)
    {
    env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
    }
    return env;
   }

  FILE   *fp1,*fp2 ;

  
  int _tmain(int argc, _TCHAR* argv[])
  {

  double xin[NMAX];
  CData  Xf;

  char key;
  int k,n,i,N1;
  double fs,T,x,f;

  
  printf ("\n For input  x(t) from file input '1' ");
  printf ("\n For manual  input  x(t) input   '0' ");
  scanf("%s",&key);
  if (key=='1' ) { goto label_05;}
  if (key=='0' ) { goto label_06;}

  label_05:

   printf("\n Open xinput.txt for reading x[n]   ");
  if ((fp1=fopen("xinput.txt","r"))==NULL)  printf ("\n File opening error ");
  fscanf(fp1,"%le",&fs);
  fscanf(fp1,"%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   fscanf (fp1,"%d",&k);
   fscanf (fp1,"%le",&xin[n]);
  }

  
  fclose(fp1);

  goto label_07;

  label_06:

  printf("\n Input fs,Hz , ( %lf ) ",fs) ;
  scanf( "%le",&fs);
  printf("\n Input number of samples of xin[n]  ") ;
  scanf( "%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   T=n*1/(N1*fs); 
   printf ("\n n=%d ; t=%le s; xin[n]= ",n,T);
   scanf ( "%le",&xin[n]);
  }

   label_07:
  printf( "\n fs=%le Hz ",fs);
  printf( "\n Nsamples= %d ",N1);

  for (n=0 ; n<N1;n++)
  {
   T=n*1/(N1*fs); 
   printf("\n n=%d x[n]=%le,  t=%le s",n,xin[n],T);
  }

  printf("\n Input 0 for continue ");
  scanf("%s",&key);
  printf("\n OK.   ");

  
  /*****************************/

  label1:
  //Xf.re= |Xf|,  Xf.im -arg(Xf) (only in this program)
  printf("\n Obtaining spectrum of x[n]");
  Xf=GetSpectrum(xin, N1);

  printf("\n Spectrum of x[n] : ");

   if ((fp2=fopen("spectrum.txt","a"))==NULL)  printf ("\n File opening error");
   fprintf(fp2,"\n   %le   %d \n",fs,N1);
   //fprintf(fp2, "\n n;\t f, Hz;\t |Xf(jw)| ;\t  arg(Xf), deg ");
  for (i=0; i<N1;i++)
  {
  f=fs*i/N1;
  //cout <<"\n f="<< f<<" Hz;|Xf(jw)|="<<Xf.re[i]<<";arg(Xf)="<<Xf.im[i]<<" deg ";
   printf ("\n f=%le Hz;|Xf(jw)|=%le;arg(Xf)=%lf deg ",  f,Xf.re[i], Xf.im[i]);
   fprintf(fp2,"\n %d   %le  %le",i ,Xf.re[i],Xf.im[i]);

  }
  printf ("\n  input   '0' for exit ");
  char key1;
  cin>>key1;

  return 0;

  }

  // stdafx.cpp : source file that includes just the standard includes
  // getspec.pch will be the pre-compiled header
  // stdafx.obj will contain the pre-compiled type information

  #include "stdafx.h"

  // TODO: reference any additional headers you need in STDAFX.H
  // and not in this file

  // stdafx.h : include file for standard system include files,
  // or project specific include files that are used frequently, but
  // are changed infrequently
  //

  #pragma once

  
  #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
  #include <stdio.h> //edit
  #include <tchar.h> //edit

  // TODO: reference additional headers your program requires here

  xinput.txt

     4.410000e+004   16

   0  1.000000e+001  
   1  7.071068e+000  
   2  6.123032e-016  
   3  -7.071068e+000  
   4  -1.000000e+001  
   5  -7.071068e+000  
   6  -1.836910e-015  
   7  7.071068e+000  
   8  1.000000e+001  
   9  7.071068e+000  
   10  3.061516e-015  
   11  -7.071068e+000  
   12  -1.000000e+001  
   13  -7.071068e+000  
   14  -4.286122e-015  
   15  7.071068e+000   

  spectrum.txt

     4.410000e+004   16

   0   7.216450e-016  -1.126199e+002
   1   7.237769e-016  -8.560129e+001
   2   1.000000e+001  9.223608e-015
   3   9.436896e-016  1.800000e+002
   4   3.373574e-016  3.464047e+001
   5   5.950643e-015  -5.660642e+001
   6   1.330312e-007  1.800000e+002
   7   7.196406e-015  -2.652734e+000
   8   3.760637e-015  -7.363206e+001
   9   3.624946e-015  -1.272209e+002
   10   2.820163e-015  1.438068e+002
   11   5.525033e-015  -1.015922e+002
   12   4.814310e-015  7.423349e+000
   13   1.110917e-014  6.612472e+001
   14   6.445273e-015  -9.246812e+001
   15   8.760391e-015  -1.326889e+002

  
  
  

  • Edited by USERPC01 Tuesday, January 23, 2018 10:29 PM adding xinput.txt

  Tuesday, January 23, 2018 10:14 PM
• 

  

  0

  Sign in to vote

  Program for signal synth. from freq. data

  // synthsig.cpp : Defines the entry point for the console application.
  //

  #include "stdafx.h"
  #include <stdio.h>
  #include <stdlib.h>
  #include <math.h>
  #include <iostream>

  #define NMAX 1000
  #define M_PI 3.1415926535897932384626433832795

  typedef struct {
  double re[NMAX];
  double im[NMAX];
  } CData;

   typedef struct  {
       double data[NMAX];
   } Array;

  
  CData  cdft(CData xn  ,int N )
   {
     
    
   CData Xk ;
   int k,n  ;
   double angle;

   for (k=0;k<N;k++)
    {
    Xk.re[k]=0;
    Xk.im[k]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*M_PI*k*n/N;  
      Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
      Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
      }

   }

   return Xk;
   }

   CData  cidft(CData idft_Xk, int N  )
   {

   int k,n ;
   double angle;
   CData  idft_xn;
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    idft_xn.re[n]=0;
    idft_xn.im[n]=0;
    for(k=0;k<N;k++)
     {
   
     angle=2*M_PI*k*n/N;
     idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
     idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
     }

   idft_xn.im[n]=idft_xn.im[n]/N;
   idft_xn.re[n]=idft_xn.re[n]/N;

   }

  return idft_xn ;
   }

  CData ScaleDFT(CData xk,int  num)
   {

   int i;
   for (i=0; i<num ; i++)
     {
      xk.re[i]=xk.re[i]/num;
      xk.im[i]=xk.im[i]/num;
     }
   return xk;
   }

  CData ScaleIQidft(CData Xk, int num)
   {
   int i;
     for (i=0;  i<num ; i++)
    {
     Xk.re[i]=Xk.re[i]*num;
     Xk.im[i]=Xk.im[i]*num;
    }

  return Xk;

  }

   CData hilbert(double *xr , int N)
   {
    int  k,n,i, no2 ;
    double angle;
    double h[NMAX];
    double dft__xn_re[NMAX];
    double dft__xn_im[NMAX];
     double dft__Xk_re[NMAX];
     double dft__Xk_im[NMAX];
       double idft__XK_re[NMAX];
       double idft__XK_im[NMAX];
     CData idft__XN;

   
    if (N==0)
    {
    printf("\n Error");
    idft__XN.re[0]=0;
    idft__XN.im[0]=0;
    return idft__XN;
    }
    //xnreal=real(xr);
    no2=(int)floor ((double)N/2);
   
    for (i=0;i<N;i++)
    {
     dft__xn_re[i]=xr[i];
     dft__xn_im[i]=0;
    }

   //fft
   // dft__Xk=dft(dft__xn,n);

  for (k=0;k<N;k++)
    {
    dft__Xk_re[k]=0;
    dft__Xk_im[k]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*M_PI*k*n/N;
      dft__Xk_re[k]+= dft__xn_re[n]*cos(angle)- dft__xn_im[n]*sin(angle);
      dft__Xk_im[k]+= dft__xn_re[n]*sin(angle)+ dft__xn_im[n]*cos(angle);
      }

   }

   //X={Xk.re,Xk.im};

   for (i=0;i<=N;i++)
    {
    h[i]=0;
    }

  
    if((N>0)&&((2*floor((double) N/2))==(N))) // n is even
    {

    h[0]=1;
     h[N/2]=1;
     for(k=1;k<N/2;k++)
     {
     h[k]=2;
     }
    }
    else //n is odd
    {
      if(N>0)
      {
       h[0]=1;
       for(k=1;k<(N+1)/2;k++)
       {
       h[k]=2;
       }
      }
    }

   for(i=0;i<N;i++)
    {
    idft__XK_re[i]= dft__Xk_re[i]*h[i];
    idft__XK_im[i]= dft__Xk_im[i]*h[i];
    }

  
   //idft__XN= idft(idft__XK, n);
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    idft__XN.re[n]=0;
    idft__XN.im[n]=0;
    for(k=0;k<N;k++)
     {
     //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
     angle=2*M_PI*k*n/N;
     idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
     idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
     }
   idft__XN.im[n]=idft__XN.im[n]/N;
   idft__XN.re[n]=idft__XN.re[n]/N;

   }

  
  return idft__XN;
   }

  CData getanalyticsignal(double *xn, int num)
   {
   CData x_n,Hilb;
   int i;
     //x=x+jy,y=0
     //hilb(x)=a+jb
     //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
     //I=x(t), Q=j*Hilb(xre(t))

  Hilb=hilbert(xn,num);
   
   for (i=0; i<num; i++)
   {
    x_n.re[i]=xn[i];
    x_n.im[i]=Hilb.im[i];
   }
   
   return x_n;
   }

  
   CData GetIQfromR_fi(CData Xinp, int numb)
   {
   int i;
   CData res;
    for (i=0;i<numb;i++)
    {
     Xinp.im[i]=Xinp.im[i]*M_PI/180;
    res.re[i]= Xinp.re[i]*cos(Xinp.im[i]);
    res.im[i]= Xinp.re[i]*sin(Xinp.im[i]);
    }

   return res;
   }

  double Arg(double re, double im)
  {
   
  return atan2(im,re); 
  }

  CData  get_mod_fi(CData x, int n, int mode   )
   {
    int i;
    CData res;
    double fi;
    for (i=0;i<n;i++)
    {
    res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
    //double treshh=1e9;
    // x.im[i]=x.im[i]*treshh;
    // x.im[i]=floor(x.im[i])/treshh;
    // x.re[i]=x.re[i]*treshh;
    // x.re[i]=floor(x.re[i])/treshh;
    res.im[i]=Arg(x.re[i],x.im[i]);
    if (mode==1)  { res.im[i]=res.im[i]*180/M_PI;  }
            
    }
    return  res;
   }

  CData GetSpectrum(double *xn, int num)
  {
   
  CData xn_iq, Xf_IQ, Xf_modfi;
  xn_iq=getanalyticsignal( xn,  num);
   
  Xf_IQ=cdft(xn_iq  ,num );
  Xf_IQ=ScaleDFT(Xf_IQ,num);
  Xf_modfi=get_mod_fi(Xf_IQ, num,1   );
   
  return Xf_modfi;
  }

  CData SynthSignal(CData Xf_modfi, int num)

  {
  CData Xf_IQ,xn_iq;
  Xf_IQ=GetIQfromR_fi( Xf_modfi, num);
  xn_iq=cidft(Xf_IQ,num  );
  xn_iq=ScaleIQidft(xn_iq,  num);
  return xn_iq;
  }

  Array getenvelope(double *x, int num)
   {
   int i;
   CData  y;
   Array env;
   y=hilbert(x,num);
   for (i=0;i<num;i++)
    {
    env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
    }
    return env;
   }

  FILE   *fp1,*fp2 ;

  
  //int main ()
  int _tmain(int argc, _TCHAR* argv[])
  {

  double xin[NMAX];
  CData  Xf;
  CData xt;
  char key;
  int k,n,i,N1;
  double fs=44100, T,x,f;

  
  printf ("\n For input  x(t) from file input '1' ");
  printf ("\n For manual  input  x(t) input   '0' ");
  scanf("%s",&key);
  if (key=='1' ) { goto label_05;}
  if (key=='0' ) { goto label_06;}

  label_05:

   printf("\n Open spectrum.txt for reading Xf   ");
  if ((fp1=fopen("spectrum.txt","r"))==NULL)  printf ("\n File opening error ");
  fscanf(fp1,"%le",&fs);
  fscanf(fp1,"%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   fscanf (fp1,"%d",&k);
   fscanf (fp1,"%le", &Xf.re[n]);
   fscanf (fp1,"%le", &Xf.im[n]);
  }

  
  fclose(fp1);

  goto label_07;

  label_06:

  printf("\n Input fs,Hz , ( %lf ) ", fs) ;
  scanf( "%le",&fs);
  printf("\n Input number of samples of xin[n]  ") ;
  scanf( "%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   
   f=fs*n/N1; 
   printf ("\n n=%d ; f=%le Hz; ",n,f);
   printf ("\n Input |X(jw)| :" );
   scanf ( "%le",&Xf.re[n]);
   printf ("\nInput arg(X(jw)), deg :" );
   scanf ( "%le",&Xf.im[n]);
  }

   label_07:
  printf( "\n fs=%le Hz ",fs);
  printf( "\n Nsamples= %d ",N1);

  for (n=0 ; n<N1;n++)
  {
    f=fs*n/N1;
   printf ("\n f=%le Hz;|Xf(jw)|=%le;arg(Xf)=%le deg", f,Xf.re[n],Xf.im[n]);
  }

  printf("\n Input 0 for continue ");
  scanf("%s",&key);
  printf("\n OK.   ");

  
  /*****************************/

  label1:
   
  printf("\n Parsing   Xf->x[n]");
  xt=SynthSignal( Xf , N1);

  printf("\n Signal (xinput.txt): ");
   
   
  if ((fp2=fopen("xinput.txt","a"))==NULL)  printf ("\n File opening error");
  fprintf(fp2,"\n   %le   %d \n",fs,N1);

  for (i=0; i<N1;i++)
  {
   
  printf ("\n  %d ;   %le ;   %le ;",  i,xt.re[i],xt.im[i]);
  //fprintf(fp2,"\n %d  %le   %le ", i,xt.re[i], xt.im[i]);
  fprintf(fp2,"\n %d  %le   ", i,xt.re[i] );
  }
  printf ("\n   input   '0' for exit ");
  std::cin>>key;
  return 0;

  }

  // stdafx.cpp : source file that includes just the standard includes
  // synthsig.pch will be the pre-compiled header
  // stdafx.obj will contain the pre-compiled type information

  #include "stdafx.h"

  // TODO: reference any additional headers you need in STDAFX.H
  // and not in this file

  // stdafx.h : include file for standard system include files,
  // or project specific include files that are used frequently, but
  // are changed infrequently
  //

  #pragma once

  
  #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
  #include <stdio.h>
  #include <tchar.h>

  // TODO: reference additional headers your program requires here

  spectrum.txt

    44100    16

   0   0   0
   1   0   0
   2   10   0
   3   0   0
   4   0   0
   5   0   0
   6   0   0
   7   0   0
   8   0   0
   9   0   0
   10   0   0
   11   0   0
   12   0   0
   13   0   0
   14  0   0
   15  0   0

  Friday, January 26, 2018 1:27 PM
• 

  

  0

  Sign in to vote

  // orthosynth.cpp : Defines the entry point for the console application.
  //

  #include "stdafx.h"

  
  #include <stdio.h>
  #include <stdlib.h>
  #include <math.h>
  #include <iostream>

  #define NMAX 1000
  #define M_PI 3.1415926535897932384626433832795

  typedef struct {
  double re[NMAX];
  double im[NMAX];
  } CData;

   typedef struct  {
       double data[NMAX];
   } Array;

  
  CData  cdft(CData xn  ,int N )
   {
     
   
   CData Xk ;
   int k,n  ;
   double angle;

   for (k=0;k<N;k++)
    {
    Xk.re[k]=0;
    Xk.im[k]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*M_PI*k*n/N;  
      Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
      Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
      }

   }

   return Xk;
   }

   CData  cidft(CData idft_Xk, int N  )
   {

   int k,n ;
   double angle;
   CData  idft_xn;
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    idft_xn.re[n]=0;
    idft_xn.im[n]=0;
    for(k=0;k<N;k++)
     {
   
     angle=2*M_PI*k*n/N;
     idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
     idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
     }

   idft_xn.im[n]=idft_xn.im[n]/N;
   idft_xn.re[n]=idft_xn.re[n]/N;

   }

  return idft_xn ;
   }

  CData ScaleDFT(CData xk,int  num)
   {

   int i;
   for (i=0; i<num ; i++)
     {
      xk.re[i]=xk.re[i]/num;
      xk.im[i]=xk.im[i]/num;
     }
   return xk;
   }

  CData ScaleIQidft(CData Xk, int num)
   {
   int i;
     for (i=0;  i<num ; i++)
    {
     Xk.re[i]=Xk.re[i]*num;
     Xk.im[i]=Xk.im[i]*num;
    }

  return Xk;

  }

   CData hilbert(double *xr , int N)
   {
    int  k,n,i, no2 ;
    double angle;
    double h[NMAX];
    double dft__xn_re[NMAX];
    double dft__xn_im[NMAX];
     double dft__Xk_re[NMAX];
     double dft__Xk_im[NMAX];
       double idft__XK_re[NMAX];
       double idft__XK_im[NMAX];
     CData idft__XN;

   
    if (N==0)
    {
    printf("\n Error");
    idft__XN.re[0]=0;
    idft__XN.im[0]=0;
    return idft__XN;
    }
    //xnreal=real(xr);
    no2=(int)floor ((double)N/2);
   
    for (i=0;i<N;i++)
    {
     dft__xn_re[i]=xr[i];
     dft__xn_im[i]=0;
    }

   //fft
   // dft__Xk=dft(dft__xn,n);

  for (k=0;k<N;k++)
    {
    dft__Xk_re[k]=0;
    dft__Xk_im[k]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*M_PI*k*n/N;
      dft__Xk_re[k]+= dft__xn_re[n]*cos(angle)- dft__xn_im[n]*sin(angle);
      dft__Xk_im[k]+= dft__xn_re[n]*sin(angle)+ dft__xn_im[n]*cos(angle);
      }

   }

   //X={Xk.re,Xk.im};

   for (i=0;i<=N;i++)
    {
    h[i]=0;
    }

  
    if((N>0)&&((2*floor((double)N/2))==(N))) // n is even
    {

    h[0]=1;
     h[N/2]=1;
     for(k=1;k<N/2;k++)
     {
     h[k]=2;
     }
    }
    else //n is odd
    {
      if(N>0)
      {
       h[0]=1;
       for(k=1;k<(N+1)/2;k++)
       {
       h[k]=2;
       }
      }
    }

   for(i=0;i<N;i++)
    {
    idft__XK_re[i]= dft__Xk_re[i]*h[i];
    idft__XK_im[i]= dft__Xk_im[i]*h[i];
    }

  
   //idft__XN= idft(idft__XK, n);
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    idft__XN.re[n]=0;
    idft__XN.im[n]=0;
    for(k=0;k<N;k++)
     {
     //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
     angle=2*M_PI*k*n/N;
     idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
     idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
     }
   idft__XN.im[n]=idft__XN.im[n]/N;
   idft__XN.re[n]=idft__XN.re[n]/N;

   }

  
  return idft__XN;
   }

  CData getanalyticsignal(double *xn, int num)
   {
   CData x_n,Hilb;
   int i;
     //x=x+jy,y=0
     //hilb(x)=a+jb
     //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
     //I=x(t), Q=j*Hilb(xre(t))

  Hilb=hilbert(xn,num);
   
   for (i=0; i<num; i++)
   {
    x_n.re[i]=xn[i];
    x_n.im[i]=Hilb.im[i];
   }
   
   return x_n;
   }

  
   CData GetIQfromR_fi(CData Xinp, int numb)
   {
   int i;
   CData res;
    for (i=0;i<numb;i++)
    {
     Xinp.im[i]=Xinp.im[i]*M_PI/180;
    res.re[i]= Xinp.re[i]*cos(Xinp.im[i]);
    res.im[i]= Xinp.re[i]*sin(Xinp.im[i]);
    }

   return res;
   }

  double Arg(double re, double im)
  {
   
  return atan2(im,re); 
  }

  CData  get_mod_fi(CData x, int n, int mode   )
   {
    int i;
    CData res;
    double fi;
    for (i=0;i<n;i++)
    {
    res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
    //double treshh=1e9;
    // x.im[i]=x.im[i]*treshh;
    // x.im[i]=floor(x.im[i])/treshh;
    // x.re[i]=x.re[i]*treshh;
    // x.re[i]=floor(x.re[i])/treshh;
    res.im[i]=Arg(x.re[i],x.im[i]);
    if (mode==1)  { res.im[i]=res.im[i]*180/M_PI;  }
            
    }
    return  res;
   }

  CData GetSpectrum(double *xn, int num)
  {
   
  CData xn_iq, Xf_IQ, Xf_modfi;
  xn_iq=getanalyticsignal( xn,  num);
   
  Xf_IQ=cdft(xn_iq  ,num );
  Xf_IQ=ScaleDFT(Xf_IQ,num);
  Xf_modfi=get_mod_fi(Xf_IQ, num,1   );
   
  return Xf_modfi;
  }

  CData SynthSignal(CData Xf_modfi, int num)

  {
  CData Xf_IQ,xn_iq;
  Xf_IQ=GetIQfromR_fi( Xf_modfi, num);
  xn_iq=cidft(Xf_IQ,num  );
  xn_iq=ScaleIQidft(xn_iq,  num);
  return xn_iq;
  }

  Array getenvelope(double *x, int num)
   {
   int i;
   CData  y;
   Array env;
   y=hilbert(x,num);
   for (i=0;i<num;i++)
    {
    env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
    }
    return env;
   }

  FILE   *fp1,*fp2 ;

  //int _tmain(int argc, _TCHAR* argv[])
   int main ()
  {

  double xin[NMAX];
  CData  Xf;
  CData xt;

  int k,n,i,N1;
  double fs=44100,T,x,f;
  int  key=0;

  printf ("\n For input  x(t) from file input '1' ");
  printf ("\n For manual  input  x(t) input   '0' ");
  std::cin>>key;
  switch(key)
  {
  case  1 :
       goto label_05;
    break;
  case  0 :
       goto label_06;
    break;
    
  }

  label_05:

   printf("\n Open spectrum.txt for reading Xf   ");
  if ((fp1=fopen("spectrum.txt","r"))==NULL)  printf ("\n File opening error ");
  fscanf(fp1,"%le",&fs);
  fscanf(fp1,"%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   fscanf (fp1,"%d",&k);
   fscanf (fp1,"%le", &Xf.re[n]);
   fscanf (fp1,"%le", &Xf.im[n]);
  }

  
  fclose(fp1);

  goto label_07;

  label_06:

  printf("\n Input fs,Hz , ( %lf ) ",fs) ;
  scanf( "%le",&fs);
  printf("\n Input number of samples of xin[n]  ") ;
  scanf( "%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   
   f=fs*n/N1; 
   printf ("\n n=%d ; f=%le Hz; ",n,f);
   printf ("\n Input |X(jw)| :" );
   scanf ( "%le",&Xf.re[n]);
   printf ("\nInput arg(X(jw)), deg :" );
   scanf ( "%le",&Xf.im[n]);
  }

   label_07:
  printf( "\n fs=%le Hz ",fs);
  printf( "\n Nsamples= %d ",N1);

  for (n=0 ; n<N1;n++)
  {
    f=fs*n/N1;
   printf ("\n f=%le Hz;|Xf(jw)|=%le;arg(Xf)=%le deg", f,Xf.re[n],Xf.im[n]);
  }

  printf("\n Input 0 for continue ");
  scanf("%s",&key);
  printf("\n OK.   ");

  
  /*****************************/

  label1:
   
  printf("\n Parsing   Xf->x[n]");
  xt=SynthSignal( Xf , N1);

  printf("\n Signal (xinput.txt): ");
   
   
  if ((fp2=fopen("orthoxn.txt","a"))==NULL)  printf ("\n File opening error");
  fprintf(fp2,"\n   %le   %d \n",fs,N1);

  for (i=0; i<N1;i++)
  {
   
  printf ("\n  %d ;   %le ;   %le  ",  i,xt.re[i],xt.im[i]);
   fprintf(fp2,"\n %d  %le   %le ", i,xt.re[i], xt.im[i]);
  //fprintf(fp2,"\n %d  %le   ", i,xt.re[i] );
  }
  printf ("\n   input   '0' for exit ");
  char key1;
  std::cin>>key1;
  std::cout << key1;

  
  return 0;

  }

  // stdafx.cpp : source file that includes just the standard includes
  // orthosynth.pch will be the pre-compiled header
  // stdafx.obj will contain the pre-compiled type information

  #include "stdafx.h"

  // TODO: reference any additional headers you need in STDAFX.H
  // and not in this file

  // stdafx.h : include file for standard system include files,
  // or project specific include files that are used frequently, but
  // are changed infrequently
  //

  #pragma once

  
  #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
  #include <stdio.h>
  #include <tchar.h>

  // TODO: reference additional headers your program requires here

  Friday, January 26, 2018 1:52 PM
• 

  

  0

  Sign in to vote

  // spectrumratio.cpp : Defines the entry point for the console application.
  //

  #include "stdafx.h"

   

  #include <stdio.h>
  #include <stdlib.h>
  #include <math.h>
  #include <iostream>

  #define NMAX 1000
  #define M_PI 3.1415926535897932384626433832795

  typedef struct {
  double re[NMAX];
  double im[NMAX];
  } CData;

   typedef struct  {
       double data[NMAX];
   } Array;

  
  CData  cdft(CData xn  ,int N )
   {
     
   
   CData Xk ;
   int k,n  ;
   double angle;

   for (k=0;k<N;k++)
    {
    Xk.re[k]=0;
    Xk.im[k]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*M_PI*k*n/N;  
      Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
      Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
      }

   }

   return Xk;
   }

   CData  cidft(CData idft_Xk, int N  )
   {

   int k,n ;
   double angle;
   CData  idft_xn;
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    idft_xn.re[n]=0;
    idft_xn.im[n]=0;
    for(k=0;k<N;k++)
     {
   
     angle=2*M_PI*k*n/N;
     idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
     idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
     }

   idft_xn.im[n]=idft_xn.im[n]/N;
   idft_xn.re[n]=idft_xn.re[n]/N;

   }

  return idft_xn ;
   }

  CData ScaleDFT(CData xk,int  num)
   {

   int i;
   for (i=0; i<num ; i++)
     {
      xk.re[i]=xk.re[i]/num;
      xk.im[i]=xk.im[i]/num;
     }
   return xk;
   }

  CData ScaleIQidft(CData Xk, int num)
   {
   int i;
     for (i=0;  i<num ; i++)
    {
     Xk.re[i]=Xk.re[i]*num;
     Xk.im[i]=Xk.im[i]*num;
    }

  return Xk;

  }

   CData hilbert(double *xr , int N)
   {
    int  k,n,i, no2 ;
    double angle;
    double h[NMAX];
    double dft__xn_re[NMAX];
    double dft__xn_im[NMAX];
     double dft__Xk_re[NMAX];
     double dft__Xk_im[NMAX];
       double idft__XK_re[NMAX];
       double idft__XK_im[NMAX];
     CData idft__XN;

   
    if (N==0)
    {
    printf("\n Error");
    idft__XN.re[0]=0;
    idft__XN.im[0]=0;
    return idft__XN;
    }
    //xnreal=real(xr);
    no2=(int)floor ((double)N/2);
   
    for (i=0;i<N;i++)
    {
     dft__xn_re[i]=xr[i];
     dft__xn_im[i]=0;
    }

   //fft
   // dft__Xk=dft(dft__xn,n);

  for (k=0;k<N;k++)
    {
    dft__Xk_re[k]=0;
    dft__Xk_im[k]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*M_PI*k*n/N;
      dft__Xk_re[k]+= dft__xn_re[n]*cos(angle)- dft__xn_im[n]*sin(angle);
      dft__Xk_im[k]+= dft__xn_re[n]*sin(angle)+ dft__xn_im[n]*cos(angle);
      }

   }

   //X={Xk.re,Xk.im};

   for (i=0;i<=N;i++)
    {
    h[i]=0;
    }

  
    if((N>0)&&((2*floor((double)N/2))==(N))) // n is even
    {

    h[0]=1;
     h[N/2]=1;
     for(k=1;k<N/2;k++)
     {
     h[k]=2;
     }
    }
    else //n is odd
    {
      if(N>0)
      {
       h[0]=1;
       for(k=1;k<(N+1)/2;k++)
       {
       h[k]=2;
       }
      }
    }

   for(i=0;i<N;i++)
    {
    idft__XK_re[i]= dft__Xk_re[i]*h[i];
    idft__XK_im[i]= dft__Xk_im[i]*h[i];
    }

  
   //idft__XN= idft(idft__XK, n);
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    idft__XN.re[n]=0;
    idft__XN.im[n]=0;
    for(k=0;k<N;k++)
     {
     //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
     angle=2*M_PI*k*n/N;
     idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
     idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
     }
   idft__XN.im[n]=idft__XN.im[n]/N;
   idft__XN.re[n]=idft__XN.re[n]/N;

   }

  
  return idft__XN;
   }

  CData getanalyticsignal(double *xn, int num)
   {
   CData x_n,Hilb;
   int i;
     //x=x+jy,y=0
     //hilb(x)=a+jb
     //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
     //I=x(t), Q=j*Hilb(xre(t))

  Hilb=hilbert(xn,num);
   
   for (i=0; i<num; i++)
   {
    x_n.re[i]=xn[i];
    x_n.im[i]=Hilb.im[i];
   }
   
   return x_n;
   }

  
   
  double Arg(double re, double im)
  {
   
  return atan2 (im,re)  ; 
  }

  CData  get_mod_fi(CData x, int n, int mode   )
   {
    int i;
    CData res;
    double fi;
    for (i=0;i<n;i++)
    {
    res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
    //double treshh=1e9;
    // x.im[i]=x.im[i]*treshh;
    // x.im[i]=floor(x.im[i])/treshh;
    // x.re[i]=x.re[i]*treshh;
    // x.re[i]=floor(x.re[i])/treshh;
    res.im[i]=Arg(x.re[i],x.im[i]);
    if (mode==1)  { res.im[i]=res.im[i]*180/M_PI;  }
            
    }
    return  res;
   }

  CData GetSpectrum(double *xn, int num)
  {
   
  CData xn_iq, Xf_IQ, Xf_modfi;
  xn_iq=getanalyticsignal( xn,  num);
   
  Xf_IQ=cdft(xn_iq  ,num );
  Xf_IQ=ScaleDFT(Xf_IQ,num);
  Xf_modfi=get_mod_fi(Xf_IQ, num,1   );
   
  return Xf_modfi;
  }

  Array getenvelope(double *x, int num)
   {
   int i;
   CData  y;
   Array env;
   y=hilbert(x,num);
   for (i=0;i<num;i++)
    {
    env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
    }
    return env;
   }
   
   
   
  double getmax( Array   k , int N)
  {
  double kmax=k.data[0];
  double k1=0;
  int i;
  for (i=1;i<N;i++)
  {
   k1=k.data[i];
    if(fabs( k1 )>fabs(kmax)) { kmax= k.data[i]; } 
    
  }
   
   return kmax;
  }
   
  double getmin( Array  k, int N)
  {
  double kmin=k.data[0];
  int i;
  for (i=0;i<N;i++)
  {
    if(fabs(k.data[i])<fabs(kmin)) { kmin= k.data[i]; } 
  }
   return kmin;
  }

  double getmid( Array  k,int N)
  {
  double kmid=0;
  int i;
  for (i=0;i<N;i++)
  {
   kmid+=   k.data[i];
  }
  kmid=kmid/N;
   return kmid;
  }

  
  FILE   *fp1,*fp2 ;

  int _tmain(int argc, _TCHAR* argv[])
  //int main ()
  {
  double delta =1e-15;
  double xin[NMAX];
  double yout[NMAX];
  CData Xf1,Xf2,Xf3;
  double Kmax=0, Kmin=0, Kmid=0;
  int key;
  int k,n,i,N1,N2;
  double fs1,fs2,T,x,f;

  
  /***********input x[n] *************/
  label_00:
  printf ("\n For input  x[n] from file input '1' ");
  printf ("\n For manual  input  x[n] input   '0' ");
   std::cin>>key;
  if (key== 1  ) { goto label_01;}
  if (key== 0  ) { goto label_02;}

  label_01:

   printf("\n Open xinput.txt for reading x[n]   ");
  if ((fp1=fopen("xinput.txt","r"))==NULL)  printf ("\n File opening error ");
  fscanf(fp1,"%le",&fs1);
  fscanf(fp1,"%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   fscanf (fp1,"%d",&k);
   fscanf (fp1,"%le",&xin[n]);
  }

  
  fclose(fp1);

  goto label_03;

  label_02:

  printf("\n Input fs input,Hz  " ) ;
  scanf( "%le",&fs1);
  printf("\n Input number of samples of xin[n]  ") ;
  scanf( "%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   T=n*1/(N1*fs1); 
   printf ("\n n=%d ; t=%le s; xin[n]= ",n,T);
   scanf ( "%le",&xin[n]);
  }

   label_03:
  printf( "\n fs1=%le Hz ",fs1);
  printf( "\n Nsamples1= %d ",N1);

  for (n=0 ; n<N1;n++)
  {
   T=n*1/(N1*fs1); 
   printf("\n n=%d x[n]=%le,  t=%le s",n,xin[n],T);
  }

  printf("\n Input 0 for continue ");
  char key1;
  std::cin>>key1;
  printf("\n OK.   ");

  /***********input y[n] *************/

  printf ("\n For input  y[n] from file input '1' ");
  printf ("\n For manual  input  y[n] input   '0' ");
  std::cin>>key;
  if (key== 1  ) { goto label_04;}
  if (key== 0  ) { goto label_05;}

  label_04:

   printf("\n Open youtput.txt for reading y[n]   ");
  if ((fp1=fopen("youtput.txt","r"))==NULL)  printf ("\n File opening error ");
  fscanf(fp1,"%le",&fs2);
  fscanf(fp1,"%d", &N2);

  for (n=0 ; n<N2;n++)
  {
   fscanf (fp1,"%d",&k);
   fscanf (fp1,"%le",&yout[n]);
  }

  
  fclose(fp1);

  goto label_06;

  label_05:

  printf("\n Input fs out,Hz , ( %lf ) ",fs1) ;
  scanf( "%le",&fs2);
  printf("\n Input number of samples of yout[n]  ") ;
  scanf( "%d",&N2);

  for (n=0 ; n<N2;n++)
  {
   T=n*1/(N2*fs2); 
   printf ("\n n=%d ; t=%le s; yout[n]= ",n,T);
   scanf ( "%le",&yout[n]);
  }

   label_06:

  
  if (N1!=N2) printf("Nin <> Nout");

  printf( "\n fs =%le Hz ",fs2);
  printf( "\n Nsamples  out = %d ",N2);

  for (n=0 ; n<N2;n++)
  {
   T=n*1/(N2*fs2); 
   printf("\n n=%d y[n]=%le,  t=%le s",n,yout[n],T);
  }

  printf("\n Input 0 for continue ");
  std::cin>>key;
  printf("\n OK.   ");

  /***********************/
  if ((fs2!=fs1)||(N1!=N2)) { printf("Nin <> Nout"); goto label_00;}

  
  /*****************************/

  label1:
   
  printf("\n Obtaining spectrum of x[n]");
  Xf1=GetSpectrum(xin, N1) ;
  printf("\n Obtaining spectrum of y[n]");
  Xf2=GetSpectrum(yout, N2) ;

  
  printf("\n  |Sy(jw)|/|Sx(jw)| : ");
   
  if ((fp2=fopen("ratio.txt","a"))==NULL)  printf ("\n File opening error");
  fprintf(fp2,"\n\n   %le   %d \n",fs1,N1);

  
  for (i=0; i<N1;i++)
  {
    if (Xf1.re[i]==0 ) { Xf3.re[i]=Xf2.re[i]/(Xf1.re[i]+delta); }
    if (Xf1.re[i]!=0 ) { Xf3.re[i]=Xf2.re[i]/(Xf1.re[i] ); }
   Xf3.im[i]=Xf2.im[i]-Xf1.im[i];

   
  printf ("\n %d ;|Sy|/|Sx|=%le ;arg(Sy(jw))-arg(Sx(jw))=%le deg ",  i,Xf3.re[i],Xf3.im[i]);
  fprintf (fp2,"\n %d   %le  %le",  i,Xf3.re[i], Xf3.im[i]);

  }
  fclose(fp2);

   

  printf ("\n  input   '0' for exit ");
  char key2;
  std::cin>>key2;

  return 0;

  }

  // stdafx.cpp : source file that includes just the standard includes
  // spectrumratio.pch will be the pre-compiled header
  // stdafx.obj will contain the pre-compiled type information

  #include "stdafx.h"

  // TODO: reference any additional headers you need in STDAFX.H
  // and not in this file

  // stdafx.h : include file for standard system include files,
  // or project specific include files that are used frequently, but
  // are changed infrequently
  //

  #pragma once

  
  #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
  #include <stdio.h>
  #include <tchar.h>

  // TODO: reference additional headers your program requires here

  Friday, January 26, 2018 2:04 PM
• 

  

  0

  Sign in to vote

  make.bat

  set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"

   
   javac THDObtainer.java
   jar cvfe THDObtainer.jar THDObtainer   *.class
   jar uf  THDObtainer.jar THDObtainer.java
  pause 0

   

  set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"
   rem java  THDObtainer
   java -jar THDObtainer.jar

  pause 0

  import java.util.Scanner;
  import java.io.*;

  
     class  THDParser
   {

   

   
   
  static int  NMAX = 1024;

    public static double[] Getfassoc(double fs , int  Numb )
  {
     double[] fvect  =new double[Numb];
    
     for(int k = 0;k<Numb;k++)
    {
       fvect[k] = fs * k / Numb;
   
    }
   return fvect;
  }
   
   
  public static double[] Gettassoc(double  fs ,    int Numb )
   {
    double[] tvect=new double[Numb];
     for(int N = 0 ; N<Numb; N++)
      {
       tvect[N] = N/fs ;
      }
   return tvect;
  }

  
    public static double  Arg(double x_re  , double x_im  )
    {
    double Arg=0;
   if (x_re > 0)  { Arg = Math.atan(x_im / x_re); }
   if ((x_re == 0) &&  (x_im > 0)) { Arg = Math.PI / 2;}
   if ((x_re == 0) &&  (x_im < 0)) { Arg = Math.PI / 2 ; }
   if ((x_re < 0) &&  (x_im >= 0)) { Arg = Math.atan(x_im / x_re) + Math.PI; }
   if ((x_re < 0) &&  (x_im < 0))  { Arg = Math.atan(x_im / x_re) - Math.PI; }
   if ((x_re == 0) &&  (x_im == 0))  {
                                    Arg = 0;
                                 /* MsgBox (" Xre=0 and Xim=0, Invalid result"); */
                                  }
   
  return Arg;
  }

  
   /*
   '36.  - 4. + 9.6568542i  - 4. + 4.i  - 4. + 1.6568542i  - 4.  - 4. - 1.6568542i
   
    '         column 7 to 8
    
    ' - 4. - 4.i  - 4. - 9.6568542i
  */

    public static double[][]  GetIQfromR_fi(  int Num  ,double [][] Xk, int mode )
  {
     //long int i ;
      
        double arg1 ;
        double[][] Xk1=new double [Num][2];
      arg1 = 1;
      if( mode == 1) { arg1 = Math.PI / 180;}       
      for(int i = 0 ;i< Num;i++)
       {
        Xk1[i][0] =  Xk[i][0]*Math.cos(arg1*Xk[i][1]);
        Xk1[i][1] =  Xk[i][0]*Math.sin(arg1*Xk[i][1]);
       }
   return Xk;   
  }

   

  
    public static double[][] CDFT(int numdft ,double [][] dft_xn )
  {
    int i, k, N;
   
   double angle ;
   double[][] dft_Xk=new double [numdft ][2];
  for( k = 0 ;k<numdft ;k++)
   {
    dft_Xk[k][0] = 0;
    dft_Xk[k][1] = 0;
     for( N = 0;N<numdft;N++ )
     {
       angle = -2*Math.PI*k*N/numdft;
     dft_Xk[k][0]+=dft_xn[N][0]*Math.cos(angle)-dft_xn[N][1]*Math.sin(angle);
     dft_Xk[k][1]+=dft_xn[N][0]*Math.sin(angle)+dft_xn[N][1]*Math.cos(angle);
     }
    }
  return dft_Xk;
  }
   
  public static double[][]  ScaleDFT(int Num , double [][]  Xk )
  {
   
    for(int i = 0 ;i<Num;i++)
    {
       Xk[i][0] = Xk[i][0]/Num;
       Xk[i][1] = Xk[i][1]/Num;
    }
   return Xk;
  }

  public static double[][] ScaleIQidft(int Num  , double[][] ScXk)
  {
    
     for (int i = 0;i<Num;i++)
     {
      ScXk[i][0] = ScXk[i][0]*Num;
      ScXk[i][1] = ScXk[i][1]*Num;
     }

  return ScXk;
  }

  public static double[][] CIDFT(int numidft , double[][] idft_Xk)
  {
    double angle;
    double[][] idft_xn=new double[numidft][2] ;
    for(int n = 0 ; n<numidft; n++)
    {
   
     idft_xn[n][0] = 0;
     idft_xn[n][1] = 0;
     for (int k=0; k<numidft; k++)
     {
      angle =2*Math.PI*k*n/numidft;
      idft_xn[n][0]+=idft_Xk[k][0]*Math.cos(angle)-idft_Xk[k][1]*Math.sin(angle);
      idft_xn[n][1]+=idft_Xk[k][0]*Math.sin(angle)+idft_Xk[k][1]*Math.cos(angle);
     }
     
      idft_xn[n][0]= idft_xn[n][0]/numidft;
      idft_xn[n][1]= idft_xn[n][1]/numidft;
   }
   
   return idft_xn;
  }

  public  static double[][] Hilbert(double[] xre , int numb  )
  {
    /* MATLAB algorithm */
     int k, i, no2;
     double[] h= new double [NMAX+1];
   
   
     double[][] dft_xn   =new  double [NMAX][2];
     double[][] dft_Xk   =new  double [NMAX][2];
     double[][] idft_Xk  =new  double [NMAX][2];
     double[][] idft_xn  =new  double [NMAX][2];

    if (numb == 0)
     {
             idft_xn[0][0]=0;
             idft_xn[0][1]=0;
      
     /*MsgBox (" Hilbert() :   Parsing error ")*/
      }

  
     no2 = (int) Math.floor(numb / 2); //
    
     for(i = 0;i<numb; i++)
     {
      dft_xn[i][0]=xre[i]; 
      dft_xn[i][1]=0; 
     }

  
     dft_Xk=CDFT(numb, dft_xn ); 
     i=0;
     for (i=0;i <= numb ;i++){ h[i]=0; }

    /*if  N is even*/
   
      if( (numb> 0)&&((2*no2)==numb) )
      {
          h[0] = 1; // 0
          h[no2] = 1;// f Niquist
           for( k = 1; k<no2;k++)
           {
              h[k]= 2;
           }
      }

     /* if N is odd */

      else {
           if(numb > 0)
            {

              h[0]  = 1;
              for( k = 1; k<=((0.5*(numb + 1))-1); k++)
              {
               h[k] = 2;
              }
         
             }
            }

  
       for( i=0;i<numb;i++)
        {
         idft_Xk[i][0]=dft_Xk[i][0]*h[i];
         idft_Xk[i][1]=dft_Xk[i][1]*h[i];
        }
   
      
       idft_xn=CIDFT(numb, idft_Xk);
     
   return idft_xn;
   
  }

  
  public static double[][] GetAnalyticSignal(double [] xnre , int Numb )
  {
   //  'I=x(t)
     // 'Q=jHilb(t)
   //  'Hilb=Hilbert()
     
     double[][] Hilb= new double[Numb][2];
     double[][] x_iq= new double[Numb][2];

     Hilb= Hilbert(xnre, Numb );
     for(int i = 0; i<Numb;i++)
     {
      x_iq[i][0] = xnre[i] ;
      x_iq[i][1] = Hilb[i][1];
     }
  return  x_iq;
  }

  public static double[] getenvelope(double [] xnre , int Numb)
     {
    
     
     double[][] Hilb= new double[Numb][2];
     double[]  env= new double[Numb] ;
     Hilb= Hilbert(xnre, Numb );

     for(int i = 0 ; i<Numb; i++)
     {
     env[i]=Math.sqrt((Hilb[i][0]*Hilb[i][0])+(Hilb[i][1]*Hilb[i][1]));
     }
  return env;

  }

  
  public static double[][] GetSpectrumIQ(double[] xn , int Num  )
   {
    double[][] x_iq= new double[Num ][2];
    double[][] Xf_iq= new double[Num ][2];
    x_iq=GetAnalyticSignal(xn , Num  ); 
    Xf_iq=CDFT(Num, x_iq);
    Xf_iq=ScaleDFT(Num,Xf_iq);
  return Xf_iq;
  }

  
  public static double[][] GetR_fi(int  Numb  , double[][] Xf_iq , int mode )
  {

      double[][] Xf_rfi=new double [Numb][2];
      for (int i = 0 ;i<Numb;i++)
     {
      Xf_rfi[i][0] =Math.sqrt((Xf_iq[i][0]*Xf_iq[i][0])+(Xf_iq[i][1]*Xf_iq[i][1]) );
      Xf_rfi[i][1] = Arg(Xf_iq[i][0], Xf_iq[i][1]);
      if( mode == 1) { Xf_rfi[i][1] = Xf_rfi[i][1]*180/Math.PI; }
     }

  return Xf_rfi;
  }

  
  public static double[][] GetSpectrum(double[] xn , int Num ,  int mode  )
  {
            double[][] x_iq= new double[Num][2];
            double[][] Xf_iq= new double[Num][2];
            double[][] Xf_rfi=new double [Num][2];

         x_iq=GetAnalyticSignal(xn, Num );
         Xf_iq=CDFT(Num,x_iq);
         Xf_iq=ScaleDFT(Num, Xf_iq);
         Xf_rfi=GetR_fi(Num, Xf_iq, mode);
    return Xf_rfi;
  }
   
   
  public static double[][] GetSpectrumfromIQ(double[][] x_iq,int Num, int mode)
  {  
         double[][] Xf_iq= new double[Num][2];
         double[][] Xf_rfi=new double [Num][2];

        Xf_iq=CDFT(Num,x_iq);
        Xf_iq=ScaleDFT(Num, Xf_iq);
        Xf_rfi=GetR_fi(Num,Xf_rfi, mode);
    return Xf_rfi;
  }

  
  public static double[][] SynthSignalfromIQ(int Num  , double[][] Xk_iq)
  {
      double[][] xn_iq= new double[Num][2];
      
            xn_iq= CIDFT(Num, Xk_iq);
            xn_iq= ScaleIQidft(Num, xn_iq);
  return xn_iq;
  }

  
  public static double[][] SynthSignalfromRFi( double[][] Xk_rfi, int Num, int mode)
  {
  double[][] Xk_iq= new double [Num][2];
  Xk_iq=GetIQfromR_fi(Num, Xk_rfi, mode );
  return SynthSignalfromIQ(Num, Xk_iq);
  }

  public static void printcomplex(double[][] xn  , int num, int col1,int col2)
  { 
      
       for (int i=0 ;i<num;i++)
          {
           
            System.out.print("\n i= ");
            System.out.print(i);
            System.out.print(" ");
            System.out.print(String.format("%.08f",xn[i][0]));
            System.out.print(" ");
            System.out.print(String.format("%.08f",xn[i][1]));
      
          }
   System.out.print("\n ");
  return;
  }

  public static void printarray(double[] xn , int  Num  , int col1   )
  {
       for (int i = 0; i<Num;i++)
        {
          System.out.print("\n ");
          System.out.print(String.format("%.08f",xn[i]));
        
        }
   System.out.print("\n ");
  return;
  }

  public static double deltat(double[] argF , double f0 , double fs , int N )
   {
    double[] fv =new double[N+1];
    int  N0;
    double deltat;
     fv=Getfassoc(fs, N);
     N0 = 0;
   for(int i = 0;i<N-1;i++)
   {
   if ((fv[i]>= f0)&&(fv[i+1]<=f0)) { N0= i; }
   }
   deltat =-(argF[N0] -argF[0])/(2*Math.PI*fv[N0]-2*Math.PI*fv[0]);
  return deltat;
  }

  
   
  public static void writetofile(double[][] xnc  ,  String FilePath , int N ) throws IOException
  {  
   File outFile1 = new File (FilePath);
   FileWriter fWriter1 = new FileWriter (outFile1);
   PrintWriter pWriter1 = new PrintWriter (fWriter1);

      for(int i = 0 ;i< N;i++)
      {
           pWriter1.println("\n ");
           pWriter1.println("  ");
    pWriter1.print(String.format("%d",i) );
           pWriter1.print("   "); 
    pWriter1.print(String.format("%e",xnc[i][0] ));
    pWriter1.print("   ");
           pWriter1.print(String.format("%e",xnc[i][1] ));
    pWriter1.print("   ");
     }
    pWriter1.close();
   return;        
    
   }

   

  
  public static void TestHilbert( )
  {

   int i ;
   int N = 8;
   double  fs  ;
   double[] fv=new double[N ];
   double[] tv=new double[N ];
   double[][] Xf=new double[N ][2];

   double[][] Yf=new double[N ][2];
   double[][] Hilb=new double[N][2];
   double[][] Idft=new double[N ][2];
   double[][] SpX=new double[N][2];
   
   double[] env=new double[N];
   double[][] Xn=new double[N][2];
   double[] Xn1=new double[N];
   
   double[] x_re=new double[N];
   double[][] Yn=new double[N][2];
   
   
   Xf[0][0] = 1;
   Xf[1][0] = 2;
   Xf[2][0] = 3;
   Xf[3][0] = 4;
   Xf[4][0] = 5;
   Xf[5][0] = 6;
   Xf[6][0] = 7;
   Xf[7][0] = 8;

   Xf[0][1] = 0;
   Xf[1][1] = 0;
   Xf[2][1] = 0;
   Xf[3][1] = 0;
   Xf[4][1] = 0;
   Xf[5][1] = 0;
   Xf[6][1] = 0;
   Xf[7][1] = 0;
   
    fs = 1000;

    printcomplex(Xf,N,1,2);
     fv=Getfassoc(fs, N);
     tv=Gettassoc(fs, N);
  System.out.print("\n f= ");
   printarray(fv, N,  0);
  System.out.print("\n t= ");
   printarray(tv, N,  0);

  
    Yn=CDFT(N,Xf);
  System.out.print("\n CDFT(XF)= ");
    printcomplex(Yn, N, 5, 6);
    Idft=CIDFT(N, Yn );
  System.out.print("\n CIDFT(Yn)= ");
    printcomplex(Idft, N, 7, 8);
      for(i=0;i<N;i++)
    {
    Xn1[i]=Idft[i][0];
    }

     Hilb=Hilbert(Xn1, N );
  System.out.print("\n Hilbert(Xn1)= ");
     printcomplex(Hilb, N, 9, 10);
     SpX=GetSpectrum(Xn1, N , 1);
  System.out.print("\n Get spectrum= ");
     printcomplex(SpX, N, 11, 12);
    
   System.out.print("\n Get envelope (Xn1) ");
     env=getenvelope(Xn1, N);
     printarray(env, N, 13);
    
     System.out.print("\nPoke new  aray ");

    //  cos (90° – x) = sin x
    // -sin ( x-90) = cos x 
    // Hilbert(sin(x)) = -cos(t)
    // Hilbert(Cos(x)) = sin(x)
   
   // For i = 0 To n - 1
   // Xi1(i + 1) = 1 * Sin(2 * M_PI * 15.625 * i / fs + 0) + 0.5 * Cos(2 * M_PI * 93.75 * i / fs + 0)
   // Next i
    
     N = 16;

    

  
    double [][] Xn2=new double[N][2] ;
    double [][] Xf2=new double[N][2] ;
     for (i= 0;i<N; i++)
      {
      Xf2[i][0] = 0;
      Xf2[i][1] = 0;
      }
     
     Xf2[1][0] = 1;
     Xf2[1][1] = -45;
     Xf2[5][0] = 0;
     Xf2[5][1] = 0;
     Xn2 =SynthSignalfromRFi(Xf2, N, 1);
    System.out.print("\n Get Xn from Xf");
    printcomplex(Xn2, N, 14, 15);
    fv=Getfassoc(fs, N);
   System.out.print("\n f[]=");
     printarray(fv, N, 16);
      System.out.print("\n SpX[]=");  
       SpX=GetSpectrumfromIQ(Xn2, N, 1);
       printcomplex(SpX, N, 17, 18);

   
   double[] Xn4=new double[N]; 
        for (i = 0 ; i<N; i++)
       {
        Xn4[i] = 0.5*Math.cos(2*Math.PI*250*i/fs+0);
        Xn4[i] +=1.5*Math.cos(2*Math.PI*125*i/fs+Math.PI/2);
        Xn4[i] +=2*Math.cos(2*Math.PI*375 *i/fs-Math.PI/2);
       }

    double [][] Xn3=new double[N][2] ;
    double [][] Xf3=new double[N][2] ;
    double [][] Yn3=new double[N][2] ;
   System.out.print("\n Spectrum (Xn4)=");    
       Xf3=GetSpectrum(Xn4, N , 1);
       printcomplex(Xf3, N, 19, 20);
    System.out.print("\n Spectrum IQ(Xn4)=");
       Yn3=GetSpectrumIQ(Xn4, N );
       printcomplex(Yn3, N, 21, 22);
   return;

  }

  
  public int  GetNumF1(double f1 ,int Numb  ,double fs )
  {
   return (int) (Math.floor(f1*Numb/fs));
  }

  public static double GetTHD_F(double[] xn , int Numb, double fs, double f1)
   {

  int numbf1;
  double[][] Xf =new double[Numb][2];
  double Urms1h;
  double Uhharm;
  int i ;

  //
  //f 0  5  10  15  20   25   30   35 Hz  N=8
  //i 0  1   2   3   4   5    6     7
  //         f1      2f1    3f1
  //      numb1=2
  //
  //f 0  5  10  15  20   25   30   35    40   45  50  55   60 Hz   N=13 imax=12
  //i 0  1   2   3   4   5    6     7    8    9   10  11   12
  //            f1           2f1             3f1           4f1
  //      numb1=3
  //

  numbf1=(int)Math.floor(f1*Numb/fs);//   2  , get associated vector of frequency
  Xf=GetSpectrum(xn, Numb,  0);
  Urms1h= Xf[numbf1][0]/Math.sqrt(2);
  i = numbf1*2;
   
  Uhharm=0;
  while (i<Numb)
  {
   if (i<Numb) { Uhharm+=(Xf[i][0]*Xf[i][0])/2 ; } 
   i=i+numbf1;
   
  }

  Uhharm =Math.sqrt(Uhharm) ;
   
  return Uhharm*100/(Urms1h);
  }

  
  public static double GetTHD_R(double[] xn , int Numb, double fs, double f1)
  {
  int numbf1;
  double[][]  Xf =new double[Numb][2];
  double Urmssigah;
  double Uhharm;
  int i ;

   
  numbf1=(int)Math.floor(f1*Numb/fs);  //   get associated vector of frequency, ORIGIN=0
  Xf=GetSpectrum(xn, Numb, 0);
  /* parser of Urms of all harmonics */

  // THD-R=Urmshharm*100/urmsig =sqrt( sum(Xrms[i] ^2))/Vrmssig;
  Urmssigah=0;
  i=numbf1;

   while (i<Numb)
   {
   if (i<Numb) { Urmssigah+=(Xf[i][0]*Xf[i][0])/2;  }
   i=i+1;
   }

  Urmssigah=Math.sqrt (Urmssigah );
  Uhharm=0;
   /* parser of Urms of higher harmonics */

   i=numbf1*2; // ' number of  bin of 2-nd harmonic 
    while (i<Numb)
   {
   if (i<Numb) { Uhharm+=(Xf[i][0]*Xf[i][0])/2;  }
   i=i+numbf1;         
   }

   Uhharm=Math.sqrt (Uhharm) ;
   
   return (Uhharm*100/Urmssigah);
  }

  public static double GetRMS(double[] xn , int Numb )
  {
  int i ;
  double Urms;
  Urms = 0;
  for (i=0;i<Numb;i++)
  {
   Urms+=xn[i]*xn[i];
  }
  Urms=Math.sqrt(Urms/Numb);

  return Urms;
  }

  
  public static double GetTHD_N(double[] xn , int Numb, double fs  , double f1 )
  {
  int numbf1;
  double[][] Xf=new double [Numb][2];
   
  double Urmssig;
  double Uhharmn;
  int i ;
  /*
  Urmssig = 0;
  for(i=0 ;i<Numb;i++)
  {
   Urmssig+=xn[i]*xn[i];
  }
  Urmssig=Math.sqrt(Urmssig/Numb);

  */
   numbf1=(int)Math.floor (f1*Numb/fs);
   Xf=GetSpectrum(xn, Numb, 0);

  i=0;
  Urmssig=0;
   while (i<Numb)
  {
  { Urmssig+=Xf[i][0]*Xf[i][0]/2;   }
  i=i+1;
  }

  Urmssig= Math.sqrt(Urmssig);

  Uhharmn=0;
  i=0;

   while (i<Numb)
  {
  if (i!= numbf1) { Uhharmn+=Xf[i][0]*Xf[i][0]/2;   }
  i=i+1;
  }

  Uhharmn= Math.sqrt(Uhharmn);
  return (Uhharmn*100/Urmssig);
  }

  public static void PrintValueDouble(String str, double s,String str1)
  {
   
  System.out.print(str);
  System.out.print(String.format("%.08f",s));
   System.out.print(str1);
  return;
  }

  public static void PrintValueInt(String str, int s,String str1)
  {
   
  System.out.print(str);
  System.out.print(String.format("%d",s));
   System.out.print(str1);
  return;
  }

  public static  double GetValueDouble(String s,Scanner scan1)
  {
   
  System.out.print(s);

  double result;
  result=(double)  scan1.nextDouble();
   
  return result;
  }

  
  public static int GetValueInt(String s,Scanner scan1)
  {
  System.out.print(s);
   
  int result=(int) scan1.nextInt();
   
  return result;
  }

  public static void TestTHDN( )
   {

  
   int N;
   Scanner sc2 = new Scanner (System.in);
  PrintValueInt("\n N=",16," ");
   N =16;// GetValueInt("\n Input N ",sc2);

   double [][] Xf =new double[N][2];
   double [][] X_iq=new  double[N ][2];
   double [][]  Xn =new  double[N ][2] ;
   double []  Xn1=new  double[N ] ;
   double [] fvect  =new double[N ] ;
   double fs;
   double f1;
   int i;

   fs =44100;// GetValueDouble("\n Input fs,Hz ",sc2);
   f1 =11025;//GetValueDouble("\n Input f [1-st harm] ,Hz ",sc2);
   fvect= Getfassoc(fs,  N);

   for(i = 0;i<N;i++)
   {
   
   //PrintValueDouble("\n f=",fvect[i]);
   Xf[i][0] =0;//GetValueDouble("\n input Um[f] " ,sc2);
   Xf[i][1] =0;//GetValueDouble("\n input fi[f] ",sc2);

   }
  PrintValueDouble("\n fs=",fs," ");
  PrintValueDouble("\n f1=",f1," ");

  PrintValueDouble("\n f=",fvect[4]," Hz");
  Xf[4][0]=1;
   PrintValueDouble("  |X[f]|=",Xf[4][0]," ");

  PrintValueDouble("\n f=",fvect[8]," Hz");
  Xf[8][0]=0.3;
   PrintValueDouble("  |X[f]|=",Xf[8][0]," ");

  PrintValueDouble("\n f=",fvect[9]," Hz");
  Xf[9][0]=0.1;
   PrintValueDouble("  |X[f]|=",Xf[9][0]," ");

   for(i = 0;i<N;i++)
   {
   
   PrintValueDouble("\n f=",fvect[i]," ");
   PrintValueDouble("  |X[f]|=",Xf[i][0], " ");
   PrintValueDouble("  fi=",Xf[i][1]," "); 

   }

  
    Xn=SynthSignalfromRFi(Xf,  N, 1);
    printcomplex(Xn, N, 5, 6);
   for (i=0;i<N;i++)
  {
  Xn1[i]=Xn[i][0];
  }
   
   double THD_N= GetTHD_N(Xn1, N, fs, f1);
   PrintValueDouble("\n THD+N= ", THD_N," %");
   double THD_R= GetTHD_R(Xn1, N, fs, f1);
    PrintValueDouble("\n THD-R= ",THD_R," %");
  double THD_F= GetTHD_F(Xn1, N, fs, f1 );
   PrintValueDouble("\n THD-F= ",THD_F," %");
  double RMS=GetRMS(Xn1, N);
   PrintValueDouble("\n U RMS=",RMS," units of ADC");
    sc2.close();
   return ;
   }

  }

  public   class  THDObtainer
   {
  public static void main(String[] args) throws IOException
   {
          THDParser.TestHilbert();
          THDParser.TestTHDN();
          }

  }

  Friday, January 26, 2018 2:08 PM
• 

  

  0

  Sign in to vote

  spectrumstat - utility for obtaining max. levels of signals on frequencies of the spectrum using short-time Fourier    transform and founder of the max. vaues

  xinputns.txt

  4,410000e+04   16  8
       
  0   1,000000e+01  
   
  1   -7,071068e+00  
   
  2   3,061617e-15  
   
  3   7,071068e+00  
   
  4   -1,000000e+01  
   
  5   7,071068e+00  
   
  6   -2,694842e-14  
   
  7   -7,071068e+00
  8   1,000000e+01  
   
  9   -7,071068e+00  
   
  10   3,061617e-15  
   
  11   7,071068e+00  
   
  12   -1,000000e+01  
   
  13   7,071068e+00  
   
  14   -2,694842e-14  
   
  15   -7,071068e+00   

  4,410000e+04   16

  maxspectr.txt
   
  0   3,020133e-15  0,000000e+00 
   
  1   1,330312e-07  5,512500e+03 
   
  2   6,447902e-15  1,102500e+04 
   
  3   1,000000e+01  1,653750e+04 
   
  4   4,642237e-15  2,205000e+04 
   
  5   5,151483e-15  2,756250e+04 
   
  6   2,334110e-15  3,307500e+04 
   
  7   7,772354e-15  3,858750e+04 

  make.bat

  set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"

   
   javac Spectrumstat.java
   jar cvfe Spectrumstat.jar Spectrumstat   *.class
    jar uf  Spectrumstat.jar Spectrumstat.java xinputns.txt
  pause 0

  run.bat

  set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"
   
   java -jar  Spectrumstat.jar  
  pause 0

  //javac Spectrum.java
  //java Spectrum

  import java.util.Scanner;
  import java.io.*;

  
   public class Spectrumstat
  {

        

    public static double[][]   cdft(double[][] xn  ,int N )
      {
    
     int j;
     double[][] Xk=new  double[N][2];
     int k,n  ;
     double angle;

  
   for (k=0;k<N;k++)
    {
    Xk[k][0]=0;
    Xk[k][1]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*Math.PI*k*n/N;  
      Xk[k][0]+= xn[n][0]*Math.cos(angle)-xn[n][1]*Math.sin(angle);
      Xk[k][1]+= xn[n][0]*Math.sin(angle)+xn[n][1]*Math.cos(angle);
      }
   }

   return Xk;
   }

  
   public static double[][] ScaleDFT(double[][] xk,int  num)
    {
   
    for (int i=0; i<num ; i++)
     {
      xk[i][0] =xk[i][0] /num;
      xk[i][1] =xk[i][1] /num;
     }
   return xk;
   }

   public static double[][] hilbert(double[] xr  )
        {
        
        double angle;
        double[]  h=new double[xr.length+1]  ;
        double[]  dftxn_re=new double[xr.length]  ;
        double[]  dftxn_im=new double[xr.length]  ;
        double[]  dftXk_re=new double[xr.length] ;
        double[]  dftXk_im=new double[xr.length] ;
        double[] idftXK_re=new double[xr.length] ;
        double[] idftXK_im=new double[xr.length] ;
        double[][] idftXN=new double[xr.length][2] ;

   
    if (xr.length==0)
    {
    System.out.println("\n Error");
   
    idftXN[0][0] =0;
    idftXN[0][1]=0;
    return idftXN;
    }
    //xnreal=real(xr);
   
   
    for (int i=0;i<xr.length;i++)
    {
     dftxn_re[i]=xr[i];
     dftxn_im[i]=0;
    }

   //fft
   // dftXk=dft(dftxn,n);

   for (int k=0;k<xr.length;k++)
    {
    dftXk_re[k]=0;
    dftXk_im[k]=0;
    for(int n=0;n<xr.length;n++)
      {
      angle=-2*Math.PI*k*n/xr.length;
      dftXk_re[k]+= dftxn_re[n]*Math.cos(angle)- dftxn_im[n]*Math.sin(angle);
      dftXk_im[k]+= dftxn_re[n]*Math.sin(angle)+ dftxn_im[n]*Math.cos(angle);
      }

   }

   //X={Xk.re,Xk.im};
   
   for (int i=0;i<=xr.length;i++)
    {
    h[i]=0;
    }
  //int no2=(int) Math.floor (xr.length/2);

    if((xr.length>0)&&((2* Math.floor(xr.length/2))==(xr.length))) // n is even
    {

    h[0]=1;
     h[xr.length/2]=1;
     for(int k=1;k<xr.length/2;k++)
     {
     h[k]=2;
     }
    }
    else //n is odd
    {
      if(xr.length>0)
      {
       h[0]=1;
       for(int k=1;k<(xr.length+1)/2;k++)
       {
       h[k]=2;
       }
      }
    }

   for(int i=0;i<xr.length;i++)
    {
    idftXK_re[i]= dftXk_re[i]*h[i];
    idftXK_im[i]= dftXk_im[i]*h[i];
    }

  
   //idftXN= idft(idftXK, n);
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (int n=0;n<xr.length;n++)
    {
    idftXN[n][0] =0;
    idftXN[n][1] =0;
    for(int k=0;k<xr.length;k++)
     {
    
     angle=2*Math.PI*k*n/xr.length;
     idftXN[n][0] +=idftXK_re[k]*Math.cos(angle)-idftXK_im[k]*Math.sin(angle);
     idftXN[n][1] +=idftXK_re[k]*Math.sin(angle)+idftXK_im[k]*Math.cos(angle);
     }
   idftXN[n][0] =idftXN[n][0] /xr.length;
   idftXN[n][1] =idftXN[n][1] /xr.length;

   }

  
    return idftXN;
   }

     public static  double[][] getanalyticsignal(double [] xn )
       {
      
    int i;
      
          double[][] Hilb =new double[xn.length][2];
          double[][] x_n  =new double[xn.length][2];

   Hilb=hilbert(xn );
   
    for (i=0; i<xn.length; i++)
    {
    x_n[i][0] =xn[i];
     x_n[i][1] =Hilb[i][1] ;
    }
   
   return x_n;
   }

   
    public static  double Arg(double re, double im)
         {
   return Math.atan2(im,re)  ; 
   }

   public static double[][]  get_mod_fi(double[][] x, int n  )
     {
     
     double[][] res=new double[ n][2] ;
     double fi;
     for (int i=0;i<n;i++)
     {
     res[i][0] =Math.sqrt(x[i][0] *x[i][0] +x[i][1] *x[i][1] );
     res[i][1] =Arg(x[i][0] ,x[i][1] );
       { res[i][1] =res[i][1] *180/Math.PI;  }        
     }
     return  res;
   }

  
   public static double[][]  GetSpectrum(   double [] xn  )
   {

    
    
           double[][] xn_iq=new double[xn.length][2];
    double[][] Xf_IQ =new double[xn.length][2] ;
    double[][] Xf_modfi =new double[xn.length][2]; 
   
    xn_iq =getanalyticsignal( xn );
    Xf_IQ =cdft(xn_iq, xn.length);
    Xf_IQ =ScaleDFT(Xf_IQ,xn.length);
    Xf_modfi =get_mod_fi(Xf_IQ, xn.length );
   
   
          
            
   return   Xf_modfi ;
   }

  
   public static double[][][]  GetSpectrumarray(   double [] xn , int num,int nseq)
  {
    //int np=(int) num/nseq;
    int npages=(int) Math.floor (num/nseq);

    double [] xnseq =new double [nseq];
    
    double [][]  Xf = new double  [nseq][2];
    double [][][] Xfarray= new double [npages][nseq][2];

  
   for (int i=0;i<npages; i++) //0,1,32
   {
    for (int j=0;j<nseq;j++)
     {
     xnseq[j]=xn[i+j];     //if  nseq=16, xnseq[j]=xn[k],k= 0...15,16...31
     // if ((i+j)<num) &&( (nseq*npages)>(num/nseq) ) xnseq[j]=0;
     }

    Xf=GetSpectrum(xnseq );
    for (int j=0;j<nseq; j++) //0,1,32
    {
    Xfarray[i][j][0]= Xf[j][0];
    Xfarray[i][j][1]=Xf[j][1];
    }

   }

  
  return Xfarray;
  }

  public static double get_max(double a,double b )
  {
  if (a>b) {return a;}

  return b;
  }

  
   public static double[] GetPeakSpectrumComponents(double [] xn , int num, int nseq)
  {
   //int np=(int)num/nseq;
   int npages=(int) Math.floor (num/nseq);
   double [][][] Xfarray =new double [npages][nseq][2];
   double [] Xfpeak=new double[nseq];
   Xfarray=GetSpectrumarray( xn , num, nseq);

   for (int j=0;j<nseq;j++) // for each Xn [frequency](page)
   {
   Xfpeak[j]=0;
   for (int i=0;i<npages;i++) //for pages of spectrum
     {
     Xfpeak[j]=get_max(Xfpeak[j], Xfarray[i][j][0] );
     }

   }

  return Xfpeak;

  }

    public static void main(String[] args) throws IOException
   {

                      Scanner sc = new Scanner(System.in);

     
                   int NMAX=32767;
                   int N=0,n,i, ns=0;
                  double[] xn =new double[NMAX] ;
                  
                    double fs=0,T;
                           
           System.out.print ("\n For input  x[n] from file input '1' ");
    System.out.print ("\n For manual  input  x[n] input   '0' ");
    int key=sc.nextInt();
           
    switch( key)
          {

   

    case 1:
   
             System.out.print ("\n  Open xinput.txt \n");
             String fileNameDefined = "xinputns.txt";
             File file1 = new File(fileNameDefined);
             
          try{
              // -read from filePooped with Scanner class
              Scanner inputStream  = new Scanner(file1);
              fs = (double) inputStream.nextDouble();
              N =inputStream.nextInt();
              ns =inputStream.nextInt(); 
               System.out.print ("\n ns=");
               System.out.print (ns);
             for (n=0 ; n<N ;n++)
      {
               inputStream.nextInt();
            
      xn[n]= inputStream.nextDouble();

      }

             
              
              inputStream.close();

  
          } catch (FileNotFoundException e) {

              e.printStackTrace();
          }

  
    break;

    case 0:

     System.out.print ("\n Input fs input,Hz  ");
     fs = (double) sc.nextDouble();
     System.out.print ("\n Input number of samples of xin[n]  ");
        N =sc.nextInt();
              System.out.print ("\n Input number of samples in CDFT sequence  ");
        ns=sc.nextInt();
               xn =new double[N] ;
     for (n=0 ; n<N ;n++)
     {
     T=n*1/(N *fs );
     System.out.print ("\n n=");
     System.out.print (n);
     System.out.print (" t=");
     System.out.print (T);
     System.out.print (" s; xn[n]=  " ); 
     xn[n]= sc.nextDouble();
     }
      break;
    default: 
        break;
   }
    System.out.print ("\n  fs =");
     System.out.print (fs);
    System.out.print (" Hz");
     System.out.print ("\n Nsamples1=");
     System.out.print (N);

   for (n=0 ; n<N ;n++)
   {
    T=n*1/(N *fs );
     System.out.print ("\n n=");
     System.out.print (n);
     System.out.print (" x[n]=");
     System.out.print (xn[n]);
     System.out.print ("  t= ");
     System.out.print (T);
     System.out.print (" s");
   }
   
     System.out.println ("\n Input 0 for continue ");
    sc.nextLine();
     System.out.print ("\n OK");

   

                    System.out.println("Array of Xn is: ");
      for (  i=0;i<N;i++)
                    {
                     System.out.print("\n");
                     System.out.print("N= ");
                     System.out.print(i);                
                     System.out.print (" ; xn[n]= ");
                     System.out.print(xn[i]);
                    }

                   double Xfmax[]  =new  double[ns]  ;
                  System.out.println("\n Parsing spectrum ");
                  Xfmax=GetPeakSpectrumComponents(xn, N, ns );
                   //Xf=GetSpectrum(xn, N);

  
                
              System.out.println("\n Array of Xfmax is: ");
      for (i=0;i<ns ;i++)
                    {
                     System.out.print("\n");
                     System.out.print("ns= ");
                     System.out.print(i);                

                      System.out.print (" ; max |Xf[k]|= ");
                      System.out.print(Xfmax[i]   );
                      
                      System.out.print (" ; f=");
                      System.out.print ( fs*i/ns );
                      System.out.print(" Hz");

                    }

      File outFile = new File ("maxspectr.txt");
      FileWriter fWriter = new FileWriter (outFile);
      PrintWriter pWriter = new PrintWriter (fWriter);

  
      pWriter.println("\n");
      pWriter.print(String.format("%e",fs ));
      pWriter.print("   ");
      pWriter.print(String.format("%d",N ) );
                     
      pWriter.println("\n");

  
     for (i=0; i<ns ;i++)
     { 

   
           pWriter.println("\n ");
    pWriter.print(String.format("%d",i) );
           pWriter.print("   "); 
    pWriter.print(String.format("%e",Xfmax[i] ));
    pWriter.print("  ");
           pWriter.print (String.format("%e", fs*i/ns) );
           pWriter.print("  ");    
     }
   
   
      pWriter.close();

                  System.out.println("\n End of sequence ");

                     sc.close();
   }

  }

   

  Friday, January 26, 2018 3:01 PM
• 

  

  0

  Sign in to vote

  You may convert Java programs into VC++ console applications (if the problem of work with big 2-dimension arrays is solved ).

  Friday, January 26, 2018 3:04 PM
• 

  

  0

  Sign in to vote

  Program for obtaining envelope of sinal sequence using Hilbert transform .

  // getenv.cpp : Defines the entry point for the console application.
  //

  #include "stdafx.h"

  #include <stdio.h>
  #include <stdlib.h>
  #include <math.h>
  #include <iostream>

  
  #define NMAX 1000
  #define M_PI 3.1415926535897932384626433832795

  typedef struct {
  double re[NMAX];
  double im[NMAX];
  } CData;

   typedef struct  {
       double data[NMAX];
   } Array;

  
  CData  cdft(CData xn  ,int N )
   {
     
    
   CData Xk ;
   int k,n  ;
   double angle;

   for (k=0;k<N;k++)
    {
    Xk.re[k]=0;
    Xk.im[k]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*M_PI*k*n/N;  
      Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
      Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
      }

   }

   return Xk;
   }

   CData  cidft(CData idft_Xk, int N  )
   {

   int k,n ;
   double angle;
   CData  idft_xn;
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    idft_xn.re[n]=0;
    idft_xn.im[n]=0;
    for(k=0;k<N;k++)
     {
   
     angle=2*M_PI*k*n/N;
     idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
     idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
     }

   idft_xn.im[n]=idft_xn.im[n]/N;
   idft_xn.re[n]=idft_xn.re[n]/N;

   }

  return idft_xn ;
   }

  CData ScaleDFT(CData xk,int  num)
   {

   int i;
   for (i=0; i<num ; i++)
     {
      xk.re[i]=xk.re[i]/num;
      xk.im[i]=xk.im[i]/num;
     }
   return xk;
   }

  CData ScaleIQidft(CData Xk, int num)
   {
   int i;
     for (i=0;  i<num ; i++)
    {
     Xk.re[i]=Xk.re[i]*num;
     Xk.im[i]=Xk.im[i]*num;
    }

  return Xk;

  }

   CData hilbert(double *xr , int N)
   {
    int  k,n,i, no2 ;
    double angle;
    double h[NMAX];
    double dft__xn_re[NMAX];
    double dft__xn_im[NMAX];
     double dft__Xk_re[NMAX];
     double dft__Xk_im[NMAX];
       double idft__XK_re[NMAX];
       double idft__XK_im[NMAX];
     CData idft__XN;

   
    if (N==0)
    {
    printf("\n Error");
    idft__XN.re[0]=0;
    idft__XN.im[0]=0;
    return idft__XN;
    }
    //xnreal=real(xr);
    no2=(int)floor ((double)N/2);
   
    for (i=0;i<N;i++)
    {
     dft__xn_re[i]=xr[i];
     dft__xn_im[i]=0;
    }

   //fft
   // dft__Xk=dft(dft__xn,n);

  for (k=0;k<N;k++)
    {
    dft__Xk_re[k]=0;
    dft__Xk_im[k]=0;
    for(n=0;n<N;n++)
      {
      angle=-2*M_PI*k*n/N;
      dft__Xk_re[k]+= dft__xn_re[n]*cos(angle)- dft__xn_im[n]*sin(angle);
      dft__Xk_im[k]+= dft__xn_re[n]*sin(angle)+ dft__xn_im[n]*cos(angle);
      }

   }

   //X={Xk.re,Xk.im};

   for (i=0;i<=N;i++)
    {
    h[i]=0;
    }

  
    if((N>0)&&((2*floor((double)N/2))==(N))) // n is even
    {

    h[0]=1;
     h[N/2]=1;
     for(k=1;k<N/2;k++)
     {
     h[k]=2;
     }
    }
    else //n is odd
    {
      if(N>0)
      {
       h[0]=1;
       for(k=1;k<(N+1)/2;k++)
       {
       h[k]=2;
       }
      }
    }

   for(i=0;i<N;i++)
    {
    idft__XK_re[i]= dft__Xk_re[i]*h[i];
    idft__XK_im[i]= dft__Xk_im[i]*h[i];
    }

  
   //idft__XN= idft(idft__XK, n);
   
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    idft__XN.re[n]=0;
    idft__XN.im[n]=0;
    for(k=0;k<N;k++)
     {
     //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
     angle=2*M_PI*k*n/N;
     idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
     idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
     }
   idft__XN.im[n]=idft__XN.im[n]/N;
   idft__XN.re[n]=idft__XN.re[n]/N;

   }

  
  return idft__XN;
   }

  CData getanalyticsignal(double *xn, int num)
   {
   CData x_n,Hilb;
   int i;
     //x=x+jy,y=0
     //hilb(x)=a+jb
     //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
     //I=x(t), Q=j*Hilb(xre(t))

  Hilb=hilbert(xn,num);
   
   for (i=0; i<num; i++)
   {
    x_n.re[i]=xn[i];
    x_n.im[i]=Hilb.im[i];
   }
   
   return x_n;
   }

  
   
  double Arg(double re, double im)
  {
   
  return atan2 (im,re)  ; 
  }

  CData  get_mod_fi(CData x, int n, int mode   )
   {
    int i;
    CData res;
    double fi;
    for (i=0;i<n;i++)
    {
    res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
    //double treshh=1e9;
    // x.im[i]=x.im[i]*treshh;
    // x.im[i]=floor(x.im[i])/treshh;
    // x.re[i]=x.re[i]*treshh;
    // x.re[i]=floor(x.re[i])/treshh;
    res.im[i]=Arg(x.re[i],x.im[i]);
    if (mode==1)  { res.im[i]=res.im[i]*180/M_PI;  }
            
    }
    return  res;
   }

  CData GetSpectrum(double *xn, int num)
  {
   
  CData xn_iq, Xf_IQ, Xf_modfi;
  xn_iq=getanalyticsignal( xn,  num);
   
  Xf_IQ=cdft(xn_iq  ,num );
  Xf_IQ=ScaleDFT(Xf_IQ,num);
  Xf_modfi=get_mod_fi(Xf_IQ, num,1   );
   
  return Xf_modfi;
  }

  Array getenvelope(double *x, int num)
   {
   int i;
   CData  y;
   Array env;
   y=hilbert(x,num);
   for (i=0;i<num;i++)
    {
    env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
    }
    return env;
   }
   
   
   
  double getmax( Array   k , int N)
  {
  double kmax=k.data[0];
  double k1=0;
  int i;
  for (i=1;i<N;i++)
  {
   k1=k.data[i];
    if(fabs( k1 )>fabs(kmax)) { kmax= k.data[i]; } 
    
  }
   
   return kmax;
  }
   
  double getmin( Array  k, int N)
  {
  double kmin=k.data[0];
  int i;
  for (i=0;i<N;i++)
  {
    if(fabs(k.data[i])<fabs(kmin)) { kmin= k.data[i]; } 
  }
   return kmin;
  }

  double getmid( Array  k,int N)
  {
  double kmid=0;
  int i;
  for (i=0;i<N;i++)
  {
   kmid+=   k.data[i];
  }
  kmid=kmid/N;
   return kmid;
  }

  
  FILE   *fp1,*fp2 ;

  int _tmain(int argc, _TCHAR* argv[])
  //int main ()
  {
  double delta =1e-15;
  double xin[NMAX];
  double yout[NMAX];
  Array  Xenv1  ;
  double Kmax=0, Kmin=0, Kmid=0;
  int key;
  int k,n,i,N1,N2;
  double fs1,fs2,T,x,f;

  
  /***********input x[n] *************/
  label_00:
  printf ("\n For input  x[n] from file input '1' ");
  printf ("\n For manual  input  x[n] input   '0' ");
  std::cin>>key;
  if (key== 1  ) { goto label_01; }
  if (key== 0 ) { goto label_02; }

  label_01:

   printf("\n Open xinput.txt for reading x[n]   ");
  if ((fp1=fopen("xinput.txt","r"))==NULL)  printf ("\n File opening error ");
  fscanf(fp1,"%le",&fs1);
  fscanf(fp1,"%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   fscanf (fp1,"%d",&k);
   fscanf (fp1,"%le",&xin[n]);
  }

  
  fclose(fp1);

  goto label_03;

  label_02:

  printf("\n Input fs input,Hz  " ) ;
  scanf( "%le",&fs1);
  printf("\n Input number of samples of xin[n]  ") ;
  scanf( "%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   T=n*1/(N1*fs1); 
   printf ("\n n=%d ; t=%le s; xin[n]= ",n,T);
   scanf ( "%le",&xin[n]);
  }

   label_03:
  printf( "\n fs1=%le Hz ",fs1);
  printf( "\n Nsamples1= %d ",N1);

  for (n=0 ; n<N1;n++)
  {
   T=n*1/(N1*fs1); 
   printf("\n n=%d x[n]=%le,  t=%le s",n,xin[n],T);
  }

  printf("\n Input 0 for continue ");
  std::cin>>key;
  printf("\n OK.   ");

   

  
  /*****************************/

  label1:
   
  printf("\n Obtaining envelope of x[n]");
  Xenv1=getenvelope(xin, N1);

  printf("\n |Hilbert(x[n])| : ");
   
  if ((fp2=fopen("envelope.txt","a"))==NULL)  printf ("\n File opening error");
  fprintf(fp2,"\n\n   %le   %d \n",fs1, N1);

  
  for (i=0; i<N1;i++)
  {
   
   
  printf ("\n %d   |Hilbert(x[n])|= %le  ",  i, Xenv1.data[i] );
  fprintf (fp2,"\n %d   %le  ",  i, Xenv1.data[i]);

  }
  fclose(fp2);

  printf ("\n  input   '0' for continue ");
  std::cin>>key;

  Kmax=getmax(Xenv1,N1);
  Kmin=getmin(Xenv1,N1);
  Kmid=getmid(Xenv1,N1);
  printf ("\n  Kmax= %le ",   Kmax);
  printf ("\n  Kmin= %le ",   Kmin);
  printf ("\n  Kmid= %le ",   Kmid);

  if ((fp2=fopen("coeffs.txt","a"))==NULL)  printf ("\n File opening error");
  fprintf(fp2,"\n\n   %le   %d \n",fs1,N1);
  fprintf (fp2,"\n  Kmax= %le ",   Kmax);
  fprintf (fp2,"\n  Kmin= %le ",   Kmin);
  fprintf (fp2,"\n  Kmid= %le ",   Kmid);
  fclose(fp2);

  printf ("\n  input   '0' for exit ");
  std::cin>>key;

  return 0;

  }

  // stdafx.cpp : source file that includes just the standard includes
  // getenv.pch will be the pre-compiled header
  // stdafx.obj will contain the pre-compiled type information

  #include "stdafx.h"

  // TODO: reference any additional headers you need in STDAFX.H
  // and not in this file

  // stdafx.h : include file for standard system include files,
  // or project specific include files that are used frequently, but
  // are changed infrequently
  //

  #pragma once

  
  #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
  #include <stdio.h>
  #include <tchar.h>

  // TODO: reference additional headers your program requires here

  Friday, January 26, 2018 4:19 PM
• 

  

  0

  Sign in to vote

  Program for obtaining of the Bessel's  functions for obtainer of the PWM and FM spectrums of  signals.

  // bessel.cpp : Defines the entry point for the console application.
  //

  #include "stdafx.h"

  
  #include <stdio.h>
  #include <stdlib.h>
  #include <iostream>
  #include <math.h>
  #include "bessel.h"

  
  int main()

  {
  int n;
  double x;
  double result1, result2,result3,result4;

  std::cout<<"Input order n : "<<"\n";
  std::cin>>n;
  std::cout<<"Input x : "<<"\n";
  std::cin>>x;
  result1=jn(n,x);
  std::cout<<"jn(n,x)= "<<result1<<"\n";
  result2=bessj( n, x );
  std::cout<<"jn(n,x)= "<<result2<<"\n";

  
  result3=yn(n,x);
  std::cout<<"yn(n,x)= "<<result3<<"\n";
  result4=bessy(  n,  x );
  std::cout<<"yn(n,x)= "<<result4<<"\n";

  std::cout<<"Input '0' for exit";
  char key;
  std::cin>>key;
  return 0;
  }

  //bessel.h (add into project using Add Item-> .h file,then  put the code into file)

  #include <stdio.h>
  #include <stdlib.h>
  #include <math.h>

  
  #define ACC 40.0
  #define BIGNO 1.0e10
  #define BIGNI 1.0e-10

  /*
  http://www.astro.rug.nl/~gipsy/sub/bessel.c
  */

  static double bessj0( double x )
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate Bessel function of first kind and order  */
  /*          0 at input x                                      */
  /*------------------------------------------------------------*/
  {
     double ax,z;
     double xx,y,ans,ans1,ans2;

     if ((ax=fabs(x)) < 8.0) {
        y=x*x;
        ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
           +y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
        ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
           +y*(59272.64853+y*(267.8532712+y*1.0))));
        ans=ans1/ans2;
     } else {
        z=8.0/ax;
        y=z*z;
        xx=ax-0.785398164;
        ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
           +y*(-0.2073370639e-5+y*0.2093887211e-6)));
        ans2 = -0.1562499995e-1+y*(0.1430488765e-3
           +y*(-0.6911147651e-5+y*(0.7621095161e-6
           -y*0.934935152e-7)));
        ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
     }
     return ans;
  }

  static double bessj1( double x )
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate Bessel function of first kind and order  */
  /*          1 at input x                                      */
  /*------------------------------------------------------------*/
  {
     double ax,z;
     double xx,y,ans,ans1,ans2;

     if ((ax=fabs(x)) < 8.0) {
        y=x*x;
        ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
           +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
        ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
           +y*(99447.43394+y*(376.9991397+y*1.0))));
        ans=ans1/ans2;
     } else {
        z=8.0/ax;
        y=z*z;
        xx=ax-2.356194491;
        ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
           +y*(0.2457520174e-5+y*(-0.240337019e-6))));
        ans2=0.04687499995+y*(-0.2002690873e-3
           +y*(0.8449199096e-5+y*(-0.88228987e-6
           +y*0.105787412e-6)));
        ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
        if (x < 0.0) ans = -ans;
     }
     return ans;
  }

  double bessj( int n, double x )
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate Bessel function of first kind and order  */
  /*          n at input x                                      */
  /* The function can also be called for n = 0 and n = 1.       */
  /*------------------------------------------------------------*/
  {
     int    j, jsum, m;
     double ax, bj, bjm, bjp, sum, tox, ans;

  
     if (n < 0)
     {
        double   dblank;
        printf("Error");
      //  setdblank_c( &dblank );
        return( 0 );
     }
     ax=fabs(x);
     if (n == 0)
        return( bessj0(ax) );
     if (n == 1)
        return( bessj1(ax) );

     if (ax == 0.0)
        return 0.0;
     else if (ax > (double) n) {
        tox=2.0/ax;
        bjm=bessj0(ax);
        bj=bessj1(ax);
        for (j=1;j<n;j++) {
           bjp=j*tox*bj-bjm;
           bjm=bj;
           bj=bjp;
        }
        ans=bj;
     } else {
        tox=2.0/ax;
        m=2*((n+(int) sqrt(ACC*n))/2);
        jsum=0;
        bjp=ans=sum=0.0;
        bj=1.0;
        for (j=m;j>0;j--) {
           bjm=j*tox*bj-bjp;
           bjp=bj;
           bj=bjm;
           if (fabs(bj) > BIGNO) {
              bj *= BIGNI;
              bjp *= BIGNI;
              ans *= BIGNI;
              sum *= BIGNI;
           }
           if (jsum) sum += bj;
           jsum=!jsum;
           if (j == n) ans=bjp;
        }
        sum=2.0*sum-bj;
        ans /= sum;
     }
     return  x < 0.0 && n%2 == 1 ? -ans : ans;
  }

  
  static double bessy0( double x )
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate Bessel function of second kind and order */
  /*          0 at input x.                                     */
  /*------------------------------------------------------------*/
  {
     double z;
     double xx,y,ans,ans1,ans2;

     if (x < 8.0) {
        y=x*x;
        ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
           +y*(10879881.29+y*(-86327.92757+y*228.4622733))));
        ans2=40076544269.0+y*(745249964.8+y*(7189466.438
           +y*(47447.26470+y*(226.1030244+y*1.0))));
        ans=(ans1/ans2)+0.636619772*bessj0(x)*log(x);
     } else {
        z=8.0/x;
        y=z*z;
        xx=x-0.785398164;
        ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
           +y*(-0.2073370639e-5+y*0.2093887211e-6)));
        ans2 = -0.1562499995e-1+y*(0.1430488765e-3
           +y*(-0.6911147651e-5+y*(0.7621095161e-6
           +y*(-0.934945152e-7))));
        ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
     }
     return ans;
  }

  static double bessy1( double x )
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate Bessel function of second kind and order */
  /*          1 at input x.                                     */
  /*------------------------------------------------------------*/
  {
     double z;
     double xx,y,ans,ans1,ans2;

     if (x < 8.0) {
        y=x*x;
        ans1=x*(-0.4900604943e13+y*(0.1275274390e13
           +y*(-0.5153438139e11+y*(0.7349264551e9
           +y*(-0.4237922726e7+y*0.8511937935e4)))));
        ans2=0.2499580570e14+y*(0.4244419664e12
           +y*(0.3733650367e10+y*(0.2245904002e8
           +y*(0.1020426050e6+y*(0.3549632885e3+y)))));
        ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x);
     } else {
        z=8.0/x;
        y=z*z;
        xx=x-2.356194491;
        ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
           +y*(0.2457520174e-5+y*(-0.240337019e-6))));
        ans2=0.04687499995+y*(-0.2002690873e-3
           +y*(0.8449199096e-5+y*(-0.88228987e-6
           +y*0.105787412e-6)));
        ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
     }
     return ans;
  }

  
  double bessy( int n, double x )
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate Bessel function of second kind and order */
  /*          n for input x. (n >= 0)                           */
  /* Note that for x == 0 the functions bessy and bessk are not */
  /* defined and a blank is returned.                           */
  /*------------------------------------------------------------*/
  {
     int j;
     double by,bym,byp,tox;

  
     if (n < 0 || x == 0.0)
     {
        double   dblank;
      printf("Error"); // setdblank_c( &dblank );
        return(0 );
     }
     if (n == 0)
        return( bessy0(x) );
     if (n == 1)
        return( bessy1(x) );

     tox=2.0/x;
     by=bessy1(x);
     bym=bessy0(x);
     for (j=1;j<n;j++) {
        byp=j*tox*by-bym;
        bym=by;
        by=byp;
     }
     return by;
  }

  static double bessi0( double x )
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate modified Bessel function In(x) and n=0.  */
  /*------------------------------------------------------------*/
  {
     double ax,ans;
     double y;

  
     if ((ax=fabs(x)) < 3.75) {
        y=x/3.75,y=y*y;
        ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
           +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
     } else {
        y=3.75/ax;
        ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
           +y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
           +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
           +y*0.392377e-2))))))));
     }
     return ans;
  }

  
  static double bessi1( double x)
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate modified Bessel function In(x) and n=1.  */
  /*------------------------------------------------------------*/
  {
     double ax,ans;
     double y;

  
     if ((ax=fabs(x)) < 3.75) {
        y=x/3.75,y=y*y;
        ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
           +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
     } else {
        y=3.75/ax;
        ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
           -y*0.420059e-2));
        ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
           +y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
        ans *= (exp(ax)/sqrt(ax));
     }
     return x < 0.0 ? -ans : ans;
  }

  
  double bessi( int n, double x)
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate modified Bessel function In(x) for n >= 0*/
  /*------------------------------------------------------------*/
  {
     int j;
     double bi,bim,bip,tox,ans;

  
     if (n < 0)
     {
        double   dblank;
       printf("Error");// setdblank_c( &dblank );
        return( 0);
     }
     if (n == 0)
        return( bessi0(x) );
     if (n == 1)
        return( bessi1(x) );

  
     if (x == 0.0)
        return 0.0;
     else {
        tox=2.0/fabs(x);
        bip=ans=0.0;
        bi=1.0;
        for (j=2*(n+(int) sqrt(ACC*n));j>0;j--) {
           bim=bip+j*tox*bi;
           bip=bi;
           bi=bim;
           if (fabs(bi) > BIGNO) {
              ans *= BIGNI;
              bi *= BIGNI;
              bip *= BIGNI;
           }
           if (j == n) ans=bip;
        }
        ans *= bessi0(x)/bi;
        return  x < 0.0 && n%2 == 1 ? -ans : ans;
     }
  }

  
  static double bessk0( double x )
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate modified Bessel function Kn(x) and n=0.  */
  /*------------------------------------------------------------*/
  {
     double y,ans;

     if (x <= 2.0) {
        y=x*x/4.0;
        ans=(-log(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420
           +y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
           +y*(0.10750e-3+y*0.74e-5))))));
     } else {
        y=2.0/x;
        ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
           +y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
           +y*(-0.251540e-2+y*0.53208e-3))))));
     }
     return ans;
  }

  
  static double bessk1( double x )
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate modified Bessel function Kn(x) and n=1.  */
  /*------------------------------------------------------------*/
  {
     double y,ans;

     if (x <= 2.0) {
        y=x*x/4.0;
        ans=(log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144
           +y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1
           +y*(-0.110404e-2+y*(-0.4686e-4)))))));
     } else {
        y=2.0/x;
        ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619
           +y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2
           +y*(0.325614e-2+y*(-0.68245e-3)))))));
     }
     return ans;
  }

  double bessk( int n, double x )
  /*------------------------------------------------------------*/
  /* PURPOSE: Evaluate modified Bessel function Kn(x) and n >= 0*/
  /* Note that for x == 0 the functions bessy and bessk are not */
  /* defined and a blank is returned.                           */
  /*------------------------------------------------------------*/
  {
     int j;
     double bk,bkm,bkp,tox;

  
     if (n < 0 || x == 0.0)
     {
        double   dblank;
       printf("Error");// setdblank_c( &dblank );
        return( 0 );
     }
     if (n == 0)
        return( bessk0(x) );
     if (n == 1)
        return( bessk1(x) );

     tox=2.0/x;
     bkm=bessk0(x);
     bk=bessk1(x);
     for (j=1;j<n;j++) {
        bkp=bkm+j*tox*bk;
        bkm=bk;
        bk=bkp;
     }
     return bk;
  }

  
  #undef ACC
  #undef BIGNO
  #undef BIGNI

  // stdafx.cpp : source file that includes just the standard includes
  // bessel.pch will be the pre-compiled header
  // stdafx.obj will contain the pre-compiled type information

  #include "stdafx.h"

  // TODO: reference any additional headers you need in STDAFX.H
  // and not in this file

  // stdafx.h : include file for standard system include files,
  // or project specific include files that are used frequently, but
  // are changed infrequently
  //

  #pragma once

  
  #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers
  #include <stdio.h>
  #include <tchar.h>
  //#include "bessel.h"

  
  // TODO: reference additional headers your program requires here

  
  

  • Edited by USERPC01 Friday, January 26, 2018 6:57 PM bug fix

  Friday, January 26, 2018 6:54 PM
• 

  

  0

  Sign in to vote

  project for obtaining IR coefficients of FIR from AFR data   with #include <vector>,  #include <complex>, #include <cmath>

  //afrtoht.cpp

  // afrtoht.cpp : main project file.

  #include "stdafx.h"
  #include <stdio.h>
  /*
  #include <tchar.h>
   #include <iostream>
  #include <stdlib.h>
  #include <cmath>
  */
  #include <iostream>
  #include <fstream>
  #include <complex>
  #include <vector>
  #include <set>
  #include <cmath>
  #include <limits>

  using namespace std;

  #define M_PI 3.1415926535897932384626433832795
  //int _tmain(int argc, _TCHAR* argv[])
  //{
  // return 0;
  //}

   //afrtohtcoefs.cpp
  /* Define complex double data type */
   typedef complex<double> dcomp;

   //std::vector<std::vector<dcomp> >
   typedef std::vector<dcomp> cplxvector ;

  #define  NMAX 512
  //FILE *fp, *fp1;
   
  int k,n;
  int w,t,N;
  double h[NMAX];
  double h_im[NMAX];
  double Re_H[NMAX];
  double Im_H[NMAX];
  double Re_x[NMAX];
  double Im_x[NMAX];
  double abs_H[NMAX];
  double arg_H[NMAX];
  double f[NMAX];
  float angle;
  double   K, fs,T,finp,Tmax;
  char str;
  char key;

  /*
  float heavi( double t)
  {
  if (t<0 ) return 0;
  if (t>=0) return 1;
  }

  float heavi1( int n)
  {
  if (n<0 ) return 0;
  if (n>=0) return 1;
  }*/

  double getarg(float Re_H,float Im_H)
  {

  return atan2(Im_H,Re_H) * 180 / M_PI;
  }

  cplxvector CIDFT(cplxvector idftXk, int N)
  {

   
   double angle;
   cplxvector  idftxn;
   dcomp data;

  //N=sizeof(xk)/sizof(double);
   //ifft
    for (int n=0;n<N;n++)
    {
    data=complex<double>(0.0,0.0);
    
   
    for(int k=0;k<N;k++)
     {
      angle=2*M_PI*k*n/N;
      data+=idftXk.at(k) *  exp(complex<double> (0.0,angle));
     // data+=idft_Xk.at(k)* cmath.exp( (0.0,angle));
     //data+=idft_Xk.at(k)*  exp( 1i*angle );
     //idft_xn.re[n]+= idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
     // idft_xn.im[n]+= idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
      }

  data/=N;
  idftxn.push_back( data);

   }
  return idftxn;
  }

  
   int   main( )
  {

   fstream finp;
   std::cout <<  "\n For input  |H(jw)|, fi(w) from file input '1' " ;
   std::cout <<  "\n For manual  input |H(jw)|, fi(w) input   '0' " ;

      std::cin>>key;
    ofstream fout ("hcoefs1.txt" , ios::ate | ios::app);
       fout.close();
     if (key=='0' ) { goto label_02;}
     if (key=='1' ) { goto label_01;}

  label_01:;
   std::cout << "\n Open afrcoefs.txt for reading AFR coefficients ";
   
  // ifstream finp("afrcoefs.txt", ios::in );
   finp.open("afrcoefs.txt", ios::in );
   
    finp>>fs ;
    finp>>N ;
  Tmax=N*1/fs;
  for (k=0 ; k<N;k++)
  {
  // f[k]=k*fs/N;
   
    //fout<<"\n n="<<k<<" |H(jw)|="<<abs_H[k]<< "; arg(H(jw))= "<< arg_h[k]<<"; f="<<f[k]<<" Hz"  ;
    finp>> n ;
    finp>>abs_H[k] ;
    finp>>arg_H[k] ;
   
  //  finp>>f[k] ;
   
   arg_H[k]*= M_PI/180;
  }
  finp.close();
  fout.close();
   
  goto label_03;

  
  label_02:;

  std::cout<<"\n Input fs ,Hz " ;
  std::cin>>fs;
  std::cout<<" Input  number of points " ;
  std::cin>>N ;
  Tmax=N*1/fs;
  std::cout<<"\n Tmax ="<<Tmax<<" s" ;
  for (k=0 ; k<N;k++)
  {
   f[k]=k*fs/N;

  //fout<<"\n n="<<k<<"  |H(jw)|=" <<  abs_H[k]<<" ; arg(H(jw))="<< arg_h[k]<<"; f="<<f[k]<<" Hz" ;
  //fout<<"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
     std::cout<<"\n F="<<f[k]<<" Hz" ;
     std::cout<<"\n k="<<k<<";    |H(jw)|[k]=" ;
     std::cin>> abs_H[k] ;
     std::cout<<" k="<<k<<" ;   arg(H(jw))[k]=" ;
      std::cin>>arg_H[k] ;
    //abs_H[k]=1;
    //arg_H[k]=0;
  arg_H[k]*= M_PI/180;

  }

  
  label_03:
   
  for (k=0 ; k<N;k++)
  {
  f[k]=k*fs/N;
   std::cout << "\n n="<<k<<";|H(jw)|= "<<abs_H[k]<< " ;arg(H(jw))="<<arg_H[k]*180/M_PI<<"  deg;f="<< f[k]<<" Hz";
   
  }

  
  std::cout<< "\n Input 0 for continue " ;
  std::cin.get();
  std::cout<<"\n OK. Parsing ..." ;
   
  // H(jw)=integr(0;inf;h(tau)*exp(-j*w*tau) ; d tau)
  //K=2*3.1415926/fs;
   /*
  DFT
  X[k]=sum(n=0;N-1; x[n]*exp(-2*pi*j*k*n/N) ), k=0,...,N-1
   IDFT
  x[n]=(1/N)*sum(k=0;N-1;X[k]*exp(2*pi*k*n/N) ), n=0,...,N-1
   100*10^-6*500
  */

  //label1:;

  
  //cidft
  //printf("\n Open hcoefs1.txt for writing coeffs h[n] ");
  

  
   
   
   cplxvector  idft_xn, idft_Xk;
   dcomp data;

   for(int k=0;k<N;k++)
    {
     Re_H[k]=abs_H[k]*cos(arg_H[k]);
     Im_H[k]=abs_H[k]*sin(arg_H[k]);
     data =complex<double> (Re_H[k] ,Im_H[k]) ; 
    // std::cout<<"\n k="<<k<< "  Xk[k]=  " <<data  ; /*real(idft_xn.at(n) ) */   ;
     idft_Xk.push_back(data); 
   }

  
   std::cout <<"\n IDFT\n";
   idft_xn= CIDFT(idft_Xk,N);

  fstream fout1 ("hcoefs1.txt" , ios::ate | ios::app);
    for (int n=0;n<N;n++)
    {
    std::cout<<"\n n="<<n<< "  h[n]=  " <<real(idft_xn.at(n) )  ;
    fout1 << "\n"<<n <<"  "<< real(idft_xn.at(n) )   ; 
    }
     fout1.close();

   //fclose (fp1);
   std::cout<< "\n  Re(h(n)=h[n]), use it to load into FIR  "   ;
  //label3:;
  std::cout<<"\n Input 0 for exit " ;
  std::cin.get();
  return 0;
  }

  //stdafx.cpp

  // stdafx.cpp : source file that includes just the standard includes
  // afrtoht.pch will be the pre-compiled header
  // stdafx.obj will contain the pre-compiled type information

  #include "stdafx.h"

  //stdafx.h

  // stdafx.h : include file for standard system include files,
  // or project specific include files that are used frequently, but
  // are changed infrequently
  //

  #pragma once

  
   #define WIN32_LEAN_AND_MEAN  // Exclude rarely-used stuff from Windows headers

  
  // TODO: reference additional headers your program requires here

  //afrcoefs.txt

  44100
  4
  0  1  0
  1  1  0
  2  0  0
  3  0  0

  //hcoefs1.txt

  
  0  0.5
  1  0.25
  2  0
  3  0.25

  
  

  • Edited by USERPC01 Tuesday, February 6, 2018 6:22 PM fix bugs

  Tuesday, February 6, 2018 6:20 PM
• 

  

  0

  Sign in to vote

  Example of using <vector>, <complex> for  Hilbert transform

  (fix bugs)

  // cplxarray.cpp : main project file.

  #include "stdafx.h"
   
  #include <iostream>
  //#include <tgmath>//for the type generate marcros.

  #include <complex>
  #include <vector>
  #include <set>
  #include <cmath>
  #include <limits>
  /*
  #include <cstdio>
  #include <cstdlib>
  #include <cmath>
  #include <limits>
  #include <iostream>
   */
  using namespace std;

  /* Define complex double data type */
   typedef complex<double> dcomp;

   //std::vector<std::vector<dcomp> >
   typedef std::vector<dcomp> cplxvector ;

  #define NMAX 2048

  #define M_PI 3.1415926535897932384626433832795
  /*
  typedef struct {

  double re[NMAX];
  double im[NMAX];
  } CData;

   typedef struct  {
       double data[NMAX];
   } Array;

  */

  cplxvector cdft( cplxvector xn  ,int N )
   {
      
   
   cplxvector Xk ;
   dcomp data;
   int k,n  ;
   double angle;
   
   for (k=0;k<N;k++)
    {
    //Xk.re[k]=0;
    //Xk.im[k]=0;
   data=(0.0,0.0);
    for(n=0;n<N;n++)
      {
       angle=-2*M_PI*k*n/N;  
       //Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
       //Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
    //data+=xn.at(n)*cmath.exp( (0.0,angle));

    data+=xn.at(n)* exp( (0.0,angle));
      }

   }

   return Xk;
   }

   cplxvector  cidft(cplxvector idft_Xk, int N  )
   {

   int k,n ;
   double angle;
   cplxvector  idft_xn;
   dcomp data; 
   //N=sizeof(xk)/sizof(double);
   //ifft
    for (n=0;n<N;n++)
    {
    data=(0.0,0.0);
    
   
    for(k=0;k<N;k++)
     {
   
      angle=2*M_PI*k*n/N;
   //creal(z1),
   //cimag(z1),
  // data+=idft_Xk.at(k)* cmath.exp( (0.0,angle));
   data+=idft_Xk.at(k)*  exp( (0.0,angle));
     //data+=idft_Xk.at(k)*  exp( 1i*angle );
     //idft_xn.re[n]+= idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
     // idft_xn.im[n]+= idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
      }

  data/=N;
  idft_xn.push_back( data);
  // idft_xn.im[n]=idft_xn.im[n]/N;
  //  idft_xn.re[n]=idft_xn.re[n]/N;

   }

  return idft_xn ;
   }

  
   cplxvector  ScaleDFT( cplxvector  xk, int  num)
   {

   int i;
   for (i=0; i<num ; i++)
     {
       xk.at(i)/= num;
    //  xk.im[i]=xk.im[i]/num;
     }
   return xk;
   }

   cplxvector  ScaleIQidft( cplxvector  Xk, int num)
   {
   int i;
      for (i=0;  i<num ; i++)
    {
    Xk.at(i)*=num;
     //Xk.re[i]=Xk.re[i]*num;
     //Xk.im[i]=Xk.im[i]*num;
    }

  return Xk;

  }
   
    cplxvector  hilbert(double *xr , int N)
   {
    int  k,  i, no2 ;
    //double angle;
    double h[NMAX];
    cplxvector dft__xn;
    //double dft__xn_re[NMAX];
    //double dft__xn_im[NMAX];
    cplxvector dft__Xk;
     //double dft__Xk_re[NMAX];
     //double dft__Xk_im[NMAX];
    cplxvector idft__XK;
      // double idft__XK_re[NMAX];
      // double idft__XK_im[NMAX];
      cplxvector  idft__XN;

   
    if (N==0)
    {
    printf("\n Error");
     idft__XN.push_back((0,0)) ;
   
    return idft__XN;
    }
    //xnreal=real(xr);
    no2=(int) floor ((double)N/2);
   // dcomp  tmpdata;

    for (i=0;i<N;i++)
    {
     // tmpdata=(xr[i],0) ;
       dft__xn.push_back((xr[i],0))  ;
     //dft__xn_re[i]=xr[i];
        //dft__xn_im[i]=0;
    }

   //fft
   // dft__Xk=dft(dft__xn,n);
  dft__Xk=cdft( dft__xn  , N );

   //X={Xk.re,Xk.im};

   for (i=0;i<=N;i++)
    {
    h[i]=0;
    }

  
    if((N>0)&&(( 2*floor((double)N/2) )==(N))) // n is even
    {

    h[0]=1;
     h[N/2]=1;
     for(k=1;k<N/2;k++)
     {
     h[k]=2;
      }
    }
    else //n is odd
    {
       if(N>0)
      {
        h[0]=1;
       for(k=1;k<(N+1)/2;k++)
       {
        h[k]=2;
       }
       }
    }

   for(i=0;i<N;i++)
    {
    idft__XK.at(i)*=h[i];
    //idft__XK_re[i]= dft__Xk_re[i]*h[i];
    //idft__XK_im[i]= dft__Xk_im[i]*h[i];
    }

  
   //idft__XN= idft(idft__XK, n);
   
    //N=sizeof(xk)/sizof(double);
   //ifft
  idft__XN=cidft(idft__XK,  N  );

  
   /*
    for (n=0;n<N;n++)
    {
    idft__XN.re[n]=0;
    idft__XN.im[n]=0;
    for(k=0;k<N;k++)
     {
      //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
     angle=2*M_PI*k*n/N;
     idft__XN.re[n]+=idft__XK_re[k]*cos(angle)-idft__XK_im[k]*sin(angle);
      idft__XN.im[n]+=idft__XK_re[k]*sin(angle)+idft__XK_im[k]*cos(angle);
      }
   idft__XN.im[n]=idft__XN.im[n]/N;
    idft__XN.re[n]=idft__XN.re[n]/N;

   }
  */
   
  return idft__XN;
   }
   
   
   cplxvector  getanalyticsignal(double *xn, int num)
   {
    cplxvector  x_n,Hilb;
   int i;
     //x=x+jy,y=0
     //hilb(x)=a+jb
     //z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
      //I=x(t), Q=j*Hilb(xre(t))

  Hilb=hilbert(xn,num);
   dcomp data1;
   double Ximg;
   for (i=0; i<num; i++)
    {
     data1=Hilb.at(i) ;
     Ximg=data1.imag();
     data1=(xn[i],Ximg );
     x_n.push_back(data1);
    //x_n.re[i]=xn[i];
    //x_n.im[i]=Hilb.im[i];
   }
   
   return x_n;
   }

   /*
   
  double Arg(double re, double im)
  {
   
  return   atan2 (im,re)  ;
  }
  */

  cplxvector  get_mod_fi(cplxvector x, int n, int mode   )
    {
    int i;
    cplxvector res;
    double fi;
    double module;
    for (i=0;i<n;i++)
    {
     module=abs(x.at(i));
     fi=arg(x.at(i));
    //res.re[i]=//pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
    //double treshh=1e9;
    // x.im[i]=x.im[i]*treshh;
    // x.im[i]=floor(x.im[i])/treshh;
    // x.re[i]=x.re[i]*treshh;
    // x.re[i]=floor(x.re[i])/treshh;
    //res.im[i]=Arg(x.re[i],x.im[i]);
     if (mode==1)  { fi=fi*180/M_PI;  }
      res.push_back((module,fi));
  std::cout<< "\n f="<<  " Hz;|Xf(jw)|=" << module << ";arg(Xf)="<< fi<<"  deg " ;
    }
     return  res;
   }

   
   
  cplxvector  GetSpectrum(double *xn, int num)
  {
   
  cplxvector  xn_iq, Xf_IQ, Xf;
  xn_iq=getanalyticsignal( xn,  num);
   
  // Xf_IQ=cdft(xn_iq  ,num );
   //Xf_IQ=ScaleDFT(Xf_IQ,num);
   // Xf =get_mod_fi(Xf_IQ, num,1   );
   
  return Xf;
  }

  
  /*

  
  Array getenvelope(double *x, int num)
   {
   int i;
    CData  y;
   Array env;
   y=hilbert(x,num);
    for (i=0;i<num;i++)
    {
    env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);
    }
     return env;
   }
  */

  
  //FILE   *fp1,*fp2 ;

  
  int main()
  {

  double xin[NMAX];
  cplxvector  Xf;

  char key;
  int   n,i,N1;
  //int k;
  double fs,T, f;
  //double x;

  
  //printf ("\n For input  x(t) from file input '1' ");
  printf ("\n For manual  input  x(t) input   '0' ");
  std::cin>>key;
  //if (key=='1' ) { goto label_05;}
  if (key=='0' ) { goto label_06;}
  /*

  //label_05:

   printf("\n Open xinput.txt for reading x[n]   ");
  if ((fp1=fopen("xinput.txt","r"))==NULL)  printf ("\n File opening error ");
  fscanf(fp1,"%le",&fs);
  fscanf(fp1,"%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   fscanf (fp1,"%d",&k);
    fscanf (fp1,"%le",&xin[n]);
  }

  
  fclose(fp1);

  goto label_07;

  */

  label_06:

  printf("\n Input fs,Hz , ( %lf ) ",fs) ;
  std::cin>>fs; //scanf( "%le",&fs);
   printf("\n Input number of samples of xin[n]  ") ;
   std::cin>>N1; //scanf( "%d",&N1);

  for (n=0 ; n<N1;n++)
  {
   T=n*1/(N1*fs);
   std::cout<<"\n n="<<n<<"; t="<<T<<" s;   xin[n]= " ;
   std::cin>>xin[n];
   }

  // label_07:;

  std::cout<< "\n fs=" <<fs << " Hz " ;
  std::cout<< "\n Nsamples= "<<N1 ;

  for (n=0 ; n<N1;n++)
  {
   T=n*1/(N1*fs);
    std::cout<< "\n n="<< n << "; x[n]=" <<xin[n] << ",  t="<<T<<" s" ;
  }

  std::cout<<"\n Input 0 for continue " ;
  std::cin.get();
  std::cout<<"\n OK.   " ;

  
  /*****************************/

  //label1:;
  //Xf.re= |Xf|,  Xf.im -arg(Xf) (only in this program)
   std::cout<<"\n Obtaining spectrum of x[n]" ;
  Xf=GetSpectrum(xin, N1);

   std::cout<< "\n Spectrum of x[n] : " ;

   //if ((fp2=fopen("spectrum.txt","a"))==NULL)  printf ("\n File opening error");
   //fprintf(fp2,"\n   %le   %d \n",fs,N1);
   //fprintf(fp2, "\n n;\t f, Hz;\t |Xf(jw)| ;\t  arg(Xf), deg ");
  for (i=0; i<N1;i++)
  {
  f=fs*i/N1;
  //cout <<"\n f="<< f<<" Hz;|Xf(jw)|="<<Xf.re[i]<<";arg(Xf)="<<Xf.im[i]<<" deg ";
  // std::cout<< "\n f="<<f<<" Hz;|Xf(jw)|=" <<   real(Xf.at(i))<< ";arg(Xf)="<< imag(Xf.at(i))<<"  deg " ;
   //fprintf(fp2,"\n %d   %le  %le",i ,Xf.re[i],Xf.im[i]);

  }

  std::cout<< "\n  input   '0' for exit " ;
   
  std::cin.get();

  return 0;

  }

  Tuesday, February 6, 2018 6:25 PM
• 

  

  0

  Sign in to vote

  // RectAIMSpectr.cpp : 

  #include "stdafx.h"
  #include <stdio.h>
  #include <stdlib.h>
  #include <iostream>
  #include <math.h>

  using namespace std;
  #define M_PI 3.1415926535897932384626433832795

  
  double tau0, deltatau, A0;
  double f0, Fm,f1;
  double w0,T0,m,m1,Wm;

   

  
  int k,n;

  #define kmax 60 //change
  #define nmax 60 //change 

  double U0;
  double UFm[nmax], Ffm[nmax];
  double fi[kmax], Ufi[kmax];
  double fk[2*nmax+1][kmax+1], Ufk[2*nmax+1][kmax+1];
  double Q0;

  
  double Q;
  int kmax1;

  void InputData3()
  {

  
   //Q0=T0/tau0;
   //T0=1/f0
   //tau0=T0/Q0;

  std::cout<<"\nInput A0,V : (1) ";
  std::cin>>A0;
  std::cout<<"\nInput f1,Hz : (16384) ";
  std::cin>>f1;

  std::cout<<"\nInput Q  : (10) ";
  std::cin>>Q;
   std::cout<<"\nInput kmax  : (60) ";
  std::cin>>kmax1;
   
  return;
  }
  /*
   ________                    _______
  |             |                    |
  |             |___________|

  */
  void getFreqRectQ(void)
  {

  double U0=A0/Q;
  double timp=1/(f1*Q);
  double w1=2*M_PI*f1;
   
  double Ti=1/f1;
  std::cout <<"\n            f0="<<0<<"Hz ; U0=" <<U0<<" V  ";

  for(k=1;k<=kmax1; k++)
  {
  //Ufi~ A0*|sinc(k*pi/q)|  =A0*sin (k*pi/q) /(k*pi/q)
  //Ufi[k]=2*A0*sin(k*M_PI/Q)/(k*M_PI);
  //fi[k]=k*f1;

  //std::cout <<"\n k="<<k<<";  f[k]="<<fi[k]<<"Hz ; Um[k]=" <<Ufi[k]<<" V  ";

  double Um;
  //Um=fabs( ( 2*A0 /(k*M_PI))*sin(k*M_PI/Q));
   
  //std::cout <<"\n k="<<k<<";  f[k]="<<k*f1<<"Hz ; Um[k]=" <<Um<<" V  ";
   Um=(2*timp/Ti)*A0*(sin(k*2*M_PI*f1*timp/2)/(k*2*M_PI*f1*timp/2));
   Um=fabs(Um);
   std::cout <<"\n k="<<k<<";  f[k]="<<k*f1<<"Hz ; Um[k]=" <<Um<<" V  ";

  }

   std::cin.get() ;

  
  return;
  }

  /********************************************/
  /*

  S(t)=A0+sum(from k=1;to k=m;{ Ak*cos(k*wi*t)})=
  =(taui*Ui/Ti)+(2*taui/Ti*Ui)*sum(from k=1;to k=m;{
  (sin(k*wi*taui/2)/(k*wi*taui/2))*cos(kwit))

  uaim(t)=(taui*Ui/Ti)+(taui/Ti)*Ui*maim*cos(W*t)+
  +(2*taui/Ti)*Ui*sum(from k=1;to k=m;{
  (sin(k*wi*taui/2)/(k*wi*taui/2))*(1+maim*cos(W*t))*cos(k*wi*t) })

  Uk*cos(k*wi*t)+0.5*maim*Uk*cos(k*wi-W)*+0.5*maim*Uk*cos(k*wi+W)*t
  Uk=(2*taui/Ti)*Ui*(sin(k*wi*taui/2)/(k*wi*taui/2))

  f=0 Um[0]= taui*Ui/Ti
  f=W/2pi Um[W]=taui*miam*Ui/Ti

  f=k*wi-W Um=0.5*maim*Uk
  f=k*wi    Uk=Uk
  f=k*wi+W Um=0.5*maim*Uk
  ...
  Uenv_level[2*pi/taui]=0

  fsr=Whi

  */

  double maim;

  void InputData4()
  {

  
   //Q0=T0/tau0;
   //T0=1/f0
   //tau0=T0/Q0;

  std::cout<<"\nInput A0,V : (1) ";
  std::cin>>A0;
  std::cout<<"\nInput f1,Hz : (16384) ";
  std::cin>>f1;

  std::cout<<"\nInput Q  : (10) ";
  std::cin>>Q;
  std::cout<<"\nInput Fmod ,Hz : (1000) ";
  std::cin>>Fm;
  std::cout<<"\nInput maim,%  : (10) ";
  std::cin>>maim;
  maim/=100;
   std::cout<<"\nInput kmax  : (60) ";
  std::cin>>kmax1;
   
  return;
  }

  
  void getFreqRectAIM(void)
  {

   
  double timp=1/(f1*Q);
  double w1=2*M_PI*f1;
  double Ti=1/f1;

   
   
   
   double U0=timp*A0/Ti;
   std::cout <<"\n  ;  f[0]="<<0<<"Hz ; Um[k]=" << U0<<" V  ";
   double Fmod=Fm;
   double Ufm=timp*maim*A0/Ti;
   std::cout <<"\n  ;  f[Fm]="<<Fmod<<"Hz ; Um[k]=" << Ufm<<" V  ";
   double Uk,Ulsbk,Uusbk;
   double fk,flsbk,fusbk;
    std::cout <<"\n  ;  f[2*pi/taui  /( 2*pi)]="<< 1/timp<<"Hz ; Um[1/tau]=" << 0<<" V  ";
   for (int k=1;k<kmax1;k++){
   
   
   Uk=(2*timp/Ti)*A0*(sin(k*2*M_PI*f1*timp/2)/(k*2*M_PI*f1*timp/2));
   Uk=fabs(Uk);
   Ulsbk=0.5*maim*Uk;
   Uusbk=0.5*maim*Uk;
   fk=k*f1;
   flsbk=fk-Fm;
   fusbk=fk+Fm;
   
   std::cout <<"\n  ;  f[lsbk]="<<flsbk<<"Hz ; Um[flsbk]=" <<Ulsbk<<" V  ";
   std::cout <<"\n  ;  f[k]=   "<<fk<<" Hz   ; Um[fk]=" <<Uk<<" V  ";
   std::cout <<"\n  ;  f[usbk]="<<fusbk<<"Hz ; Um[fusbk]=" <<Uusbk<<" V  ";
   
   
   }

   std::cin.get() ;

  return;
  }

  int main( )
  {
   
   int key;
   std::cout<<"\n Input 1 for rect. signal\n Input 2 for AIM ";
   std::cin>>key;
   
   switch(key)
   {

  
   case 1:
           std::cout<<"\n  Rect. signal with Q   ";
          InputData3();
          getFreqRectQ();
         break;
   case 2:
       InputData4();
      getFreqRectAIM( );
     break;
   }

    std::cout<<"\n Input 0 for exit ";
    std::cin.get() ;
  return 0;
  }

  Tuesday, February 6, 2018 9:12 PM