none
Obtainer of Hilbert transform (C++) RRS feed

  • Question

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


    //obtainer of Hilbert transform . Using scilab algorithm,
    // translated into C++  language
    //Zaporozhye, 23.10.2015


    double xnre[2048],xnim[2048];
    double Xkre[2048], Xkim[2048];
    double XNRE[2048],XNIM[2048];
    double XKRE[2048], XKIM[2048];

    void getfft(  int N1)
    {

    int k,n,N;
    double angle;

    N=N1;
    //N=sizeof(xr)/sizof(double);
    //fft
     for (k=0;k<N;k++)
     {
     Xkre[k]=0;
     Xkim[k]=0;
     for(n=0;n<N;n++)
      {
      angle=2*M_PI*k*n/N;
      Xkre[k]=Xkre[k]+xnre[n]*cos(angle);
      Xkim[k]=Xkim[k]-xnre[n]*sin(angle);

      //XkRE[n]=XkRE[n]+XnRE[n]*cos(angle)+XnIM[n]*sin(angle);
      //XkIM[n]=XkIM[n]-XnRE[n]*sin(angle)-XKIM[n]*cos(angle);


      }

     }

    return  ;
    }

    void doifft(   int N1)
    {

    int k,n,N;
    double angle;

    N=N1;
    //N=sizeof(xk)/sizof(double);
    //ifft
     for (n=0;n<N;n++)
     {
     XNRE[n]=0;
     XNIM[n]=0;
     for(k=0;k<N;k++)
      {
    //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
      angle=2*M_PI*k*n/N;
    //a=xkre
    //b=xkim
    //c=cos
    //d=sin


      XNRE[n]=XNRE[n]+XKRE[k]*cos(angle)-XKIM[k]*sin(angle);
      XNIM[n]=XNIM[n]+XKRE[k]*sin(angle)+XKIM[k]*cos(angle);
      }
    XNIM[n]=XNIM[n]/N;
    XNRE[n]=XNRE[n]/N;
     }

    return  ;
    }


    void hilbert(double xr[2048], int N)
    {
    int n,k,i,no2;
    double h[2048];

    //n=sizeof(xr);
    n=N;
    if (n==0)
    {
    //free(xr);
    printf("\n Error");
    return;

    }

    //xnreal=real(xr);
    no2=int (n/2);

    for (i=0;i<n;i++)
    {
    xnre[i]=xr[i];

    }


    //fft

    getfft(n);
    //X={Xkre,Xkim};


    for(i=0;i<n;i++)
    {
    printf("\n fft(m)[ %d ]= %le +j* %le ",i, Xkre[i],Xkim[i]);
    }


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

    h[i]=0;
    }


    if((n>0)&&((2*floor(n/2))==(n))) // n is even
    {
    printf("\n n is even ");
     h[0]=1;
    // h[n/2]=1;
     h[n/2]=1;
     
     for(k=1;k<n/2;k++)
     {
     h[k]=2;
     }
    }
     
    else //n is odd
    {
    printf("\n n is odd ");
     if(n>0)
     {
     h[0]=1;
     
     for(k=1;k<(n+1)/2;k++)
     {
     h[k]=2;
     }
    }

    }

     //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
    // (a+bi)*c=ac+bc*i
    //Xktw[i]=x[i]*h[i];
    for(i=0;i<n;i++)
    {
    XKRE[i]=Xkre[i]*h[i];
    XKIM[i]=Xkim[i]*h[i];
    }

    for (i=0;i<n;i++)
    {
    printf("\n twiddle [ %d ]=%le +j %le ",i,XKRE[i],XKIM[i] );
    }


    doifft( n);
    printf("\n y=hilbert(xr) ");
    for (i=0;i<n;i++)
    {
    printf("\n y[ %d ]=%le +j %le  ",i,XNRE[i],XNIM[i] );
    }

    return ;
    }


    int main(void)
    {
    double m[8]={1,2,3,4,5,6,7,8};
    int num=8;
    int j;
    printf("\n*************\n n=%d  \n m=\n",num);
     
    for (j=0;j<num;j++)
    {

    printf(" \n m[%d]= %le ",j,m[j]);
    }

    hilbert(m, num);

    getch();
     // m[]={1,2,3,4,5,6,7,8};
      num=7;
    printf("\n*************\n n=%d  \n m=\n",num);
    for (j=0;j<num;j++)
    {

    printf("\n m[%d]= %le ",j,m[j]);
    }

    hilbert(m, num);

    return 0;
    }

    Friday, October 23, 2015 3:58 AM

All replies

  • *************
     n=8
     m=

     m[0]= 1.000000e+000
     m[1]= 2.000000e+000
     m[2]= 3.000000e+000
     m[3]= 4.000000e+000
     m[4]= 5.000000e+000
     m[5]= 6.000000e+000
     m[6]= 7.000000e+000
     m[7]= 8.000000e+000
     fft(m)[ 0 ]= 3.600000e+001 +j* 0.000000e+000
     fft(m)[ 1 ]= -4.000000e+000 +j* 9.656854e+000
     fft(m)[ 2 ]= -4.000000e+000 +j* 4.000000e+000
     fft(m)[ 3 ]= -4.000000e+000 +j* 1.656854e+000
     fft(m)[ 4 ]= -4.000000e+000 +j* -3.918740e-015
     fft(m)[ 5 ]= -4.000000e+000 +j* -1.656854e+000
     fft(m)[ 6 ]= -4.000000e+000 +j* -4.000000e+000
     fft(m)[ 7 ]= -4.000000e+000 +j* -9.656854e+000
     n is even
     twiddle [ 0 ]=3.600000e+001 +j 0.000000e+000
     twiddle [ 1 ]=-8.000000e+000 +j 1.931371e+001
     twiddle [ 2 ]=-8.000000e+000 +j 8.000000e+000
     twiddle [ 3 ]=-8.000000e+000 +j 3.313708e+000
     twiddle [ 4 ]=-4.000000e+000 +j -3.918740e-015
     twiddle [ 5 ]=-0.000000e+000 +j -0.000000e+000
     twiddle [ 6 ]=-0.000000e+000 +j -0.000000e+000
     twiddle [ 7 ]=-0.000000e+000 +j -0.000000e+000
     y=hilbert(xr)
     y[ 0 ]=1.000000e+000 +j 3.828427e+000
     y[ 1 ]=2.000000e+000 +j -1.000000e+000
     y[ 2 ]=3.000000e+000 +j -1.000000e+000
     y[ 3 ]=4.000000e+000 +j -1.828427e+000
     y[ 4 ]=5.000000e+000 +j -1.828427e+000
     y[ 5 ]=6.000000e+000 +j -1.000000e+000
     y[ 6 ]=7.000000e+000 +j -1.000000e+000
     y[ 7 ]=8.000000e+000 +j 3.828427e+000

    *************
     n=7
     m=

     m[0]= 1.000000e+000
     m[1]= 2.000000e+000
     m[2]= 3.000000e+000
     m[3]= 4.000000e+000
     m[4]= 5.000000e+000
     m[5]= 6.000000e+000
     m[6]= 7.000000e+000
     fft(m)[ 0 ]= 2.800000e+001 +j* 0.000000e+000
     fft(m)[ 1 ]= -3.500000e+000 +j* 7.267825e+000
     fft(m)[ 2 ]= -3.500000e+000 +j* 2.791157e+000
     fft(m)[ 3 ]= -3.500000e+000 +j* 7.988522e-001
     fft(m)[ 4 ]= -3.500000e+000 +j* -7.988522e-001
     fft(m)[ 5 ]= -3.500000e+000 +j* -2.791157e+000
     fft(m)[ 6 ]= -3.500000e+000 +j* -7.267825e+000
     n is odd
     twiddle [ 0 ]=2.800000e+001 +j 0.000000e+000
     twiddle [ 1 ]=-7.000000e+000 +j 1.453565e+001
     twiddle [ 2 ]=-7.000000e+000 +j 5.582314e+000
     twiddle [ 3 ]=-7.000000e+000 +j 1.597704e+000
     twiddle [ 4 ]=-0.000000e+000 +j -0.000000e+000
     twiddle [ 5 ]=-0.000000e+000 +j -0.000000e+000
     twiddle [ 6 ]=-0.000000e+000 +j -0.000000e+000
     y=hilbert(xr)
     y[ 0 ]=1.000000e+000 +j 3.102238e+000
     y[ 1 ]=2.000000e+000 +j -1.279048e+000
     y[ 2 ]=3.000000e+000 +j -7.974734e-001
     y[ 3 ]=4.000000e+000 +j -2.051434e+000
     y[ 4 ]=5.000000e+000 +j -7.974734e-001
     y[ 5 ]=6.000000e+000 +j -1.279048e+000
     y[ 6 ]=7.000000e+000 +j 3.102238e+000

    Friday, October 23, 2015 3:59 AM
  • //getarg.h

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

    double Arg(double xre, double xim )
    {
    if (xre>0) { return  atan(xim/xre); }     // 1,4 quadrant  0...+/-89.9999 deg (271-359.9999)
    if ((xre==0)&&(xim>0)) { return M_PI/2; } //   90 deg
    if ((xre==0)&&(xim<0)) { return -M_PI/2;} // -90 deg,    (+270 deg)
    if ((xre<0)&&(xim>=0)) { return atan(xim/xre)+M_PI ; }  //2 quadrant 90.0001...179.9999 deg
    if ((xre<0)&&(xim<0))  { return atan(xim/xre)-M_PI ; }  //3 quadrant 180.00001...269.9999 deg
    if ((xre==0)&&(xim=0)) {                  //|X|=0
                             printf("\n Arg parsing error ,(0+j*0) \n");
                             return 0  ;
                           }
    }

    //dftlib.h

     
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    #include "getarg.h"

    #define NMAX 1000

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


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

     
    CVector  dft(CVector xn  ,int N )
    {
       
    int j; 
    CVector 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;
       //dft_Xk->re[k]=dft_Xk->re[k]+dft_xn->re[n]*cos(angle);
       //dft_Xk->im[k]=dft_Xk->im[k]+dft_xn->re[n]*sin(angle);

       //  if angle with '-'
       //a=xnre
       //b=xnim
       //c=cos
       //d=sin
       //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
       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;
    }

    CVector DoScaleDFT(CVector 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;
    }

    CVector  idft(CVector idft_Xk, int N  )
    {

    int k,n ;
    double angle;
    CVector  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++)
      {
    //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
      angle=2*M_PI*k*n/N;
    //a=xkre
    //b=xkim
    //c=cos
    //d=sin


      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 ;
    }

    CVector pokecomplex(double *xnre, double *xnim, int num)
    {
    int i;
    CVector X;
    for(i=0;i<num;i++) 
    {
    X.re[i]=xnre[i];
    X.im[i]=xnim[i]; 
    }
     
    return X; 
    }

    CVector  get_r_fi(CVector x, int n, int mode   )
    {
     int i;
     CVector 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; 
    }

    CVector ScaleIQidft(CVector 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;

    }

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

    return res;
    }


    CVector GetVectorSum(CVector RFI , int N)
    {
    CVector S ;
    int i;
     for (i=0;i<N;i++)
     {
     S.re[i]=((RFI.re[i]*RFI.re[i]+RFI.im[i]*RFI.im[i]),0.5); //module
     S.im[i]=Arg(RFI.re[i],RFI.im[i]);                        // phase
     }
    return S;
    }

    CVector  get_fassvector( int n, double fs  )
    {
    CVector res;
    int k; 
    for (k=0;k<n;k++) 

    res.re[k]=fs*k/n;
    res.im[k]=0;
    }
    return res; 
    }


    void printcomplex(CVector xn , int n)
    {
     int i;
       printf("\n   " );
     for (i=0;i<n;i++)
     {
        printf("\n %d   %le + j *(%le) ",i,xn.re[i],xn.im[i]); 
        }    
        printf("\n   " );
    return;
    }

    void printarray(DArray xn , int n)
    {
     int i;
     printf("\n   " );
     for (i=0;i<n;i++)
     {
            printf("\n %d   %le  ",i,xn.data[i] ); 
            }    
            printf("\n " );
    return;
    }
    void writetofile(CVector xn , int n,  char *fname, char *message,int  mode      )
    {
    FILE *fp;
    int i;
     
     fp=fopen(fname,"a");
     
      fprintf(fp,"\n ; %s ;   " , message );
       fprintf(fp,"\n ; ; ; ;   " );
       if (mode==0) fprintf(fp,"\n n ; Re(x) ; +j* Im(x)  ;  " );
       if (mode==1) fprintf(fp,"\n n ; |X|   ; exp(j* Arg(x) );   " );
     for (i=0;i<n;i++)
     {
         fprintf(fp,"\n    %d ;  %le ;  %le ;   ",i,xn.re[i],xn.im[i] ); 
        }    
        fprintf(fp,"\n ; ; ; ;  " );
        fclose(fp);
        printf("\n data saved to file   %s    ", fname);
    return;
    }

    void writetofile_darray(DArray xn , int n,  char *fname, char *message     )
    {
    FILE *fp;
    int i;
     
           fp=fopen(fname,"a");
      fprintf(fp,"\n ; %s ;  ; " , message );
       fprintf(fp,"\n ; ;  " );
       fprintf(fp,"\n n ; Array ;   " );
       
     for (i=0;i<n;i++)
     {
            fprintf(fp,"\n    %d ;  %le ;    ",i,xn.data[i]  ); 
            }    
            fprintf(fp,"\n ; ;    " );
        fclose(fp);
        printf("\n data saved to file   %s    ", fname);
    return;
    }


    CVector hilbert(CVector xr , int n)
    {
     int  k,i,no2;
     double h[NMAX];
     CVector dft__xn,dft__Xk,idft__XK,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=floor (n/2);
     for (i=0;i<n;i++)
     {
      dft__xn.re[i]=xr.re[i];
      dft__xn.im[i]=0;
     }

    //fft
     dft__Xk=dft(dft__xn,n);
    //X={Xk.re,Xk.im};

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


     if((n>0)&&((2*floor(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);

    return idft__XN;
    }


    CVector getanalyticsignal(CVector xn, int num)
    {
    CVector 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.re[i];
     x_n.im[i]=Hilb.im[i];
    }
     
    return x_n;
    }

    DArray getenvelope(CVector x, int num)
    {
    int i;
    CVector  y;
    DArray 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;
    }

    Thursday, October 29, 2015 2:06 PM
  • //spectrum.h

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    //#include "dftlib.h"

    CVector GetSpectrumIQ(CVector xn, int num )
    {
    CVector xn1,Xk;
    xn1=getanalyticsignal(  xn,   num);
    Xk=dft(xn1, num);
    Xk=DoScaleDFT(Xk, num);

    return Xk;
    }


    CVector GetSpectrum(CVector xn, int num )
    {
    CVector xn1,Xk, result;
    xn1=getanalyticsignal(  xn,   num);
    Xk=dft(xn1, num);
     //printcomplex(Xk,num);
    Xk=DoScaleDFT(Xk, num);
    result=get_r_fi(Xk, num ,1 );
    return result;
    }

    Thursday, October 29, 2015 2:07 PM
  • //getsignal1.h

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    //#include "dftlib.h"


    //t = 0:1/fs:1;
    //1/10000
    //Ts=1/fs

    //hilb (sin(x)) =-cos(x)
    //hilb(cos(x)) = sin(x)
    //hilb (exp(it)=-i*exp(it)
    //hilb(exp(-it))=i*exp(-it)
    //x = 2.5cos(2*pi*20.3*t)+sin(2*pi*72.1*t)+cos(2*pi*100.1*t);
    //y = hilbert(x);
    //printf("\n xn=  "); 
     //xn=pokecomplex(m,n,num); 
     //xn=pokecomplex(inph,quad,num);
    //double inph[NMAX];
    //double quad[NMAX];
     
     //double m[NMAX]={1,2,3,4,5,6,7,8};
     //double n[NMAX]={0,0,0,0,0,0,0,0};


    CVector do_initsignal(double *fs, int *num)
    {

    double deg=M_PI/180;

    CVector res;

    int i;
    *fs=1000 ;
    *num=*fs  ; 
    int Nmax=*num; 
    //inph[i]=m[i];
     
     
     double ampl1=1;
    double freq1=10;
    double phase1=90;
    double ampl2=0.5;
    double freq2=72;
    double phase2=35;

     
    printf("\n Input A1 (0.05) ");
    scanf("%le",&ampl1);
    printf("\n Input f1 (10) ");
    scanf("%le",&freq1);
    printf("\n Input phase1 ,deg (90) ");
    scanf("%le",&phase1);
    printf("\n Input A2 (0.1)");
    scanf("%le",&ampl2);
    printf("\n Input f1 (72) ");
    scanf("%le",&freq2);
    printf("\n Input phase1 ,deg (35) ");
    scanf("%le",&phase2); 
     
     
     
     
     
    for (i=0;i<Nmax;i++)
    {


    double arg1= (2*M_PI*i*freq1/ *fs)+(phase1*deg);
    double arg2=(2*M_PI*i*freq2/  *fs)+(phase2*deg);
    res.re[i]=ampl1*cos( arg1 )+ampl2*cos( arg2);
     
    res.im[i]=0;
    }

    getch();
     
      
     return res;
    }


    • Edited by USERPC01 Thursday, October 29, 2015 2:08 PM
    Thursday, October 29, 2015 2:07 PM
  • //getsignal2.h

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    //#include "dftlib.h"

    CVector do_initsignal(double *fs, int *num)
    {

    double deg=M_PI/180;

    CVector res;

    int i;
    *fs=1000 ;
    *num=*fs  ; 
    int Nmax=*num; 

     
    double ampl1;
    double freq1;
    double phase1;
    double ampl2;
    double freq2;
    double phase2;


     
     
     
     
     
    for (i=0;i<Nmax;i++)
    {
     
    ampl1=1;
    ampl2=1; 
    freq1=5;
    freq2=100;

    double arg1=(2*M_PI*i*freq1/ *fs);
    double arg2=(2*M_PI*i*freq2/ *fs);
    res.re[i]=ampl2*(1+ampl1*cos(arg1))*cos(arg2);
     
    res.im[i]=0;
    }

    getch();
     
      
     return res;
    }

    Thursday, October 29, 2015 2:08 PM
  • //testspectrum.cpp

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
     #include "dftlib.h"
    #include "getsignal1.h"
    #include "spectrum.h"


    int main(void)
    {
    int num;
    CVector Xk,xn,xn1,hilb,hexp,fass;
    double fs;

      xn=do_initsignal(&fs, &num);
      writetofile( xn , num, "signal.csv", "xre(t) signal ",0 );
      fass=get_fassvector( num,  fs  );
      printf("\n f[n]=  ");
      //printcomplex(fass,num);
      writetofile(  fass , num, "fassoc.csv", "f[n], (Hz) ",1 );
     
      getch();

     Xk=GetSpectrum(  xn,   num );
     //printcomplex(Xk,num);
     writetofile( Xk , num, "dftas_exp1.csv", "dft ,  |X(jw)|,fi(jw)  ",1 );

     getch();

     Xk=GetSpectrumIQ(  xn,   num );
     //printcomplex(Xk,num);
     writetofile( Xk , num, "dftanalyt2.csv", "dft of  analytic (/num) signal  , Xi[k]  ,Xq[k] ",0 );

    return 0;

    }

    Thursday, October 29, 2015 2:09 PM
  • //#include <iostream>
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    #include "dftlib1.h"
    #include "getsignal1.h"

    int main(void)
    {
    double deg=M_PI/180;
    int num;
     double inph[NMAX];
    double quad[NMAX];
     
    int i,j;
    double *xkr,*xki;
    CVector Xk,xn,xn1,hilb,hexp,fass;
    double fs; 
       getch();
       printf ("\n init  xn   ");  
       xn=do_initsignal(&fs, &num);
       printf("\n xn=  ");
      // printcomplex(xn,num);
       writetofile( xn , num, "signal.csv", "xre(t) signal ",0 );
        
      fass=get_fassvector( num,  fs  );
      printf("\n f[n]=  ");
      writetofile(  fass , num, "fassoc.csv", "f[n], (Hz) ",1 );
      getch();
     
      //analytic signal
          
      xn1=getanalyticsignal(  xn,   num);
     
      writetofile( xn1 , num, "analytic.csv", "analytic signal I(t),Q(t)  ",1 );   
      getch(); 

      printf("\n Xk=dft(xn) ");
      Xk= dft(xn1, num);
      getch();

      Xk=DoScaleDFT(Xk, num);
     
      writetofile(Xk , num, "dftanalyt2.csv", "dft of  analytic (/num) signal  , Xi[k]  ,Xq[k] ",0 );
      getch();

     hexp=get_r_fi(Xk, num ,1 );
    // printcomplex( hexp , num);
     writetofile( hexp , num, "dftas_exp2.csv", "dft ,  |X(jw)|,fi(jw)  ",1 );


     
    return 0;
    }

    Thursday, October 29, 2015 2:09 PM
  • //testdft4.cpp

    //#include <iostream>
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    #include "dftlib.h"
    #include "getsignal1.h"

    int main(void)
    {
    double deg=M_PI/180;
    int num;
     double inph[NMAX];
    double quad[NMAX];
     
    int i,j;
    double *xkr,*xki;
    CVector Xk,xn,xn1,hilb,hexp,fass;
    double fs; 
       getch();
       printf ("\n init  xn   ");  
       xn=do_initsignal(&fs, &num);
       printf("\n xn=  ");
       printcomplex(xn,num);
       writetofile( xn , num, "signal.csv", "xre(t) signal ",0 );
       
     
      fass=get_fassvector( num,  fs  );
      printf("\n f[n]=  ");
      printcomplex(fass,num);
      writetofile(  fass , num, "fassoc.csv", "f[n], (Hz) ",1 );
     
      getch();
     
     
      printf("\n y=hilbert(xr) ");
      hilb=hilbert(xn,num);
      printcomplex( hilb , num);
      writetofile(  xn , num, "hilbert1.csv", "y=hilbert(xr(t)) ",0 );
      getch();
      
      printf("\n y=hilbert(xr) ,exp ");
      hexp=get_r_fi(hilb, num ,1 );
      printcomplex( hexp , num);
      writetofile(  xn , num, "hilbexp1.csv", "y=hilbert(xr(t)) exp ",0 );
    //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))
     
     //analytic signal
     
      for (i=0; i<num; i++)
    {
    inph[i]=xn.re[i] ;

    quad[i]=hilb.im[i];
    }
      
      xn1=pokecomplex(inph,quad,num);
      printcomplex(xn1 , num);
      writetofile( xn1 , num, "analytic.csv", "analytic signal I(t),Q(t)  ",1 );
     
     
      
     getch();
     
     
     printf("\n Xk=dft(xn) ");
     Xk= dft(xn1, num);
    // printcomplex( Xk , num);
    // writetofile(Xk , num, "dftanalyt1.csv", "dft of  analytic 1 signal  , Xi[k]  ,Xq[k] ",0 );
     getch();
    // hexp=get_r_fi(Xk, num ,1 );
     
    // printcomplex( hexp , num);
     
    // writetofile( hexp , num, "dftas_exp1.csv", "dft ,  |X(jw)|,fi(jw)  ",1 );
     
     getch();
     
     xn1=idft(Xk, num  );
     printcomplex( xn1 , num);
     writetofile( xn1 , num, "idft_dft.csv", "idft. I[n],Q[n]    ",0 );
     
      hexp=get_r_fi(hilb, num ,1 );
      printcomplex( hexp , num);
      writetofile(  xn , num, "idftdexp.csv", "y=idft(dft((z[n]))), R,fi ",1 );
     getch();
     
     
     
     
     for (i=0; i<num ; i++)
     {
      Xk.re[i]=Xk.re[i]/num;
      Xk.im[i]=Xk.im[i]/num;
     }
     printcomplex( Xk , num);
     writetofile(Xk , num, "dftanalyt2.csv", "dft of  analytic (/num) signal  , Xi[k]  ,Xq[k] ",0 );
     getch();

     hexp=get_r_fi(Xk, num ,1 );
     printcomplex( hexp , num);
     writetofile( hexp , num, "dftas_exp2.csv", "dft ,  |X(jw)|,fi(jw)  ",1 );
     
     
     printf("\n xn=idft(Xk), Xk scaled to 1/num  ");
     
      for (i=0;  i<num ; i++)
     {
      Xk.re[i]=Xk.re[i]*num;
      Xk.im[i]=Xk.im[i]*num;
     }
     
     xn1=idft(Xk, num  );
     printcomplex( xn1 , num);
     writetofile( xn1 , num, "idft_dft_sc.csv", "idft. I[n],Q[n] ,scaled 1/num   ",0 );
     
      hexp=get_r_fi(hilb, num ,1 );
      printcomplex( hexp , num);
      writetofile(  xn , num, "idftdexpsc.csv", "y=idft(dft((z[n]))), R,fi ",1 );
     
     getch();
     /*
     printf("\n y=hilbert(xr) ");
     hilb=hilbert(xn,num);
     printcomplex( hilb , num);
     hexp=get_r_fi(hilb, num  );
     printcomplex( hexp , num);
     */
     
    return 0;
    }

    Thursday, October 29, 2015 2:10 PM
  • //testdft2.cpp

     
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    #include "dftlib.h"


    int main(void)
    {
     
    int num=8;
     
    int i,j;
    double *xkr,*xki;
    CVector Xk,xn,xn1,hilb,hexp,fass;

     double m[NMAX]={1,2,3,4,5,6,7,8};
     double n[NMAX]={0,0,0,0,0,0,0,0};
     double inph[NMAX];
     double quad[NMAX];

    double fs=200 ;
    //num=fs  ; 
    //t = 0:1/fs:1;
    //1/10000
    //Ts=1/fs
     
    for (i=0;i<num;i++)
    {
    //x = 2.5cos(2*pi*20.3*t)+sin(2*pi*72.1*t)+cos(2*pi*100.1*t);
    //y = hilbert(x);

    //hilb (sin(x)) =-cos(x)
    //hilb(cos(x)) = sin(x)
    //hilb (exp(it)=-i*exp(it)
    //hilb(exp(-it))=i*exp(-it)
    inph[i]=m[i];

    //inph[i]= cos(2*M_PI*i*10 /fs+(90/360)*M_PI/180)+0.5*sin(2*M_PI*i*26/fs);
    //printf("\n inph(%d)=%le  ",i,inph[i]  );
    quad[i]=0;
    }
     
    getch();

    fass=get_fassvector( num,  fs  );
    printf("\n f[n]=  ");
    printcomplex(fass,num);
     writetofile(  fass , num, "fassoc.csv", "f[n], (Hz) ",1 ); 

     printf("\n xn=  "); 
     //xn=pokecomplex(m,n,num); 
      xn=pokecomplex(inph,quad,num); 
      printcomplex(xn,num);

     writetofile(  xn , num, "signal.csv", "xre(t) signal ",0 );
     getch();
     
     
     printf("\n y=hilbert(xr) ");
     hilb=hilbert(xn,num);
     printcomplex( hilb , num);
     writetofile(  xn , num, "hilbert1.csv", "y=hilbert(xr(t)) ",0 );
     getch();
     printf("\n y=hilbert(xr) ,exp ");
      hexp=get_r_fi(hilb, num ,1 );
      printcomplex( hexp , num);
      writetofile(  xn , num, "hilbexp1.csv", "y=hilbert(xr(t)) exp ",0 );
    //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))
     
     //analytic signal
     
      for (i=0; i<num; i++)
    {
    inph[i]=inph[i] ;

    quad[i]=hilb.im[i];
    }
      
      xn1=pokecomplex(inph,quad,num);
      printcomplex(xn1 , num);
      writetofile( xn1 , num, "analytic.csv", "analytic signal I(t),Q(t)  ",1 );
     
     
      
     getch();
     
     
     printf("\n Xk=dft(xn) ");
     Xk= dft(xn1, num);
     printcomplex( Xk , num);
     writetofile(Xk , num, "dftanalyt1.csv", "Xk=dft(xn) of  analytic  signal  , Xi[k]  ,Xq[k] ",0 );
     getch();

     hexp=get_r_fi(Xk, num ,1 );
     printcomplex( hexp , num);
     writetofile( hexp , num, "dftas_exp1.csv", "dft(xn)   |X(jw)|,fi(jw)  ",1 );
     
     getch();
     
      printf("\n xn=idft(Xk), Xk   ");
     xn1=idft(Xk, num  );
     printcomplex( xn1 , num);
     writetofile( xn1 , num, "idft_dft.csv", "xn1=idft.(Xk )  ",0 );
     hexp=get_r_fi(xn1, num ,1 );
      printcomplex( hexp , num);
      writetofile(   hexp , num, "idftdexp.csv", "xn1exp=idft(Xk), R,fi ",1 );
     getch();

      printf("\n Xk=Xk/N   ");
     for (i=0; i<num ; i++)
     {
      Xk.re[i]=Xk.re[i]/num;
      Xk.im[i]=Xk.im[i]/num;
     }
      printcomplex( Xk , num);
     writetofile(Xk , num, "dftanalyt2.csv", "dft (Xk/N) , Xi[k]  ,Xq[k] ",0 );
     getch();
     //srec=ifft(sp/n)
      printf("\n idft(Xk=Xk/N )  ");
     xn1=idft(Xk, num  );
     printcomplex( xn1 , num);
     writetofile( xn1 , num, "idft_dftsc.csv", "idft(Xk/N)     ",0 );
     //print srec
      hexp=get_r_fi(xn1, num ,1 );
      printcomplex( hexp , num);
      writetofile(  hexp , num, "idftdexpsc.csv", "y=idft(Xk/N), R,fi ",1 );


     getch();
    printf("\n scaling xn*N  ");
    //srec=srec*n
       for (i=0;  i<num ; i++)
     {
      xn1.re[i]=xn1.re[i]*num;
      xn1.im[i]=xn1.im[i]*num;
     }

    printcomplex( xn1 , num);
     writetofile( xn1 , num, "idftxknsc .csv", "xn1=idft(Xk/N)*N   ",0 );
     
      hexp=get_r_fi(hilb, num ,1 );
      printcomplex( hexp , num);
      writetofile(  hexp , num, "idftxkexpn.csv", "xnexp=idft(Xk/N)*N, R,fi ",1 );

    getch();

    return 0;
    }

    f[n]=

    0   0.000000e+000 + j *(0.000000e+000)
    1   2.500000e+001 + j *(0.000000e+000)
    2   5.000000e+001 + j *(0.000000e+000)
    3   7.500000e+001 + j *(0.000000e+000)
    4   1.000000e+002 + j *(0.000000e+000)
    5   1.250000e+002 + j *(0.000000e+000)
    6   1.500000e+002 + j *(0.000000e+000)
    7   1.750000e+002 + j *(0.000000e+000)

    data saved to file   fassoc.csv
    xn=

    0   1.000000e+000 + j *(0.000000e+000)
    1   2.000000e+000 + j *(0.000000e+000)
    2   3.000000e+000 + j *(0.000000e+000)
    3   4.000000e+000 + j *(0.000000e+000)
    4   5.000000e+000 + j *(0.000000e+000)
    5   6.000000e+000 + j *(0.000000e+000)
    6   7.000000e+000 + j *(0.000000e+000)
    7   8.000000e+000 + j *(0.000000e+000)

    data saved to file   signal.csv

     y=hilbert(xr)

     0   1.000000e+000 + j *(3.828427e+000)
     1   2.000000e+000 + j *(-1.000000e+000)
     2   3.000000e+000 + j *(-1.000000e+000)
     3   4.000000e+000 + j *(-1.828427e+000)
     4   5.000000e+000 + j *(-1.828427e+000)
     5   6.000000e+000 + j *(-1.000000e+000)
     6   7.000000e+000 + j *(-1.000000e+000)
     7   8.000000e+000 + j *(3.828427e+000)

     data saved to file   hilbert1.csv

     0   1.000000e+000 + j *(3.828427e+000)
     1   2.000000e+000 + j *(-1.000000e+000)
     2   3.000000e+000 + j *(-1.000000e+000)
     3   4.000000e+000 + j *(-1.828427e+000)
     4   5.000000e+000 + j *(-1.828427e+000)
     5   6.000000e+000 + j *(-1.000000e+000)
     6   7.000000e+000 + j *(-1.000000e+000)
     7   8.000000e+000 + j *(3.828427e+000)

     data saved to file   analytic.csv

     Xk=dft(xn)

     0   3.600000e+001 + j *(-8.881784e-016)
     1   -8.000000e+000 + j *(1.931371e+001)
     2   -8.000000e+000 + j *(8.000000e+000)
     3   -8.000000e+000 + j *(3.313708e+000)
     4   -4.000000e+000 + j *(-5.329071e-015)
     5   -6.245221e-015 + j *(5.324951e-015)
     6   -1.776357e-014 + j *(-8.095170e-015)
     7   1.536965e-014 + j *(2.146569e-014)

     data saved to file   dftanalyt1.csv

     0   3.600000e+001 + j *(-1.413580e-015)
     1   2.090501e+001 + j *(1.125000e+002)
     2   1.131371e+001 + j *(1.350000e+002)
     3   8.659138e+000 + j *(1.575000e+002)
     4   4.000000e+000 + j *(-1.800000e+002)
     5   8.207185e-015 + j *(1.395476e+002)
     6   1.952117e-014 + j *(-1.555004e+002)
     7   2.640079e-014 + j *(5.439695e+001)

     data saved to file   dftas_exp1.csv

     0   3.956874e+000 + j *(7.536119e+001)
     1   2.236068e+000 + j *(-2.656505e+001)
     2   3.162278e+000 + j *(-1.843495e+001)
     3   4.398084e+000 + j *(-2.456546e+001)
     4   5.323828e+000 + j *(-2.008673e+001)
     5   6.082763e+000 + j *(-9.462322e+000)
     6   7.071068e+000 + j *(-8.130102e+000)
     7   8.868870e+000 + j *(2.557360e+001)

     data saved to file   idftdexp.csv

     Xk=Xk/N

     0   4.500000e+000 + j *(-1.110223e-016)
     1   -1.000000e+000 + j *(2.414214e+000)
     2   -1.000000e+000 + j *(1.000000e+000)
     3   -1.000000e+000 + j *(4.142136e-001)
     4   -5.000000e-001 + j *(-6.661338e-016)
     5   -7.806527e-016 + j *(6.656188e-016)
     6   -2.220446e-015 + j *(-1.011896e-015)
     7   1.921206e-015 + j *(2.683211e-015)

     data saved to file   dftanalyt2.csv

     idft(Xk=Xk/N )

     0   1.250000e-001 + j *(4.785534e-001)
     1   2.500000e-001 + j *(-1.250000e-001)
     2   3.750000e-001 + j *(-1.250000e-001)
     3   5.000000e-001 + j *(-2.285534e-001)
     4   6.250000e-001 + j *(-2.285534e-001)
     5   7.500000e-001 + j *(-1.250000e-001)
     6   8.750000e-001 + j *(-1.250000e-001)
     7   1.000000e+000 + j *(4.785534e-001)

     data saved to file   idft_dftsc.csv

     0   4.946093e-001 + j *(7.536119e+001)
     1   2.795085e-001 + j *(-2.656505e+001)
     2   3.952847e-001 + j *(-1.843495e+001)
     3   5.497605e-001 + j *(-2.456546e+001)
     4   6.654785e-001 + j *(-2.008673e+001)
     5   7.603453e-001 + j *(-9.462322e+000)
     6   8.838835e-001 + j *(-8.130102e+000)
     7   1.108609e+000 + j *(2.557360e+001)

     data saved to file   idftdexpsc.csv

    idft_dftsc.csv

     ; ; ; ; 
     ; idft(Xk/N)      ;  
     ; ; ; ;  
     n ; Re(x) ; +j* Im(x)  ; 
        0 ;  1.250000e-001 ;  4.785534e-001 ;  
        1 ;  2.500000e-001 ;  -1.250000e-001 ;  
        2 ;  3.750000e-001 ;  -1.250000e-001 ;  
        3 ;  5.000000e-001 ;  -2.285534e-001 ;  
        4 ;  6.250000e-001 ;  -2.285534e-001 ;  
        5 ;  7.500000e-001 ;  -1.250000e-001 ;  
        6 ;  8.750000e-001 ;  -1.250000e-001 ;  
        7 ;  1.000000e+000 ;  4.785534e-001 ; 

    idftxkexpn.csv

     ; ; ; ; 
     ; xnexp=idft(Xk/N)*N, R,fi  ;  
     ; ; ; ;  
     n ; |X|   ; exp(j* Arg(x) );  
        0 ;  3.956874e+000 ;  7.536119e+001 ;  
        1 ;  2.236068e+000 ;  -2.656505e+001 ;  
        2 ;  3.162278e+000 ;  -1.843495e+001 ;  
        3 ;  4.398084e+000 ;  -2.456546e+001 ;  
        4 ;  5.323828e+000 ;  -2.008673e+001 ;  
        5 ;  6.082763e+000 ;  -9.462322e+000 ;  
        6 ;  7.071068e+000 ;  -8.130102e+000 ;  
        7 ;  8.868870e+000 ;  2.557360e+001 ;  
     ; ; ; ; 

    idftxknsc.csv
     ; ; ; ; 
     ; xn1=idft(Xk/N)*N    ;  
     ; ; ; ;  
     n ; Re(x) ; +j* Im(x)  ; 
        0 ;  1.000000e+000 ;  3.828427e+000 ;  
        1 ;  2.000000e+000 ;  -1.000000e+000 ;  
        2 ;  3.000000e+000 ;  -1.000000e+000 ;  
        3 ;  4.000000e+000 ;  -1.828427e+000 ;  
        4 ;  5.000000e+000 ;  -1.828427e+000 ;  
        5 ;  6.000000e+000 ;  -1.000000e+000 ;  
        6 ;  7.000000e+000 ;  -1.000000e+000 ;  
        7 ;  8.000000e+000 ;  3.828427e+000 ;  
     ; ; ; ; 

    • Edited by USERPC01 Thursday, October 29, 2015 9:55 PM
    Thursday, October 29, 2015 2:11 PM
  • // patched version of dftlib.h with correct rfi->iq, add synth procedures

     
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    //#include "getarg.h"
    //#include "cvector.h"

    #define NMAX 1000
     
     
    typedef struct  {
        double re[NMAX];
        double im[NMAX];
    } CVector;

    typedef struct  {
        double data[NMAX];
       
    } DArray;
     
     
    double Arg(double xre, double xim )
    {
    if (xre>0) { return  atan(xim/xre); }     // 1,4 quadrant  0...+/-89.9999 deg (271-359.9999)
    if ((xre==0)&&(xim>0)) { return M_PI/2; } //   90 deg
    if ((xre==0)&&(xim<0)) { return -M_PI/2;} // -90 deg,    (+270 deg)
    if ((xre<0)&&(xim>=0)) { return atan(xim/xre)+M_PI ; }  //2 quadrant 90.0001...179.9999 deg
    if ((xre<0)&&(xim<0))  { return atan(xim/xre)-M_PI ; }  //3 quadrant 180.00001...269.9999 deg
    if ((xre==0)&&(xim=0)) {                  //|X|=0
                             printf("\n Arg parsing error ,(0+j*0) \n");
                             return 0  ;
                           }
    }

    CVector  dft(CVector xn  ,int N )
    {
       
    int j; 
    CVector 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;
       //dft_Xk->re[k]=dft_Xk->re[k]+dft_xn->re[n]*cos(angle);
       //dft_Xk->im[k]=dft_Xk->im[k]+dft_xn->re[n]*sin(angle);

       //  if angle with '-'
       //a=xnre
       //b=xnim
       //c=cos
       //d=sin
       //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
       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;
    }

    CVector DoScaleDFT(CVector 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;
    }

    CVector  idft(CVector idft_Xk, int N  )
    {

    int k,n ;
    double angle;
    CVector  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++)
      {
    //(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
      angle=2*M_PI*k*n/N;
    //a=xkre
    //b=xkim
    //c=cos
    //d=sin


      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 ;
    }

    CVector pokecomplex(double *xnre, double *xnim, int num)
    {
    int i;
    CVector X;
    for(i=0;i<num;i++) 
    {
    X.re[i]=xnre[i];
    X.im[i]=xnim[i]; 
    }
     
    return X; 
    }

    CVector  get_r_fi(CVector x, int n, int mode   )
    {
     int i;
     CVector 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; 
    }

    CVector ScaleIQidft(CVector 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;

    }

    CVector GetIQfromR_fi(CVector RFI, int numb, int mode )
    {
    int i;
    CVector res;
    double rad=1;
    if (mode==1) rad=M_PI/180;
     for (i=0;i<numb;i++)
     {
      //10   ; 1.000000e+000 ;  9.000000e+001 ;
     res.re[i]= RFI.re[i]*cos(rad* RFI.im[i] );
     res.im[i]= RFI.re[i]*sin(rad* RFI.im[i] );
     }

    return res;
    }

    CVector GetVectorSum(CVector RFI , int N)
    {
    CVector S ;
    int i;
     for (i=0;i<N;i++)
     {
     S.re[i]=((RFI.re[i]*RFI.re[i]+RFI.im[i]*RFI.im[i]),0.5); //module
     S.im[i]=Arg(RFI.re[i],RFI.im[i]);                        // phase
     }
    return S;
    }

    CVector  get_tassvector( int n, double fs  )
    {
    CVector res;
    int k; 
    for (k=0;k<n;k++) 

    res.re[k]=k/fs;
    res.im[k]=0;
    }
    return res; 
    }

    CVector  get_fassvector( int n, double fs  )
    {
    CVector res;
    int k; 
    for (k=0;k<n;k++) 

    res.re[k]=fs*k/n;
    res.im[k]=0;
    }
    return res; 
    }


    void printcomplex(CVector xn , int n)
    {
     int i;
       printf("\n   " );
     for (i=0;i<n;i++)
     {
        printf("\n %d   %le + j *(%le) ",i,xn.re[i],xn.im[i]); 
        }    
        printf("\n   " );
    return;
    }

    void printarray(DArray xn , int n)
    {
     int i;
     printf("\n   " );
     for (i=0;i<n;i++)
     {
            printf("\n %d   %le  ",i,xn.data[i] ); 
            }    
            printf("\n " );
    return;
    }
    void writetofile(CVector xn , int n,  char *fname, char *message,int  mode      )
    {
    FILE *fp;
    int i;
     
     fp=fopen(fname,"a");
     
      fprintf(fp,"\n ; %s ;   " , message );
       fprintf(fp,"\n ; ; ; ;   " );
       if (mode==0) fprintf(fp,"\n n ; Re(x) ; +j* Im(x)  ;  " );
       if (mode==1) fprintf(fp,"\n n ; |X|   ; exp(j* Arg(x) );   " );
     for (i=0;i<n;i++)
     {
         fprintf(fp,"\n    %d ;  %le ;  %le ;   ",i,xn.re[i],xn.im[i] ); 
        }    
        fprintf(fp,"\n ; ; ; ;  " );
        fclose(fp);
        printf("\n data saved to file   %s    ", fname);
    return;
    }

    void writetofile_darray(DArray xn , int n,  char *fname, char *message     )
    {
    FILE *fp;
    int i;
     
           fp=fopen(fname,"a");
      fprintf(fp,"\n ; %s ;  ; " , message );
       fprintf(fp,"\n ; ;  " );
       fprintf(fp,"\n n ; Array ;   " );
       
     for (i=0;i<n;i++)
     {
            fprintf(fp,"\n    %d ;  %le ;    ",i,xn.data[i]  ); 
            }    
            fprintf(fp,"\n ; ;    " );
        fclose(fp);
        printf("\n data saved to file   %s    ", fname);
    return;
    }


    CVector hilbert(CVector xr , int n)
    {
     int  k,i,no2;
     double h[NMAX];
     CVector dft__xn,dft__Xk,idft__XK,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=floor (n/2);
     for (i=0;i<n;i++)
     {
      dft__xn.re[i]=xr.re[i];
      dft__xn.im[i]=0;
     }

    //fft
     dft__Xk=dft(dft__xn,n);
    //X={Xk.re,Xk.im};

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


     if((n>0)&&((2*floor(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);

    return idft__XN;
    }


    CVector getanalyticsignal(CVector xn, int num)
    {
    CVector 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.re[i];
     x_n.im[i]=Hilb.im[i];
    }
     
    return x_n;
    }

    DArray getenvelope(CVector x, int num)
    {
    int i;
    CVector  y;
    DArray 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;
    }
    CVector GetSpectrumIQ(CVector xn, int num )
    {
    CVector xn1,Xk;
    xn1=getanalyticsignal(  xn,   num);
    Xk=dft(xn1, num);
    Xk=DoScaleDFT(Xk, num);

    return Xk;
    }


    CVector GetSpectrum(CVector xn, int num )
    {
    CVector xn1,Xk, result;
    xn1=getanalyticsignal(  xn,   num);
    Xk=dft(xn1, num);
     //printcomplex(Xk,num);
    Xk=DoScaleDFT(Xk, num);
    result=get_r_fi(Xk, num ,1 );
    return result;
    }

    CVector SynthSignalfromIQ(CVector idftXk, int Num)
    {

    CVector idftxn;
    idftXk=ScaleIQidft(idftXk, Num);
    idftxn=idft(  idftXk,  Num  );
    return idftxn;

    }
     
     CVector SynthSignalfromRFi(CVector idftXk, int Num, int mode)
    {

    CVector idftxn;

    idftXk=GetIQfromR_fi(idftXk, Num, mode);
    idftxn=SynthSignalfromIQ(idftXk,  Num);
    return idftxn;

    }


    • Edited by USERPC01 Thursday, October 29, 2015 10:28 PM
    Thursday, October 29, 2015 4:28 PM
  • //testspectrum.cpp

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
     #include "dftlib.h"
    #include "getsignal1.h"

    int main(void)
    {
    int num;
    CVector Xk,xn,xn1,hilb,hexp,fass;
    double fs;

      xn=do_initsignal(&fs, &num);
      writetofile( xn , num, "signal.csv", "xre(t) signal ",0 );
      fass=get_fassvector( num,  fs  );
      printf("\n f[n]=  ");
      //printcomplex(fass,num);
      writetofile(  fass , num, "fassoc.csv", "f[n], (Hz) ",1 );
     
      getch();

     Xk=GetSpectrum(  xn,   num );
     //printcomplex(Xk,num);
     writetofile( Xk , num, "dftas_exp1.csv", "dft ,  |X(jw)|,fi(jw)  ",1 );

     getch();


     Xk=GetSpectrumIQ(  xn,   num );
     //printcomplex(Xk,num);
     writetofile( Xk , num, "dftanalyt2.csv", " analytic (/num) signal  , Xi[k]  ,Xq[k] ",0 );

    return 0;

    }

    Thursday, October 29, 2015 4:29 PM
  • //testsynth.cpp

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    //#include "cvector.h"
    #include "dftlib.h"
    #include "getsignal1.h"
    //#include "spectrum.h"
    //#include "synthsig.h"
     
    int main(void)
    {
    int num;
    CVector Xk,xn,xn1,hilb,hexp,fass; 
    double fs;

      xn=do_initsignal(&fs, &num);
      writetofile( xn , num, "signal.csv", "xre(t) signal ",0 );
      fass=get_fassvector( num,  fs  );
      printf("\n f[n]=  ");
      //printcomplex(fass,num);
      writetofile(  fass , num, "fassoc.csv", "f[n], (Hz) ",1 );
     
      getch();

      Xk=GetSpectrumIQ(  xn,   num );
      //printcomplex(Xk,num);
      writetofile( Xk , num, "dftanalyt2.csv", "dft of  analytic (/num) signal  , Xi[k]  ,Xq[k] ",0 );

     
       xn1=SynthSignalfromIQ(Xk, num);
       writetofile( xn1 , num, "signalidftiq.csv", "x (t)  I, Q ",1 );

    Xk=get_r_fi(Xk, num, 1  );
     
    xn1=SynthSignalfromRFi(Xk, num,1);
      
     
       writetofile( xn1 , num, "signalidft.csv", "xre(t) signal ",0 );
    return 0;

    }

    Thursday, October 29, 2015 4:29 PM
  • //getenvel.cpp

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    #include "dftlib.h"
    #include "getsignal2.h"

    int main(void)
    {
    double deg=M_PI/180;
    int num;

     
    int i,j;
    double *xkr,*xki;
    CVector xn,xn1,hilb,hexp,fass;
    DArray  X;
    double fs; 
       getch();
       printf ("\n init  xn   ");  
       xn=do_initsignal(&fs, &num);
       printf("\n xn=  ");
      // printcomplex(xn,num);
       writetofile( xn , num, "signal.csv", "xre(t) signal ",0 );
        
      fass=get_fassvector( num,  fs  );
      printf("\n f[n]=  ");
      writetofile(  fass , num, "fassoc.csv", "f[n], (Hz) ",1 );
      getch();
     
      //analytic signal
          
      xn1=getanalyticsignal(  xn,   num);
     
      writetofile( xn1 , num, "analytic.csv", "analytic signal I(t),Q(t)  ",1 );   
      getch(); 

      printf("\n y=envelope(xn) ");
      X= getenvelope( xn,  num);
     //  printarray(X,num);
       writetofile_darray( X , num, "envelope.csv", "analytic signal      "  );
      getch();

     

    return 0;
    }

    Thursday, October 29, 2015 4:31 PM
  • //getsignal2.h

    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <conio.h>
    //#include "dftlib.h"

    CVector do_initsignal(double *fs, int *num)
    {

    double deg=M_PI/180;

    CVector res;

    int i;
    *fs=1000 ;
    *num=*fs  ; 
    int Nmax=*num; 

     
    double ampl1;
    double freq1;
    double phase1;
    double ampl2;
    double freq2;
    double phase2;


     
     
     
     
     
    for (i=0;i<Nmax;i++)
    {
     
    ampl1=1;
    ampl2=1; 
    freq1=5;
    freq2=100;

    double arg1=(2*M_PI*i*freq1/ *fs);
    double arg2=(2*M_PI*i*freq2/ *fs);
    res.re[i]=ampl2*(1+ampl1*cos(arg1))*cos(arg2);
     
    res.im[i]=0;
    }

    getch();
     
      
     return res;
    }

    Thursday, October 29, 2015 4:31 PM
  • http://www.mathworks.com/help/signal/ug/analytic-signal-for-cosine.html

    https://en.wikipedia.org/wiki/Hilbert_transform

    VBA   script for Excel (ALT+F11) 

    'Option Explicit
     Const M_PI As Double = 3.141592654
     Const NMAX As Integer = 32767

     Sub Getfassoc(ByVal fs As Double, ByRef fvect() As Double, ByVal Numb As Integer)
       Dim k As Integer
        For k = 0 To Numb - 1
         fvect(k + 1) = fs * k / Numb
        Next k
      End Sub
     
     
      Sub Gettassoc(ByVal fs As Double, ByRef tvect() As Double, ByVal Numb As Integer)
       Dim n As Integer
        For n = 0 To Numb - 1
         tvect(n + 1) = n / fs
        Next n
      End Sub


      Function Arg(ByVal x_re As Double, ByVal x_im As Double) As Double
       'Dim x_re,x_im As Double
       'Dim arg As Double
       Dim message As String
      If (x_re > 0) Then
                 Arg = Atn(x_im / x_re)
                 End If
      If (x_re = 0) And (x_im > 0) Then
                 Arg = M_PI / 2
                 End If
      If (x_re = 0) And (x_im < 0) Then
                 Arg = -M_PI / 2
                 End If
      If (x_re < 0) And (x_im >= 0) Then
                 Arg = Atn(x_im / x_re) + M_PI
                 End If
      If (x_re < 0) And (x_im < 0) Then
                 Arg = Atn(x_im / x_re) - M_PI
                 End If
      If (x_re = 0) And (x_im = 0) Then
               '  MsgBox (" Xre=0 and Xim=0, Invalid result")
                 Arg = 0
                 End If

     End Function
     '36.  - 4. + 9.6568542i  - 4. + 4.i  - 4. + 1.6568542i  - 4.  - 4. - 1.6568542i
      ' - 4. - 4.i  - 4. - 9.6568542i

     Sub GetIQfromR_fi(ByVal Num As Integer, ByRef modXk() As Double, ByRef argXk() As Double, _
     ByRef Xkre() As Double, ByRef Xkim() As Double, ByVal mode As Integer)
       Dim i As Integer
       Dim rad As Double
      
        rad = 1
         If mode = 1 Then
           rad = M_PI / 180
         End If
        
         For i = 1 To Num
          Xkre(i) = modXk(i) * Cos(rad * argXk(i))
          Xkim(i) = modXk(i) * Sin(rad * argXk(i))
         Next i
        
       End Sub

     


     Sub DFT(ByVal numdft As Integer, ByRef dft_xnre() As Double, ByRef dft_xnim() As Double, ByRef dft_Xkre() As Double, ByRef dft_Xkim() As Double)

      Dim i, k, n As Integer
     
     Dim angle As Double

     For k = 0 To numdft - 1
      dft_Xkre(k + 1) = 0
      dft_Xkim(k + 1) = 0
       For n = 0 To numdft - 1
        angle = -2 * M_PI * k * n / numdft
       dft_Xkre(k + 1) = dft_Xkre(k + 1) + dft_xnre(n + 1) * Cos(angle) - dft_xnim(n + 1) * Sin(angle)
       dft_Xkim(k + 1) = dft_Xkim(k + 1) + dft_xnre(n + 1) * Sin(angle) + dft_xnim(n + 1) * Cos(angle)
        Next n
      Next k

    End Sub
     
     Sub ScaleDFT(ByVal Num As Integer, ByRef ScXkre() As Double, ByRef ScXkim() As Double)

    Dim i  As Integer
      For i = 1 To Num
         ScXkre(i) = ScXkre(i) / Num
         ScXkim(i) = ScXkim(i) / Num
      Next i
     End Sub
     Sub ScaleIQidft(ByVal Num As Integer, ByRef ScXk_re() As Double, ByRef ScXk_im() As Double)
       Dim i  As Integer
       For i = 1 To Num
        ScXk_re(i) = ScXk_re(i) * Num
        ScXk_im(i) = ScXk_im(i) * Num
       Next i
      End Sub

     Sub IDFT(ByVal numidft As Integer, ByRef idft_Xkre() As Double, _
     ByRef idft_Xkim() As Double, ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
     Dim k, n    As Integer
      For n = 0 To numidft - 1
       'idft_item(n) = numidft - 1
       idft_xnre(n + 1) = 0
       idft_xnim(n + 1) = 0
        For k = 0 To numidft - 1
        angle = 2 * M_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)
        Next k
       
       idft_xnre(n + 1) = idft_xnre(n + 1) / numidft
       idft_xnim(n + 1) = idft_xnim(n + 1) / numidft

     Next n
     
      'idft_xnre(n), idft_xnim(n), numidft
     End Sub

     

     Sub Hilbert(ByRef xre() As Double, ByVal Numb As Integer, ByRef Hilbertre() As Double, ByRef Hilbertim() As Double)
      'MATLAB algorithm
       Dim k, i, no2 As Integer
       Dim h(1 To NMAX) As Double
     
      ' Dim dft_xnre(1 To NMAX)     As Double
       Dim dft_xnim(1 To NMAX)     As Double
       Dim dft_Xkre(1 To NMAX)    As Double
       Dim dft_Xkim(1 To NMAX)    As Double
       Dim idft_Xkre(1 To NMAX)   As Double
       Dim idft_Xkim(1 To NMAX)   As Double
       Dim idft_xnre(1 To NMAX)   As Double
       Dim idft_xnim(1 To NMAX)   As Double

      If (Numb = 0) Then
         Hilbertre(1) = 0
         Hilbertim(1) = 0
       MsgBox (" Hilbert() :   Parsing error ")
        End If


      no2 = Int(Numb / 2)
     
      
      
      
       For i = 1 To Numb
     
          dft_xnim(i) = 0
       Next i

       Call DFT(Numb, xre, dft_xnim, dft_Xkre, dft_Xkim)

       i = 0
         Do While i <= Numb

            h(i + 1) = 0
             i = i + 1
        Loop

      'if  N is even
     
        If (Numb > 0) And (2 * no2) = Numb Then
            h(1) = 1
            h(no2 + 1) = 1
            For k = 2 To no2
               h(k) = 2
            Next k

       'if N is odd

        ElseIf (Numb > 0) Then

             h(1) = 1
              For k = 2 To (0.5 * (Numb + 1))
                h(k) = 2
              Next k
          
         End If


         For i = 1 To Numb
          idft_Xkre(i) = dft_Xkre(i) * h(i)
          idft_Xkim(i) = dft_Xkim(i) * h(i)
         Next i
     
        
         Call IDFT(Numb, idft_Xkre, idft_Xkim, idft_xnre, idft_xnim)
       
       For i = 1 To Numb
       Hilbertre(i) = idft_xnre(i)
       Hilbertim(i) = idft_xnim(i)
       Next i
     
      End Sub

     Sub GetAnalyticSignal(ByRef xnre() As Double, ByVal Numb As Integer, ByRef x_re() As Double, ByRef x_im() As Double)

      'I=x(t)
       'Q=jHilb(t)
       'Hilb=Hilbert()
       Dim i As Integer
       Dim Hilbertre(1 To NMAX) As Double
       Dim Hilbertim(1 To NMAX) As Double
       Call Hilbert(xnre, Numb, Hilbertre, Hilbertim)
       For i = 1 To Numb
        x_re(i) = xnre(i)
        x_im(i) = Hilbertim(i)
       Next i
      End Sub

     Sub getenvelope(ByRef xnre() As Double, ByVal Numb As Integer, ByRef envelope() As Double)
       Dim i As Integer
       Dim Hilbertre(1 To NMAX) As Double
       Dim Hilbertim(1 To NMAX) As Double
       Call Hilbert(xnre, Numb, Hilbertre, Hilbertim)
       For i = 1 To Numb
       envelope(i) = ((Hilbertre(i) ^ 2) + (Hilbertim(i) ^ 2)) ^ 0.5
       Next i
      End Sub


     Sub GetSpectrumIQ(ByRef xn() As Double, ByVal Num As Integer, ByRef Xkre() As Double, ByRef Xkim() As Double)
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double

      Call GetAnalyticSignal(xn, Num, x_re, x_im)
      Call DFT(Num, x_re, x_im, Xkre, Xkim)
       Call ScaleDFT(Num, Xkre, Xkim)

     End Sub
     Sub GetR_fi(ByVal Numb As Integer, ByRef Xkre() As Double, ByRef Xkim() As Double, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Integer)
       Dim i As Integer
       For i = 1 To Numb
        modXk(i) = (Xkre(i) ^ 2 + Xkim(i) ^ 2) ^ 0.5
        argXk(i) = Arg(Xkre(i), Xkim(i))

      If mode = 1 Then
        argXk(i) = argXk(i) * 180 / M_PI
       End If

      Next i

     End Sub


    Sub GetSpectrum(ByRef xn() As Double, ByVal Num As Integer, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Integer)
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double
      Dim Xkre(1 To NMAX) As Double
      Dim Xkim(1 To NMAX) As Double
          Call GetAnalyticSignal(xn, Num, x_re, x_im)
          Call DFT(Num, x_re, x_im, Xkre, Xkim)
          Call ScaleDFT(Num, Xkre, Xkim)
          Call GetR_fi(Num, Xkre, Xkim, modXk, argXk, mode)
     
      End Sub
     
     
      Sub GetSpectrumfromIQ(ByRef x_re() As Double, ByRef x_im() As Double, _
      ByVal Num As Integer, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Integer)
      
      
      Dim Xkre(1 To NMAX) As Double
      Dim Xkim(1 To NMAX) As Double
          
          Call DFT(Num, x_re, x_im, Xkre, Xkim)
          Call ScaleDFT(Num, Xkre, Xkim)
          Call GetR_fi(Num, Xkre, Xkim, modXk, argXk, mode)
     
      End Sub


       Sub SynthSignalfromIQ(ByVal Num As Integer, ByRef Xkre() As Double, ByRef Xkim() As Double, _
       ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
       
               
              'Call ScaleIQidft(Num, Xkre, Xkim)
              Call IDFT(Num, Xkre, Xkim, idft_xnre, idft_xnim)
              Call ScaleIQidft(Num, idft_xnre, idft_xnim)
       End Sub

       Sub SynthSignalfromRFi(ByRef modXk() As Double, ByRef argXk() As Double, _
        ByRef xnre() As Double, ByRef xnim() As Double, ByVal Num As Integer, ByVal mode As Integer)
        Dim Xk_re(1 To NMAX) As Double
        Dim Xk_im(1 To NMAX) As Double
         Call GetIQfromR_fi(Num, modXk, argXk, Xk_re, Xk_im, mode)
         Call SynthSignalfromIQ(Num, Xk_re, Xk_im, xnre, xnim)
       End Sub

    Sub printcomplex(ByRef xnre() As Double, ByRef xnim() As Double, _
     ByVal Num As Integer, ByVal Col1 As Integer, ByVal Col2 As Integer)
         For i = 1 To Num
              Cells(i + 1, Col1).Value = xnre(i)
              Cells(i + 1, Col2).Value = xnim(i)
            Next i
     End Sub

    Sub printarray(ByRef xnre() As Double, ByVal Num As Integer, ByVal Col1 As Integer)
         For i = 1 To Num
            Cells(i + 1, Col1).Value = xnre(i)
            
            Next i
     End Sub


     Sub writetofile(ByRef xnre() As Double, ByRef xnim() As Double, _
     ByRef FilePath As String, ByVal n As Integer)
      Dim FilePath As String
      Dim i As Integer

     Open FilePath For Output As #1
     
         'Write #1,message
         'Write #1,"\n ; ; ; ;   "
         ' If  mode=0  Then
             'Write #1,"\n n ; Re(x) ; +j* Im(x)      ;   "
             '  End If
         
          'If  mode=1 Then
             ' Write #1,"\n n ; |X|   ; exp(j* Arg(x) );   "
             '  End If

        For i = 1 To n
             Write #1, i
             Write #1, xnre(i)
             Write #1, xnim(i)
         Next i
            
         Write #1, "\n ; ; ; ;  "
         Close #1
         MsgBox (" data saved to file")
       
     End Sub

     Sub TestHilbert()
    Dim fs As Double
     Dim fv(1 To NMAX) As Double
     Dim tv(1 To NMAX) As Double
     Dim Xi(1 To NMAX) As Double
     Dim Xq(1 To NMAX) As Double

     Dim Yi(1 To NMAX) As Double
     Dim Yq(1 To NMAX) As Double
     Dim Hi(1 To NMAX) As Double
     Dim Hq(1 To NMAX) As Double
     Dim Ii(1 To NMAX) As Double
     Dim Iq(1 To NMAX) As Double
     Dim modX(1 To NMAX) As Double
     Dim argX(1 To NMAX) As Double
     Dim env(1 To NMAX) As Double
     Dim modX1(1 To NMAX) As Double
     Dim argX1(1 To NMAX) As Double
     Dim Xi1(1 To NMAX) As Double
     Dim i As Integer
     Dim n As Integer
    Dim x_re(1 To NMAX) As Double
    Dim Yre(1 To NMAX) As Double
    Dim Yim(1 To NMAX) As Double
    n = 8
     Xi(1) = 1
     Xi(2) = 2
     Xi(3) = 3
     Xi(4) = 4
     Xi(5) = 5
     Xi(6) = 6
     Xi(7) = 7
     Xi(8) = 8


     Xq(1) = 0
     Xq(2) = 0
     Xq(3) = 0
     Xq(4) = 0
     Xq(5) = 0
     Xq(6) = 0
     Xq(7) = 0
     Xq(8) = 0
     fs = 1000
      Call printcomplex(Xi, Xq, n, 1, 2)
      Call Getfassoc(fs, fv, n)
      Call Gettassoc(fs, tv, n)
      Call printcomplex(tv, fv, n, 3, 4)
     
     Call DFT(n, Xi, Xq, Yi, Yq)
     
     Call printcomplex(Yi, Yq, n, 5, 6)
     
       Call IDFT(n, Yi, Yq, Ii, Iq)
        Call printcomplex(Ii, Iq, n, 7, 8)
       Call Hilbert(Xi, n, Hi, Hq)
        Call printcomplex(Hi, Hq, n, 9, 10)
        Call GetSpectrum(Xi, n, modX, argX, 1)
         Call printcomplex(modX, argX, n, 11, 12)
      
     
       Call getenvelope(Xi, n, env)
       Call printarray(env, n, 13)
     
     
      
       n = 16
     
     
       For i = 0 To n - 1
        argX(i + 1) = 0
        modX(i + 1) = 0
        Next i
      
       modX(2) = 1
       argX(2) = 0
       modX(6) = 0
       argX(6) = 0
      
      '  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
      
       Call SynthSignalfromRFi(modX, argX, Xi, Xq, n, 1)
       Call printcomplex(Xi, Xq, n, 14, 15)
      
      
          Call Getfassoc(fs, fv, n)
          Call printarray(fv, n, 16)
         
         Call GetSpectrumfromIQ(Xi, Xq, n, modX1, argX1, 1)
         Call printcomplex(modX1, argX1, n, 17, 18)
         
         For i = 0 To n - 1
          Xi1(i + 1) = 0.5 * Cos(2 * M_PI * 250 * i / fs + 0) + 1 * Sin(2 * M_PI * 125 * i / fs + 0)
      
         Next i
          'Xi1(1) = 0
         
         Call GetSpectrum(Xi1, n, modX1, argX1, 1)
         Call printcomplex(modX1, argX1, n, 19, 20)
     
         Call GetSpectrumIQ(Xi1, n, Yre, Yim)
         Call printcomplex(Yre, Yim, n, 21, 22)
     
     End Sub



    • Edited by USERPC01 Friday, November 13, 2015 7:24 PM
    Thursday, November 12, 2015 4:21 PM
  • Version for  65536  bins

    Option Explicit
     Const M_PI As Double = 3.14159265358979
      'Const NMAX As Integer = 32767
      'Const NMAX As Long = 2147483647
      Const NMAX As Long = 131072
    'Const NMAX As Double = 1.79769313486231570E+308

     Sub Getfassoc(ByVal fs As Double, ByRef fvect() As Double, ByVal Numb As Long)
       Dim k As Long
        For k = 0 To Numb - 1
         fvect(k + 1) = fs * k / Numb
        Next k
      End Sub
     
     
      Sub Gettassoc(ByVal fs As Double, ByRef tvect() As Double, ByVal Numb As Long)
       Dim N As Long
        For N = 0 To Numb - 1
         tvect(N + 1) = N / fs
        Next N
      End Sub


      Function Arg(ByVal x_re As Double, ByVal x_im As Double) As Double
       'Dim x_re,x_im As Double
       'Dim arg As Double
       Dim message As String
      If (x_re > 0) Then
                 Arg = Atn(x_im / x_re)
                 End If
      If (x_re = 0) And (x_im > 0) Then
                 Arg = M_PI / 2
                 End If
      If (x_re = 0) And (x_im < 0) Then
                 Arg = -M_PI / 2
                 End If
      If (x_re < 0) And (x_im >= 0) Then
                 Arg = Atn(x_im / x_re) + M_PI
                 End If
      If (x_re < 0) And (x_im < 0) Then
                 Arg = Atn(x_im / x_re) - M_PI
                 End If
      If (x_re = 0) And (x_im = 0) Then
               '  MsgBox (" Xre=0 and Xim=0, Invalid result")
                 Arg = 0
                 End If

     End Function
     '36.  - 4. + 9.6568542i  - 4. + 4.i  - 4. + 1.6568542i  - 4.  - 4. - 1.6568542i
     
     '         column 7 to 8
     
      ' - 4. - 4.i  - 4. - 9.6568542i
     Sub GetIQfromR_fi(ByVal Num As Long, ByRef modXk() As Double, ByRef argXk() As Double, _
     ByRef Xkre() As Double, ByRef Xkim() As Double, ByVal mode As Long)
       Dim i As Long
       Dim rad As Double
      
        rad = 1
         If mode = 1 Then
           rad = M_PI / 180
         End If
        
         For i = 1 To Num
          Xkre(i) = modXk(i) * Cos(rad * argXk(i))
          Xkim(i) = modXk(i) * Sin(rad * argXk(i))
         Next i
        
       End Sub

     


     Sub DFT(ByVal numdft As Long, ByRef dft_xnre() As Double, ByRef dft_xnim() As Double, ByRef dft_Xkre() As Double, ByRef dft_Xkim() As Double)

      Dim i, k, N As Long
     
     Dim angle As Double

     For k = 0 To numdft - 1
      dft_Xkre(k + 1) = 0
      dft_Xkim(k + 1) = 0
       For N = 0 To numdft - 1
        angle = -2 * M_PI * k * N / numdft
       dft_Xkre(k + 1) = dft_Xkre(k + 1) + dft_xnre(N + 1) * Cos(angle) - dft_xnim(N + 1) * Sin(angle)
       dft_Xkim(k + 1) = dft_Xkim(k + 1) + dft_xnre(N + 1) * Sin(angle) + dft_xnim(N + 1) * Cos(angle)
        Next N
      Next k

    End Sub
     
     Sub ScaleDFT(ByVal Num As Long, ByRef ScXkre() As Double, ByRef ScXkim() As Double)

    Dim i  As Long
      For i = 1 To Num
         ScXkre(i) = ScXkre(i) / Num
         ScXkim(i) = ScXkim(i) / Num
      Next i
     End Sub
     Sub ScaleIQidft(ByVal Num As Long, ByRef ScXk_re() As Double, ByRef ScXk_im() As Double)
       Dim i  As Long
       For i = 1 To Num
        ScXk_re(i) = ScXk_re(i) * Num
        ScXk_im(i) = ScXk_im(i) * Num
       Next i
      End Sub

     Sub IDFT(ByVal numidft As Long, ByRef idft_Xkre() As Double, _
     ByRef idft_Xkim() As Double, ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
     Dim k, N    As Long
      For N = 0 To numidft - 1
       'idft_item(n) = numidft - 1
       idft_xnre(N + 1) = 0
       idft_xnim(N + 1) = 0
        For k = 0 To numidft - 1
        angle = 2 * M_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)
        Next k
       
       idft_xnre(N + 1) = idft_xnre(N + 1) / numidft
       idft_xnim(N + 1) = idft_xnim(N + 1) / numidft

     Next N
     
      'idft_xnre(n), idft_xnim(n), numidft
     End Sub

     

     Sub Hilbert(ByRef xre() As Double, ByVal Numb As Long, ByRef Hilbertre() As Double, ByRef Hilbertim() As Double)
      'MATLAB algorithm
       Dim k, i, no2 As Long
       Dim h(1 To NMAX) As Double
     
      ' Dim dft_xnre(1 To NMAX)     As Double
       Dim dft_xnim(1 To NMAX)     As Double
       Dim dft_Xkre(1 To NMAX)    As Double
       Dim dft_Xkim(1 To NMAX)    As Double
       Dim idft_Xkre(1 To NMAX)   As Double
       Dim idft_Xkim(1 To NMAX)   As Double
       Dim idft_xnre(1 To NMAX)   As Double
       Dim idft_xnim(1 To NMAX)   As Double

      If (Numb = 0) Then
         Hilbertre(1) = 0
         Hilbertim(1) = 0
       MsgBox (" Hilbert() :   Parsing error ")
        End If


      no2 = Int(Numb / 2)
     
      
      
      
       For i = 1 To Numb
     
          dft_xnim(i) = 0
       Next i

       Call DFT(Numb, xre, dft_xnim, dft_Xkre, dft_Xkim)

       i = 0
         Do While i <= Numb

            h(i + 1) = 0
             i = i + 1
        Loop

      'if  N is even
     
        If (Numb > 0) And (2 * no2) = Numb Then
            h(1) = 1
            h(no2 + 1) = 1
            For k = 2 To no2
               h(k) = 2
            Next k

       'if N is odd

        ElseIf (Numb > 0) Then

             h(1) = 1
              For k = 2 To (0.5 * (Numb + 1))
                h(k) = 2
              Next k
          
         End If


         For i = 1 To Numb
          idft_Xkre(i) = dft_Xkre(i) * h(i)
          idft_Xkim(i) = dft_Xkim(i) * h(i)
         Next i
     
        
         Call IDFT(Numb, idft_Xkre, idft_Xkim, idft_xnre, idft_xnim)
       
       For i = 1 To Numb
       Hilbertre(i) = idft_xnre(i)
       Hilbertim(i) = idft_xnim(i)
       Next i
     
      End Sub

     Sub GetAnalyticSignal(ByRef xnre() As Double, ByVal Numb As Long, ByRef x_re() As Double, ByRef x_im() As Double)

      'I=x(t)
       'Q=jHilb(t)
       'Hilb=Hilbert()
       Dim i As Long
       Dim Hilbertre(1 To NMAX) As Double
       Dim Hilbertim(1 To NMAX) As Double
       Call Hilbert(xnre, Numb, Hilbertre, Hilbertim)
       For i = 1 To Numb
        x_re(i) = xnre(i)
        x_im(i) = Hilbertim(i)
       Next i
      End Sub

     Sub getenvelope(ByRef xnre() As Double, ByVal Numb As Long, ByRef envelope() As Double)
       Dim i As Long
       Dim Hilbertre(1 To NMAX) As Double
       Dim Hilbertim(1 To NMAX) As Double
       Call Hilbert(xnre, Numb, Hilbertre, Hilbertim)
       For i = 1 To Numb
       envelope(i) = ((Hilbertre(i) ^ 2) + (Hilbertim(i) ^ 2)) ^ 0.5
       Next i
      End Sub


     Sub GetSpectrumIQ(ByRef xn() As Double, ByVal Num As Long, ByRef Xkre() As Double, ByRef Xkim() As Double)
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double

      Call GetAnalyticSignal(xn, Num, x_re, x_im)
      Call DFT(Num, x_re, x_im, Xkre, Xkim)
       Call ScaleDFT(Num, Xkre, Xkim)

     End Sub
     Sub GetR_fi(ByVal Numb As Long, ByRef Xkre() As Double, ByRef Xkim() As Double, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
       Dim i As Long
       For i = 1 To Numb
        modXk(i) = (Xkre(i) ^ 2 + Xkim(i) ^ 2) ^ 0.5
        argXk(i) = Arg(Xkre(i), Xkim(i))

      If mode = 1 Then
        argXk(i) = argXk(i) * 180 / M_PI
       End If

      Next i

     End Sub


    Sub GetSpectrum(ByRef xn() As Double, ByVal Num As Long, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double
      Dim Xkre(1 To NMAX) As Double
      Dim Xkim(1 To NMAX) As Double
          Call GetAnalyticSignal(xn, Num, x_re, x_im)
          Call DFT(Num, x_re, x_im, Xkre, Xkim)
          Call ScaleDFT(Num, Xkre, Xkim)
          Call GetR_fi(Num, Xkre, Xkim, modXk, argXk, mode)
     
      End Sub
     
     
      Sub GetSpectrumfromIQ(ByRef x_re() As Double, ByRef x_im() As Double, _
      ByVal Num As Long, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
      
      
      Dim Xkre(1 To NMAX) As Double
      Dim Xkim(1 To NMAX) As Double
          
          Call DFT(Num, x_re, x_im, Xkre, Xkim)
          Call ScaleDFT(Num, Xkre, Xkim)
          Call GetR_fi(Num, Xkre, Xkim, modXk, argXk, mode)
     
      End Sub


       Sub SynthSignalfromIQ(ByVal Num As Long, ByRef Xkre() As Double, ByRef Xkim() As Double, _
       ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
       
               
              'Call ScaleIQidft(Num, Xkre, Xkim)
              Call IDFT(Num, Xkre, Xkim, idft_xnre, idft_xnim)
              Call ScaleIQidft(Num, idft_xnre, idft_xnim)
       End Sub

       Sub SynthSignalfromRFi(ByRef modXk() As Double, ByRef argXk() As Double, _
        ByRef xnre() As Double, ByRef xnim() As Double, ByVal Num As Long, ByVal mode As Long)
        Dim Xk_re(1 To NMAX) As Double
        Dim Xk_im(1 To NMAX) As Double
         Call GetIQfromR_fi(Num, modXk, argXk, Xk_re, Xk_im, mode)
         Call SynthSignalfromIQ(Num, Xk_re, Xk_im, xnre, xnim)
       End Sub

    Sub printcomplex(ByRef xnre() As Double, ByRef xnim() As Double, _
     ByVal Num As Long, ByVal Col1 As Long, ByVal Col2 As Long)
         For i = 1 To Num
              Cells(i + 1, Col1).Value = xnre(i)
              Cells(i + 1, Col2).Value = xnim(i)
            Next i
     End Sub

    Sub printarray(ByRef xnre() As Double, ByVal Num As Long, ByVal Col1 As Long)
         For i = 1 To Num
            Cells(i + 1, Col1).Value = xnre(i)
            
            Next i
     End Sub

     Function deltat(ByRef argF() As Double, ByVal f0 As Double, ByVal fs As Double, ByVal N As Long) As Double
      Dim fv(1 To NMAX + 1) As Double
      Dim N0 As Long
      Dim i As Long
      Call Getfassoc(fs, fv, N)
     N0 = 0
     For i = 1 To N
      If (fv(i) >= f0) And (fv(i + 1) <= f0) Then
       N0 = i
       End If
       Next i
     
     deltat = -(argF(N0) - argF(0)) / (2 * M_PI * fv(N0) - 2 * M_PI * fv(0))
     End Function


     Sub writetofile(ByRef xnre() As Double, ByRef xnim() As Double, _
     ByRef FilePath As String, ByVal N As Long)
      Dim FilePath As String
      Dim i As Long

     Open FilePath For Output As #1
     
         'Write #1,message
         'Write #1,"\n ; ; ; ;   "
         ' If  mode=0  Then
             'Write #1,"\n n ; Re(x) ; +j* Im(x)      ;   "
             '  End If
         
          'If  mode=1 Then
             ' Write #1,"\n n ; |X|   ; exp(j* Arg(x) );   "
             '  End If

        For i = 1 To N
             Write #1, i
             Write #1, xnre(i)
             Write #1, xnim(i)
         Next i
            
         Write #1, "\n ; ; ; ;  "
         Close #1
         MsgBox (" data saved to file")
       
     End Sub

     Sub TestHilbert()
    Dim fs As Double
     Dim fv(1 To NMAX) As Double
     Dim tv(1 To NMAX) As Double
     Dim Xi(1 To NMAX) As Double
     Dim Xq(1 To NMAX) As Double

     Dim Yi(1 To NMAX) As Double
     Dim Yq(1 To NMAX) As Double
     Dim Hi(1 To NMAX) As Double
     Dim Hq(1 To NMAX) As Double
     Dim Ii(1 To NMAX) As Double
     Dim Iq(1 To NMAX) As Double
     Dim modX(1 To NMAX) As Double
     Dim argX(1 To NMAX) As Double
     Dim env(1 To NMAX) As Double
     Dim modX1(1 To NMAX) As Double
     Dim argX1(1 To NMAX) As Double
     Dim Xi1(1 To NMAX) As Double
     Dim i As Long
     Dim N As Long
    Dim x_re(1 To NMAX) As Double
    Dim Yre(1 To NMAX) As Double
    Dim Yim(1 To NMAX) As Double
    N = 8
     Xi(1) = 1
     Xi(2) = 2
     Xi(3) = 3
     Xi(4) = 4
     Xi(5) = 5
     Xi(6) = 6
     Xi(7) = 7
     Xi(8) = 8


     Xq(1) = 0
     Xq(2) = 0
     Xq(3) = 0
     Xq(4) = 0
     Xq(5) = 0
     Xq(6) = 0
     Xq(7) = 0
     Xq(8) = 0
     fs = 1000
      Call printcomplex(Xi, Xq, N, 1, 2)
      Call Getfassoc(fs, fv, N)
      Call Gettassoc(fs, tv, N)
      Call printcomplex(tv, fv, N, 3, 4)
     
     Call DFT(N, Xi, Xq, Yi, Yq)
     
     Call printcomplex(Yi, Yq, N, 5, 6)
     
       Call IDFT(N, Yi, Yq, Ii, Iq)
        Call printcomplex(Ii, Iq, N, 7, 8)
       Call Hilbert(Xi, N, Hi, Hq)
        Call printcomplex(Hi, Hq, N, 9, 10)
        Call GetSpectrum(Xi, N, modX, argX, 1)
         Call printcomplex(modX, argX, N, 11, 12)
      
     
       Call getenvelope(Xi, N, env)
       Call printarray(env, N, 13)
     
     
      
       N = 16
     
     
       For i = 0 To N - 1
        argX(i + 1) = 0
        modX(i + 1) = 0
        Next i
      
       modX(2) = 1
       argX(2) = -45
       modX(6) = 0
       argX(6) = 0
      
      '  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
      
       Call SynthSignalfromRFi(modX, argX, Xi, Xq, N, 1)
       Call printcomplex(Xi, Xq, N, 14, 15)
      
      
          Call Getfassoc(fs, fv, N)
          Call printarray(fv, N, 16)
         
         Call GetSpectrumfromIQ(Xi, Xq, N, modX1, argX1, 1)
         Call printcomplex(modX1, argX1, N, 17, 18)
         
         For i = 0 To N - 1
          Xi1(i + 1) = 0.5 * Cos(2 * M_PI * 250 * i / fs + 0) + 1 * Sin(2 * M_PI * 125 * i / fs + 0)
      
         Next i
          'Xi1(1) = 0
         
         Call GetSpectrum(Xi1, N, modX1, argX1, 1)
         Call printcomplex(modX1, argX1, N, 19, 20)
     
         Call GetSpectrumIQ(Xi1, N, Yre, Yim)
         Call printcomplex(Yre, Yim, N, 21, 22)
     
     End Sub

    Tuesday, November 17, 2015 9:37 PM
  • VBE Subroutines   for load/saveing data to file

    Option Explicit
     Const NMAX As Long = 65536
     
     
     '  ByRef xnre() As Double, ByRef xnim() As Double, _

     Sub ReadFromFileCplx(ByRef FilePath As String, ByRef N As Long, ByRef ind() As Long, _
      ByRef xnre() As Double, ByRef xnim() As Double, ByVal f As Integer, ByVal Comma As String)

      'Dim FilePath As String
       Dim i As Long

       Dim str As String
       Dim str1 As String
       Dim str2 As String
       Dim str3 As String
     ' Dim ind(1 To NMAX) As Long
       Dim word() As String

       str1 = ""
      Open FilePath For Input As #1
     
           Line Input #1, str1 ' Preamble          Input #1,str1
           Line Input #1, str1 ' "Function name"
          
           If f = 0 Then
           Line Input #1, str1 '   ; ; ; 8 ;
           word = Split(str1, ";")
            N = CLng(word(3))
            End If
           
             If f = 1 Then
           Input #1, str1, str1, str1, str1, str2  '   , , ,   8 ,
           
            N = CLng(str1)
            End If
           
          
           Line Input #1, str1 ' "  n ; Re(x) ; +j* Im(x)      ;   " or "\n n ; |X|   ; exp(j* Arg(x) );   "
     
      
        For i = 1 To N
          
         If f = 0 Then
        
        Line Input #1, str
        word = Split(str, ";")
        '  n1 = UBound(word)
        str1 = word(0)
        str2 = word(1)
        str3 = word(2)
        End If
       
        If f = 1 Then
        
         Input #1, str1, str2, str3
       
        End If
       
        ind(i) = CLng(Replace(str1, Comma, ","))
        If f = 0 Then
        ind(i) = ind(i) + 1
        End If
       
        xnre(i) = CDbl(Replace(str2, Comma, ","))
        xnim(i) = CDbl(Replace(str3, Comma, ","))
       
       '  Cells(i + 1, 1).Value = ind(i)
      '   Cells(i + 1, 2).Value = xnre(i)
       '  Cells(i + 1, 3).Value = xnim(i)
      
          Next i
            
        Line Input #1, str1   '  "\n ; ; ; ;  "
         Close #1
       '   MsgBox (" data loading  OK")
       
     End Sub
     
     
      Sub ReadFromFile(ByRef FilePath As String, ByRef N As Long, ByRef ind() As Long, _
      ByRef xnre() As Double, ByVal f As Integer, ByVal Comma As String)


       Dim i As Long
     
       Dim str As String
       Dim str1 As String
       Dim str2 As String
       Dim word() As String

       str1 = ""
      
      Open FilePath For Input As #1
     
           Line Input #1, str1 ' Preamble          Input #1,str1
           Line Input #1, str1 ' "Function name"
          
           If f = 0 Then   'default  format , which  uses  ';' as delim. symbol
           Line Input #1, str1 '   ; ; ; 8 ;
           word = Split(str1, ";")
            N = CLng(word(3))
            End If
           
             If f = 1 Then '  uses  ' "string1","string2", "string3"' as format , internal
           Input #1, str1, str1, str1, str1, str2  '   , , ,   8 ,
           
            N = CLng(str1)
            End If
           
          
           Line Input #1, str1 ' "  n ; Re(x) ; +j* Im(x)      ;   " or "\n n ; |X|   ; exp(j* Arg(x) );   "
     
      
        For i = 1 To N
          
         If f = 0 Then
        
        Line Input #1, str
        word = Split(str, ";")
        '  n1 = UBound(word)
        str1 = word(0)
        str2 = word(1)
       ' str3 = word(2)
        End If
       
        If f = 1 Then
        
         Input #1, str1, str2 ', str3
       
        End If
       
        ind(i) = CLng(Replace(str1, Comma, ","))
        If f = 0 Then
        ind(i) = ind(i) + 1
        End If
       
        xnre(i) = CDbl(Replace(str2, Comma, ","))
     '   xnim(i) = CDbl(Replace(str3, Comma, ","))
       
      '   Cells(i + 1, 1).Value = ind(i)
       '  Cells(i + 1, 2).Value = xnre(i)
      '   Cells(i + 1, 3).Value = xnim(i)
      
          Next i
            
        Line Input #1, str1   '  "\n ; ; ; ;  "
         Close #1
       '   MsgBox (" data loading  OK")
       
     End Sub
     
     
     
     Sub WriteToFileCplx(ByRef xnre() As Double, ByRef xnim() As Double, _
     ByRef FilePath As String, ByVal N As Long, ByVal mode As Integer)
     ' Dim FilePath As String
      Dim i As Long

     Open FilePath For Output As #1
     
          Write #1, , , , , " Data sequence     "
          Write #1, , , , , " Begin of sequence "
          Write #1, , , , N
           If mode = 0 Then
              Write #1, " n ", " Re(x) ", " +j* Im(x) "
                End If
         
           If mode = 1 Then
                Write #1, " n ", " |X|  ", " Exp(j * Arg(x)) "
                 
                End If

        For i = 1 To N
       
             
             Write #1, i, xnre(i), xnim(i)
       
             
         Next i
            
         Write #1, , , , , "  End of sequence "
         Close #1
         MsgBox ("Saved , OK")
       
     End Sub

     
     
     Sub WriteToFile(ByRef xn() As Double, _
     ByRef FilePath As String, ByVal N As Long)
     ' Dim FilePath As String
      Dim i As Long

     Open FilePath For Output As #1
     
          Write #1, , , , , " Data sequence     "
          Write #1, , , , , " Begin of sequence "
          Write #1, , , , N
        
              Write #1, , , , , " Data  "

        For i = 1 To N
             Write #1, i, xn(i)
         Next i
            
         Write #1, , , , , "  End of sequence "
         Close #1
         MsgBox ("  Saved , OK")
       
     End Sub
     
     
     
     
     
     
    Sub PrintCplx(ByRef xnre() As Double, ByRef xnim() As Double, ByVal num As Long, ByVal Col1 As Long, ByVal Col2 As Long)
      Dim i As Long
         For i = 1 To num
              Cells(i + 1, Col1).Value = xnre(i)
              Cells(i + 1, Col2).Value = xnim(i)
            Next i
     End Sub
     
     Sub GetCplxVal(ByRef xnre() As Double, ByRef xnim() As Double, ByVal num As Long, ByVal Col1 As Long, ByVal Col2 As Long)
      Dim i As Long
         For i = 1 To num
              xnre(i) = Cells(i + 1, Col1).Value
              xnim(i) = Cells(i + 1, Col2).Value
            Next i
     End Sub
     
     Sub PrintArray(ByRef xnre() As Double, ByVal num As Long, ByVal Col1 As Long)
          Dim i As Long
         For i = 1 To num
            Cells(i + 1, Col1).Value = xnre(i)
            
            Next i
     End Sub
     
     Sub GetArray(ByRef xnre() As Double, ByVal num As Long, ByVal Col1 As Long)
      Dim i As Long
          For i = 1 To num
            xnre(i) = Cells(i + 1, Col1).Value
            
            Next i
     End Sub

    Sub doInput()
           Dim ind(1 To NMAX) As Long
        Dim xre(1 To 100) As Double
        Dim xim(1 To 100) As Double
        Dim num As Long
       
       Call ReadFromFileCplx("file21.csv", num, ind, xre, xim, 1, ".")
       Call PrintCplx(xre, xim, num, 2, 3)
      
       Call WriteToFileCplx(xre, xim, "file21.csv", 8, 0)
     
     
     
     
     
    End Sub

     

    Thursday, November 19, 2015 7:27 PM
  • Example of format 0 of file (default, good for opening in Excel )

    \n
     ; xn1=idft(Xk/N)*N    ;  
     ; ; ; 8 ;  
     n ; Re(x) ; +j* Im(x)  ; 
        0 ;  1.000000e+000 ;  3.828427e+000 ;  
        1 ;  2.000000e+000 ;  -1.000000e+000 ;  
        2 ;  3.000000e+000 ;  -1.000000e+000 ;  
        3 ;  4.000000e+000 ;  -1.828427e+000 ;  
        4 ;  5.000000e+000 ;  -1.828427e+000 ;  
        5 ;  6.000000e+000 ;  -1.000000e+000 ;  
        6 ;  7.000000e+000 ;  -1.000000e+000 ;  
        7 ;  8.000000e+000 ;  3.828427e+000 ;  
     ; ; ; ; 
     

    Example of format 1 of file.txt (csv), internal for loading  data  to VBA programs, needs a macro

    ,,,," Data sequence     "
    ,,,," Begin of sequence "
    ,,,8
    " n "," Re(x) "," +j* Im(x) "
    1,1,3.828427
    2,2,-1
    3,3,-1
    4,4,-1.828427
    5,5,-1.828427
    6,6,-1
    7,7,-1
    8,8,3.828427
    ,,,,"  End of sequence "

    Thursday, November 19, 2015 7:30 PM
  • Matched filter  on  x(n)->GETIQ-> CFFT -> SCALE 1/num  - >Za->(mul) -> CIFFT - > scale * num-> y(n) - >

                                                                                                       |

                                             Za* (Re,Im) coeffs   ->-------------------*

    Option Explicit
     Const M_PI As Double = 3.14159265358979
     Const NMAX As Long = 131072

     Sub ReadFromFileCplx(ByRef FilePath As String, ByRef N As Long, ByRef ind() As Long, _
      ByRef xnre() As Double, ByRef xnim() As Double, ByVal f As Integer, ByVal Comma As String)

      'Dim FilePath As String
       Dim i As Long

       Dim str As String
       Dim str1 As String
       Dim str2 As String
       Dim str3 As String
     ' Dim ind(1 To NMAX) As Long
       Dim word() As String

       str1 = ""
      Open FilePath For Input As #1
     
           Line Input #1, str1 ' Preamble          Input #1,str1
           Line Input #1, str1 ' "Function name"
          
           If f = 0 Then
           Line Input #1, str1 '   ; ; ; 8 ;
           word = Split(str1, ";")
            N = CLng(word(3))
            End If
           
             If f = 1 Then
           Input #1, str1, str1, str1, str1, str2  '   , , ,   8 ,
           
            N = CLng(str1)
            End If
           
          
           Line Input #1, str1 ' "  n ; Re(x) ; +j* Im(x)      ;   " or "\n n ; |X|   ; exp(j* Arg(x) );   "
     
      
        For i = 1 To N
          
         If f = 0 Then
        
        Line Input #1, str
        word = Split(str, ";")
        '  n1 = UBound(word)
        str1 = word(0)
        str2 = word(1)
        str3 = word(2)
        End If
       
        If f = 1 Then
        
         Input #1, str1, str2, str3
       
        End If
       
        ind(i) = CLng(Replace(str1, Comma, ","))
        If f = 0 Then
        ind(i) = ind(i) + 1
        End If
       
        xnre(i) = CDbl(Replace(str2, Comma, ","))
        xnim(i) = CDbl(Replace(str3, Comma, ","))
       
       '  Cells(i + 1, 1).Value = ind(i)
      '   Cells(i + 1, 2).Value = xnre(i)
       '  Cells(i + 1, 3).Value = xnim(i)
      
          Next i
            
        Line Input #1, str1   '  "\n ; ; ; ;  "
         Close #1
       '   MsgBox (" data loading  OK")
       
     End Sub
     
     
      Sub ReadFromFile(ByRef FilePath As String, ByRef N As Long, ByRef ind() As Long, _
      ByRef xnre() As Double, ByVal f As Integer, ByVal Comma As String)


       Dim i As Long
     
       Dim str As String
       Dim str1 As String
       Dim str2 As String
       Dim word() As String

       str1 = ""
      
      Open FilePath For Input As #1
     
           Line Input #1, str1 ' Preamble          Input #1,str1
           Line Input #1, str1 ' "Function name"
          
           If f = 0 Then   'default  format , which  uses  ';' as delim. symbol
           Line Input #1, str1 '   ; ; ; 8 ;
           word = Split(str1, ";")
            N = CLng(word(3))
            End If
           
             If f = 1 Then '  uses  ' "string1","string2", "string3"' as format , internal
           Input #1, str1, str1, str1, str1, str2  '   , , ,   8 ,
           
            N = CLng(str1)
            End If
           
          
           Line Input #1, str1 ' "  n ; Re(x) ; +j* Im(x)      ;   " or "\n n ; |X|   ; exp(j* Arg(x) );   "
     
      
        For i = 1 To N
          
         If f = 0 Then
        
        Line Input #1, str
        word = Split(str, ";")
        '  n1 = UBound(word)
        str1 = word(0)
        str2 = word(1)
       ' str3 = word(2)
        End If
       
        If f = 1 Then
        
         Input #1, str1, str2 ', str3
       
        End If
       
        ind(i) = CLng(Replace(str1, Comma, ","))
        If f = 0 Then
        ind(i) = ind(i) + 1
        End If
       
        xnre(i) = CDbl(Replace(str2, Comma, ","))
     '   xnim(i) = CDbl(Replace(str3, Comma, ","))
       
      '   Cells(i + 1, 1).Value = ind(i)
       '  Cells(i + 1, 2).Value = xnre(i)
      '   Cells(i + 1, 3).Value = xnim(i)
      
          Next i
            
        Line Input #1, str1   '  "\n ; ; ; ;  "
         Close #1
       '   MsgBox (" data loading  OK")
       
     End Sub
     
     
     
     Sub WriteToFileCplx(ByRef xnre() As Double, ByRef xnim() As Double, _
     ByRef FilePath As String, ByVal N As Long, ByVal mode As Integer)
     ' Dim FilePath As String
      Dim i As Long

     Open FilePath For Output As #1
     
          Write #1, , , , , " Data sequence     "
          Write #1, , , , , " Begin of sequence "
          Write #1, , , , N
           If mode = 0 Then
              Write #1, " n ", " Re(x) ", " +j* Im(x) "
                End If
         
           If mode = 1 Then
                Write #1, " n ", " |X|  ", " Exp(j * Arg(x)) "
                 
                End If

        For i = 1 To N
       
             
             Write #1, i, xnre(i), xnim(i)
       
             
         Next i
            
         Write #1, , , , , "  End of sequence "
         Close #1
         MsgBox ("Saved , OK")
       
     End Sub

     
     
     Sub WriteToFile(ByRef xn() As Double, _
     ByRef FilePath As String, ByVal N As Long)
     ' Dim FilePath As String
      Dim i As Long

     Open FilePath For Output As #1
     
          Write #1, , , , , " Data sequence     "
          Write #1, , , , , " Begin of sequence "
          Write #1, , , , N
        
              Write #1, , , , , " Data  "

        For i = 1 To N
             Write #1, i, xn(i)
         Next i
            
         Write #1, , , , , "  End of sequence "
         Close #1
         MsgBox ("  Saved , OK")
       
     End Sub
     
     
     
     
     
     
    Sub PrintCplx(ByRef xnre() As Double, ByRef xnim() As Double, ByVal num As Long, ByVal Col1 As Long, ByVal Col2 As Long)
      Dim i As Long
         For i = 1 To num
              Cells(i + 1, Col1).Value = xnre(i)
              Cells(i + 1, Col2).Value = xnim(i)
            Next i
     End Sub
     
     Sub GetCplxVal(ByRef xnre() As Double, ByRef xnim() As Double, ByVal num As Long, ByVal Col1 As Long, ByVal Col2 As Long)
      Dim i As Long
         For i = 1 To num
              xnre(i) = Cells(i + 1, Col1).Value
              xnim(i) = Cells(i + 1, Col2).Value
            Next i
     End Sub
     
     Sub PrintArray(ByRef xnre() As Double, ByVal num As Long, ByVal Col1 As Long)
          Dim i As Long
         For i = 1 To num
            Cells(i + 1, Col1).Value = xnre(i)
            
            Next i
     End Sub
     
     Sub GetArray(ByRef xn() As Double, ByVal num As Long, ByVal Col1 As Long)
      Dim i As Long
          For i = 1 To num
            xn(i) = Cells(i + 1, Col1).Value
            
            Next i
     End Sub

    '*********************************************
     Sub Getfassoc(ByVal fs As Double, ByRef fvect() As Double, ByVal Numb As Long)
       Dim k As Long
        For k = 0 To Numb - 1
         fvect(k + 1) = fs * k / Numb
        Next k
      End Sub
     
     
      Sub Gettassoc(ByVal fs As Double, ByRef tvect() As Double, ByVal Numb As Long)
       Dim N As Long
        For N = 0 To Numb - 1
         tvect(N + 1) = N / fs
        Next N
      End Sub


      Function Arg(ByVal x_re As Double, ByVal x_im As Double) As Double
       'Dim x_re,x_im As Double
       'Dim arg As Double
       Dim message As String
      If (x_re > 0) Then
                 Arg = Atn(x_im / x_re)
                 End If
      If (x_re = 0) And (x_im > 0) Then
                 Arg = M_PI / 2
                 End If
      If (x_re = 0) And (x_im < 0) Then
                 Arg = -M_PI / 2
                 End If
      If (x_re < 0) And (x_im >= 0) Then
                 Arg = Atn(x_im / x_re) + M_PI
                 End If
      If (x_re < 0) And (x_im < 0) Then
                 Arg = Atn(x_im / x_re) - M_PI
                 End If
      If (x_re = 0) And (x_im = 0) Then
               '  MsgBox (" Xre=0 and Xim=0, Invalid result")
                 Arg = 0
                 End If

     End Function

     Sub GetIQfromR_fi(ByVal num As Long, ByRef modXk() As Double, ByRef argXk() As Double, _
     ByRef Xkre() As Double, ByRef Xkim() As Double, ByVal mode As Long)
       Dim i As Long
       Dim rad As Double
      
        rad = 1
         If mode = 1 Then
           rad = M_PI / 180
         End If
        
         For i = 1 To num
          Xkre(i) = modXk(i) * Cos(rad * argXk(i))
          Xkim(i) = modXk(i) * Sin(rad * argXk(i))
         Next i
        
       End Sub

     


     Sub CDFT(ByVal numdft As Long, ByRef dft_xnre() As Double, ByRef dft_xnim() As Double, ByRef dft_Xkre() As Double, ByRef dft_Xkim() As Double)

      Dim i, k, N As Long
     
     Dim angle As Double

     For k = 0 To numdft - 1
      dft_Xkre(k + 1) = 0
      dft_Xkim(k + 1) = 0
       For N = 0 To numdft - 1
        angle = -2 * M_PI * k * N / numdft
       dft_Xkre(k + 1) = dft_Xkre(k + 1) + dft_xnre(N + 1) * Cos(angle) - dft_xnim(N + 1) * Sin(angle)
       dft_Xkim(k + 1) = dft_Xkim(k + 1) + dft_xnre(N + 1) * Sin(angle) + dft_xnim(N + 1) * Cos(angle)
        Next N
      Next k

    End Sub
     
     Sub ScaleCDFT(ByVal num As Long, ByRef ScXkre() As Double, ByRef ScXkim() As Double)

    Dim i  As Long
      For i = 1 To num
         ScXkre(i) = ScXkre(i) / num
         ScXkim(i) = ScXkim(i) / num
      Next i
     End Sub
     Sub ScaleIQidft(ByVal num As Long, ByRef ScXk_re() As Double, ByRef ScXk_im() As Double)
       Dim i  As Long
       For i = 1 To num
        ScXk_re(i) = ScXk_re(i) * num
        ScXk_im(i) = ScXk_im(i) * num
       Next i
      End Sub

     Sub CIDFT(ByVal numidft As Long, ByRef idft_Xkre() As Double, _
     ByRef idft_Xkim() As Double, ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
     Dim k, N    As Long
      Dim angle As Double
     
      For N = 0 To numidft - 1
       'idft_item(n) = numidft - 1
       idft_xnre(N + 1) = 0
       idft_xnim(N + 1) = 0
        For k = 0 To numidft - 1
        angle = 2 * M_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)
        Next k
       
       idft_xnre(N + 1) = idft_xnre(N + 1) / numidft
       idft_xnim(N + 1) = idft_xnim(N + 1) / numidft

     Next N
     
      'idft_xnre(n), idft_xnim(n), numidft
     End Sub

     

     Sub Hilbert(ByRef xre() As Double, ByVal Numb As Long, ByRef Hilbertre() As Double, ByRef Hilbertim() As Double)
      'MATLAB algorithm
       Dim k, i, no2 As Long
       Dim h(1 To NMAX) As Double
     
      ' Dim dft_xnre(1 To NMAX)     As Double
       Dim dft_xnim(1 To NMAX)     As Double
       Dim dft_Xkre(1 To NMAX)    As Double
       Dim dft_Xkim(1 To NMAX)    As Double
       Dim idft_Xkre(1 To NMAX)   As Double
       Dim idft_Xkim(1 To NMAX)   As Double
       Dim idft_xnre(1 To NMAX)   As Double
       Dim idft_xnim(1 To NMAX)   As Double

      If (Numb = 0) Then
         Hilbertre(1) = 0
         Hilbertim(1) = 0
       MsgBox (" Hilbert() :   Parsing error ")
        End If


      no2 = Int(Numb / 2)
     
      
      
      
       For i = 1 To Numb
     
          dft_xnim(i) = 0
       Next i

       Call CDFT(Numb, xre, dft_xnim, dft_Xkre, dft_Xkim)

       i = 0
         Do While i <= Numb

            h(i + 1) = 0
             i = i + 1
        Loop

      'if  N is even
     
        If (Numb > 0) And (2 * no2) = Numb Then
            h(1) = 1
            h(no2 + 1) = 1
            For k = 2 To no2
               h(k) = 2
            Next k

       'if N is odd

        ElseIf (Numb > 0) Then

             h(1) = 1
              For k = 2 To (0.5 * (Numb + 1))
                h(k) = 2
              Next k
          
         End If


         For i = 1 To Numb
          idft_Xkre(i) = dft_Xkre(i) * h(i)
          idft_Xkim(i) = dft_Xkim(i) * h(i)
         Next i
     
        
         Call CIDFT(Numb, idft_Xkre, idft_Xkim, idft_xnre, idft_xnim)
       
       For i = 1 To Numb
       Hilbertre(i) = idft_xnre(i)
       Hilbertim(i) = idft_xnim(i)
       Next i
     
      End Sub

     Sub GetAnalyticSignal(ByRef xnre() As Double, ByVal Numb As Long, ByRef x_re() As Double, ByRef x_im() As Double)

      'I=x(t)
       'Q=jHilb(t)
       'Hilb=Hilbert()
       Dim i As Long
       Dim Hilbertre(1 To NMAX) As Double
       Dim Hilbertim(1 To NMAX) As Double
       Call Hilbert(xnre, Numb, Hilbertre, Hilbertim)
       For i = 1 To Numb
        x_re(i) = xnre(i)
        x_im(i) = Hilbertim(i)
       Next i
      End Sub

     Sub GetEnvelope(ByRef xnre() As Double, ByVal Numb As Long, ByRef envelope() As Double)
       Dim i As Long
       Dim Hilbertre(1 To NMAX) As Double
       Dim Hilbertim(1 To NMAX) As Double
       Call Hilbert(xnre, Numb, Hilbertre, Hilbertim)
       For i = 1 To Numb
       envelope(i) = ((Hilbertre(i) ^ 2) + (Hilbertim(i) ^ 2)) ^ 0.5
       Next i
      End Sub


     Sub GetSpectrumIQ(ByRef xn() As Double, ByVal num As Long, ByRef Xkre() As Double, ByRef Xkim() As Double)
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double

      Call GetAnalyticSignal(xn, num, x_re, x_im)
      Call CDFT(num, x_re, x_im, Xkre, Xkim)
       Call ScaleCDFT(num, Xkre, Xkim)

     End Sub
     Sub GetR_fi(ByVal Numb As Long, ByRef Xkre() As Double, ByRef Xkim() As Double, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
       Dim i As Long
       For i = 1 To Numb
        modXk(i) = (Xkre(i) ^ 2 + Xkim(i) ^ 2) ^ 0.5
        argXk(i) = Arg(Xkre(i), Xkim(i))

      If mode = 1 Then
        argXk(i) = argXk(i) * 180 / M_PI
       End If

      Next i

     End Sub


    Sub GetSpectrum(ByRef xn() As Double, ByVal num As Long, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double
      Dim Xkre(1 To NMAX) As Double
      Dim Xkim(1 To NMAX) As Double
          Call GetAnalyticSignal(xn, num, x_re, x_im)
          Call CDFT(num, x_re, x_im, Xkre, Xkim)
          Call ScaleCDFT(num, Xkre, Xkim)
          Call GetR_fi(num, Xkre, Xkim, modXk, argXk, mode)
     
      End Sub
     
     
      Sub GetSpectrumfromIQ(ByRef x_re() As Double, ByRef x_im() As Double, _
      ByVal num As Long, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Long)
      
      
      Dim Xkre(1 To NMAX) As Double
      Dim Xkim(1 To NMAX) As Double
          
          Call CDFT(num, x_re, x_im, Xkre, Xkim)
          Call ScaleCDFT(num, Xkre, Xkim)
          Call GetR_fi(num, Xkre, Xkim, modXk, argXk, mode)
     
      End Sub


       Sub SynthSignalfromIQ(ByVal num As Long, ByRef Xkre() As Double, ByRef Xkim() As Double, _
       ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
              'Call ScaleIQidft(Num, Xkre, Xkim)
              Call CIDFT(num, Xkre, Xkim, idft_xnre, idft_xnim)
              Call ScaleIQidft(num, idft_xnre, idft_xnim)
       End Sub

       Sub SynthSignalfromRFi(ByRef modXk() As Double, ByRef argXk() As Double, _
        ByRef xnre() As Double, ByRef xnim() As Double, ByVal num As Long, ByVal mode As Long)
        Dim Xk_re(1 To NMAX) As Double
        Dim Xk_im(1 To NMAX) As Double
         Call GetIQfromR_fi(num, modXk, argXk, Xk_re, Xk_im, mode)
         Call SynthSignalfromIQ(num, Xk_re, Xk_im, xnre, xnim)
       End Sub


    Sub GetMFCoesffs(ByRef xn() As Double, ByVal num As Long, ByRef Zasre() As Double, ByRef Zasim() As Double)
     
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double
       Dim Zare(1 To NMAX) As Double
       Dim Zaim(1 To NMAX) As Double
       Dim i As Long
      Call GetAnalyticSignal(xn, num, x_re, x_im)
      Call CDFT(num, x_re, x_im, Zare, Zaim)
      For i = 1 To num
      Zasre(i) = Zare(i) / num
      Zasim(i) = -Zaim(i) / num
      Next i

      End Sub


      Sub GetMatchedFilter(ByRef xn() As Double, ByVal num As Long, ByRef Zasre() As Double, ByRef Zasim() As Double, _
                        ByRef xoutre() As Double, ByRef xoutim() As Double)
     
      Dim x_re(1 To NMAX) As Double
      Dim x_im(1 To NMAX) As Double
       Dim Zare(1 To NMAX) As Double
       Dim Zaim(1 To NMAX) As Double
       'Dim xoutre (1 To NMAX)  As Double
       'Dim xoutim(1 To NMAX)   As Double
       Dim i As Long
      Call GetAnalyticSignal(xn, num, x_re, x_im)
      Call CDFT(num, x_re, x_im, Zare, Zaim)
      For i = 1 To num
      ' (a1+j*b1)*(a2+jb2)=(a1*a2-b1*b2) +j*(a1*b2+b1*a2)
      Zare(i) = (Zare(i) * Zasre(i) - Zaim(i) * Zasim(i)) / num
      Zaim(i) = (Zare(i) * Zasim(i) + Zaim(i) * Zasre(i)) / num
      Next i
      Call CIDFT(num, Zare, Zaim, xoutre, xoutim)
       For i = 1 To num
      xoutre(i) = xoutre(i) * num
      xoutim(i) = xoutim(i) * num
      Next i
     ' Parsing excl. of teta 0

      End Sub


    Sub GetRejectMatchedFilter(ByRef xn() As Double, ByVal num As Long, ByRef Zasre() As Double, ByRef Zasim() As Double, _
                        ByRef xoutre() As Double, ByRef xoutim() As Double)
     
     
       
       Dim xre(1 To NMAX) As Double
       Dim xim(1 To NMAX) As Double
       Dim yre(1 To NMAX) As Double
       Dim yim(1 To NMAX) As Double
       Dim i As Long
       Call GetAnalyticSignal(xn, num, xre, xim)
       Call GetMatchedFilter(xn, num, Zasre, Zasim, yre, yim)
      
       For i = 1 To num
       xoutre(i) = xre(i) - yre(i) ' or another  expression
       xoutim(i) = xim(i) - yim(i) ' or another  expression
      
       Next i
     
     
     
     
     
     
     
     End Sub
     
     '*************************************************

     

    Sub doInput()
           Dim ind(1 To NMAX) As Long
        Dim xre(1 To 100) As Double
        Dim xim(1 To 100) As Double
        Dim num As Long
       
       Call ReadFromFileCplx("file21.csv", num, ind, xre, xim, 1, ".")
       Call PrintCplx(xre, xim, num, 2, 3)
      
       Call WriteToFileCplx(xre, xim, "file21.csv", 8, 0)
     
     
     
     
     
    End Sub

      Sub GetSignalMF()
      Dim ModXkre(1 To NMAX) As Double
      Dim argXk(1 To NMAX)  As Double
      Dim xnre(1 To NMAX) As Double
      Dim xnim(1 To NMAX)  As Double
      Dim num As Long
       num = Cells(1, 1).Value
     
       Call GetArray(ModXkre, num, 2)
        Call GetArray(argXk, num, 3)
       
       Call SynthSignalfromRFi(ModXkre, argXk, xnre, xnim, num, 0)
      Call PrintCplx(xnre, xnim, num, 4, 5)
     
      End Sub
     
      Sub CoeffsMF()
       Dim xnre(1 To NMAX) As Double
      Dim xnim(1 To NMAX)  As Double
       Dim Zasre(1 To NMAX) As Double
       Dim Zasim(1 To NMAX) As Double
       Dim num As Long
       num = Cells(1, 1).Value
       Call GetArray(xnre, num, 4)
       'Call GetArray(xnim, num, 5)
       Call GetMFCoesffs(xnre, num, Zasre, Zasim)
       Call PrintCplx(Zasre, Zasim, num, 7, 8)
      End Sub
     
      Sub TestMF()
       Dim xnre(1 To NMAX) As Double
       Dim xnim(1 To NMAX)  As Double
       Dim Zasre(1 To NMAX) As Double
       Dim Zasim(1 To NMAX) As Double
       Dim xre(1 To NMAX) As Double
       Dim xim(1 To NMAX) As Double
       Dim env(1 To NMAX) As Double
      Dim num As Long
       num = Cells(1, 1).Value
       Call GetArray(xnre, num, 4)
       'Call GetArray(xnim, num, 5)
        Call GetArray(Zasre, num, 7)
        Call GetArray(Zasim, num, 8)
       Call GetMatchedFilter(xnre, num, Zasre, Zasim, xre, xim)
      
       Call PrintCplx(xre, xim, num, 10, 11)
       Call GetR_fi(num, xre, xim, xnre, xnim, 1)
      'Call GetEnvelope(xnre, num, env)
     
       Call PrintArray(xnre, num, 12)
      
      End Sub
     
     
      Sub TestRMF()
       Dim xnre(1 To NMAX) As Double
       Dim xnim(1 To NMAX)  As Double
       Dim Zasre(1 To NMAX) As Double
       Dim Zasim(1 To NMAX) As Double
       Dim xre(1 To NMAX) As Double
       Dim xim(1 To NMAX) As Double
       Dim env(1 To NMAX) As Double
      Dim num As Long
       num = Cells(1, 1).Value
       Call GetArray(xnre, num, 4)
       'Call GetArray(xnim, num, 5)
        Call GetArray(Zasre, num, 7)
        Call GetArray(Zasim, num, 8)
       Call GetRejectMatchedFilter(xnre, num, Zasre, Zasim, xre, xim)
      
       Call PrintCplx(xre, xim, num, 13, 14)
         Call GetR_fi(num, xre, xim, xnre, xnim, 1)
         Call PrintArray(xnre, num, 15)
      End Sub


    • Edited by USERPC01 Thursday, November 19, 2015 9:51 PM
    Thursday, November 19, 2015 7:35 PM
  • //getspectrum.cpp


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

    #define NMAX 1000


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

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


    CData  cdft(CData xn  ,int N )
     {
       
     int j;
     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 (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(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 main ()
    {

    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] : ");
     
    fprintf(fp2, "\n n;\t f, Hz;\t |Xf(jw)| ;\t  arg(Xf), deg ");
    if ((fp2=fopen("spectrum.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp2,"\n   %le   %d \n",fs,N1);

    for (i=0; i<N1;i++)
    {
    f=fs*i/N1;
    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 ");
    scanf("%s",&key);
    return 0;

    }

    Friday, December 2, 2016 11:45 PM
  • //synthsig.cpp
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>

    #define NMAX 1000


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

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


    CData  cdft(CData xn  ,int N )
     {
       
     int j;
     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 (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(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 ()
    {

    double xin[NMAX];
    CData  Xf;
    CData xt;
    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 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 ");
    scanf("%s",&key);
    return 0;

    }

    Friday, December 2, 2016 11:46 PM
  • //orthosig.cpp
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>

    #define NMAX 1000


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

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


    CData  cdft(CData xn  ,int N )
     {
       
     int j;
     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 (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(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 ()
    {

    double xin[NMAX];
    CData  Xf;
    CData xt;
    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 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 ");
    scanf("%s",&key);
    return 0;

    }

    Friday, December 2, 2016 11:47 PM
  • //getenvelope.cpp


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

    #define NMAX 1000


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

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


    CData  cdft(CData xn  ,int N )
     {
       
     int j;
     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 (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(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 main ()
    {

    double xin[NMAX];
    Array  Xenv;

    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 envelope of x[n]");
    Xenv=getenvelope(xin, N1);

    printf("\n Envelope of x[n] : ");
     
    fprintf(fp2, "\n n;\t f, Hz;\t |Xf(jw)| ;\t  arg(Xf), deg ");
    if ((fp2=fopen("envelope.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp2,"\n\n   %le   %d \n",fs,N1);

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

    }
    fclose(fp2);
    printf ("\n  input   '0' for exit ");
    scanf("%s",&key);
    return 0;

    }

    Friday, December 2, 2016 11:47 PM
  • //getenvratio.cpp
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>

    #define NMAX 1000


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

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


    CData  cdft(CData xn  ,int N )
     {
       
     int j;
     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 (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(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 main ()
    {
    double delta =1e-15;
    double xin[NMAX];
    double yout[NMAX];
    Array  Xenv1,Xenv2,Xenv3;
    double Kmax=0, Kmin=0, Kmid=0;
    char 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' ");
    scanf("%s",&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 ");
    scanf("%s",&key);
    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' ");
    scanf("%s",&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 ");
    scanf("%s",&key);
    printf("\n OK.   ");

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


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

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


    printf("\n |Hilbert(y[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++)
    {
      if (Xenv1.data[i]==0 ) { Xenv3.data[i]=Xenv2.data[i]/(Xenv1.data[i]+delta); }
      if (Xenv1.data[i]!=0 ) { Xenv3.data[i]=Xenv2.data[i]/(Xenv1.data[i] ); }


     
    printf ("\n %d  |Hilbert(y[n])|/|Hilbert(x[n])|= %le  ",  i,Xenv3.data[i]);
    fprintf (fp2,"\n %d   %le  ",  i,Xenv3.data[i]);

    }
    fclose(fp2);

    printf ("\n  input   '0' for continue ");
    scanf("%s",&key);

    Kmax=getmax(Xenv3,N1);
    Kmin=getmin(Xenv3,N1);
    Kmid=getmid(Xenv3,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 ");
    scanf("%s",&key);

    return 0;

    }

    Friday, December 2, 2016 11:48 PM
  • //spectrumratio.cpp
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>

    #define NMAX 1000


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

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


    CData  cdft(CData xn  ,int N )
     {
       
     int j;
     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 (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(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 main ()
    {
    double delta =1e-15;
    double xin[NMAX];
    double yout[NMAX];
    CData Xf1,Xf2,Xf3;
    double Kmax=0, Kmin=0, Kmid=0;
    char 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' ");
    scanf("%s",&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 ");
    scanf("%s",&key);
    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' ");
    scanf("%s",&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 ");
    scanf("%s",&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 ");
    scanf("%s",&key);

    return 0;

    }

    Friday, December 2, 2016 11:49 PM
  • //getsigratio.cpp
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>

    #define NMAX 1000


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

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


    CData  cdft(CData xn  ,int N )
     {
       
     int j;
     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 (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(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 main ()
    {
    double delta =1e-15;
    double xin[NMAX];
    double yout[NMAX];
    Array  Xr ;
    double Kmax=0, Kmin=0, Kmid=0;
    char 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' ");
    scanf("%s",&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 ");
    scanf("%s",&key);
    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' ");
    scanf("%s",&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 ");
    scanf("%s",&key);
    printf("\n OK.   ");

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


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

    label1:


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


    for (i=0; i<N1;i++)
    {
      if (xin[i]==0 ) { Xr.data[i]=yout[i]/(xin[i]+delta); }
     if (xin[i]!=0 ) { Xr.data[i]=yout[i]/(xin[i] ); }


     
    printf ("\n %d   y[n] / x[n] = %le  ",  i,Xr.data[i]);
    fprintf (fp2,"\n %d   %le  ",  i,Xr.data[i]);

    }
    fclose(fp2);

    printf ("\n  input   '0' for continue ");
    scanf("%s",&key);

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

    if ((fp2=fopen("yxcoeffs.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 ");
    scanf("%s",&key);

    return 0;

    }

    Friday, December 2, 2016 11:50 PM
  • //fir.cpp

    #include <stdio.h> // or  cstdio.h
    #include <stdlib.h> //or  cstdlib.h
    #include <math.h>   //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 < 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;
    }

     


    int main()
    {
     
    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;
    }

     

    Friday, December 2, 2016 11:51 PM