locked
Obtaining |H(jw)|, arg(H(jw)) from h(t) (h[n]) using DFT in C++ RRS feed

  • Question

  • Example of sampes, modules  for obtaining         |H(jw)|, arg(H(jw)) from h(t) (h[n])  using DFT in C++

    //ichar1.cpp

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

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


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

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

    */


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

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


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

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


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

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

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

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

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

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

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

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

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

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

    return 0;
    }

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

    }

    */

    Friday, November 25, 2016 11:32 AM

All replies

  •   Parser of one of realisation of h[n] from  H(jw), arg (H(jw)) (not unwrapped, not use GDT, may be with bugs )

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

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

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

    double getarg(float Re_H,float Im_H)
    {

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

    #define  NMAX 512
    int main()
    {
    int k,n;
    int w,t,N;
    double h[NMAX];
    double h_im[NMAX];
    double Re_H[NMAX];
    double Im_H[NMAX];
    double Re_x[NMAX];
    double Im_x[NMAX];
    double abs_H[NMAX];
    double arg_H[NMAX];
    double f[NMAX];
    float angle;
    double   K, fs,T,Tmax;
    char str;
    char key;
    printf("\n Input fs ,Hz ");
    scanf("%le",&fs);

    printf(" Input  number of points ");
    scanf("%d",&N);

    Tmax=N*1/fs;
    printf("\n Tmax =%lf s", Tmax);


    for (k=0 ; k<N;k++)
    {
    f[k]=k*fs/N; 
    abs_H[k]=1;
    arg_H[k]=0;


    printf("\n F=%le Hz; |H|[jw]=%le , arg (H)[jw]=%le ", f[k], abs_H[k], arg_H[k]  );
    }

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

    // H(jw)=integr(0;inf;h(tau)*exp(-j*w*tau) ; d tau)
    //K=2*3.1415926/fs;

    /*

    DFT
    X[k]=sum(n=0;N-1; x[n]*exp(-2*pi*j*k*n/N) ), k=0,...,N-1

    IDFT
    x[n]=(1/N)*sum(k=0;N-1;X[k]*exp(2*pi*k*n/N) ), n=0,...,N-1

    100*10^-6*500
    */


    label1:

     

    //cidft

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

     Re_x[n]=0;
     Im_x[n]=0;

      for(k=0;k<N;k++)
      {    
     
       angle=2*M_PI*k*n/N;   
            //Uoutp[w]/Uinp[w]   exp(j( fi out[w]-fi inp[w]) )
       Re_H[k]=abs_H[k]*cos(arg_H[k]);
       Im_H[k]=abs_H[k]*sin(arg_H[k]); 

        Re_x[n] = Re_x[n]+  Re_H[k]*cos(angle)-Im_H[k]*sin(angle);
        Im_x[n] = Im_x[n]+  Re_H[k]*sin(angle) + Im_H[k]*cos(angle);
      }
     Re_x[n]=Re_x[n]/N ;
     Im_x[n]=Im_x[n]/N ;
     
    // n scaled CIDFT

     //Re_x[n]=Re_x[n]*N;  
     //Im_x[n]=Im_x[n]*N;  
     
     
        printf("\n n=%d  Re(h(n))=%lf  Im(h(n))=%lf ",n, Re_x[n] ,Im_x[n]  ) ;

    }

    label3:


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

    return 0;
    }

    How to make it better to obtain h(t), k(t)=dh(t)/dt from |H(jw)|, arg((H(jw))), how to use smart phase unwrappper (is it  possible for some  anticipated state of   waveforms , for example  rectangular   ?)?

    Friday, November 25, 2016 11:55 AM
  • Simple parser of tau, fcut of RC filter

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

    /*
    X=-j*(1/(w*C))
     Z=sqrt(R^2+X^2)
     UR=U*R/Z;
     Ur=Uc=U/sqrt(2); AFR
     I = C(dU/dt),
     C(dU/dt)= -U/R
     ln|U| = - t/RC + Const
     U = e^-t/RC * e^Const
     U = e-t/RC * Const.

    UC = U(1 - e^-t/RC)


    UC = U(1 - exp(-t/RC)) integr
     t = RC, UC = U(1 - e^-1) = U(1 - 1/e) .
     UCð = Ue^(-t/tau) = U/e^(t/tau)

    K(jw)rcint=1/(jwRC+1);

    |K(jw)|=1/sqrt((wRC)^2+1)

    fi(w)=-arctg(wRC);

    tau=RC

    x(t)=1(t)

    h(t)int=1-exp(-t/tau);
     Uout(t)=Uinp(t)*k(t);
    k(t)=dh(t)/dt
    */


    double R,C,tau,fc;
    char key;
    int main ()
    {

    printf ("\n RC LPF  " );


    label1:

    printf ("\n Input R , Ohm " );
    scanf("%le",&R);
    printf ("\n Input C, uF " );
    scanf("%le", &C);
    C=C*10E-6;
    tau=R*C;
    printf (" \n tau=%le s", tau);
    fc=1/(2*M_PI*tau);

    printf (" \n fcut=%le Hz", fc);

    double f1,f2, df , f;
    int N,i;
    double argK[128];
    double absK[128];

    printf("\n Input f1 , Hz ");
    scanf("%le",&f1);
    printf("\n Input f2 , Hz ");
    scanf("%le",&f2);
    printf("\n Input delta f (step) , Hz ");
    scanf("%le",&df);
    printf("\n Input N  (< 128) ");
    scanf("%d",&N);

    f=f1;
    for(i=0;i<N;i++)
    {
    absK[i]=1/pow(((2*M_PI*f*R*C)*(2*M_PI*f*R*C)+1),0.5);
    argK[i]=-atan(2*M_PI*f*R*C)*180/M_PI;
    printf ("\nf= %le Hz; |K(jw)|=%le;fi(w)=%le deg",f,absK[i], argK[i]  );
    f=f+df;
    }


    printf ("\n Return - enter 0 \n New sequence- enter 1 :");
    scanf("%s",&key);
    if (key=='1') { goto label1 ; }
    return 0;
    }


    • Edited by USERPC01 Friday, November 25, 2016 12:21 PM patch ,fix bugs
    Friday, November 25, 2016 12:01 PM
  • Version  for buffered  RC LPF with n-th order (temporary in this program buffers with K(jw)=1*exp(j0) are ideal , problem of their impedances is not obtained  in this program ).

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

    /*
    X=-j*(1/(w*C))
     Z=sqrt(R^2+X^2)
     UR=U*R/Z;
     Ur=Uc=U/sqrt(2); AFR
     I = C(dU/dt),
     C(dU/dt)= -U/R
     ln|U| = - t/RC + Const
     U = e^-t/RC * e^Const
     U = e-t/RC * Const.

    UC = U(1 - e^-t/RC)


    UC = U(1 - exp(-t/RC)) integr
     t = RC, UC = U(1 - e^-1) = U(1 - 1/e) .
     UCр = Ue^(-t/tau) = U/e^(t/tau)

    K(jw)rcint=1/(jwRC+1);

    |K(jw)|=1/sqrt((wRC)^2+1)

    fi(w)=-arctg(wRC);

    tau=RC

    x(t)=1(t)

    h(t)int=1-exp(-t/tau);
     Uout(t)=Uinp(t)*k(t);
    k(t)=dh(t)/dt
    */


    double R,C,tau,fc;
    char key;
    int main ()
    {

    printf ("\n  'n -RC with ideal buffer'  LPF  with n-th order  " );


    label1:

    printf ("\n Input R , Ohm " );
    scanf("%le",&R);
    printf ("\n Input C, uF " );
    scanf("%le", &C);
    C=C*10E-6;
    tau=R*C;
    printf (" \n tau=%le s", tau);
    fc=1/(2*M_PI*tau);

    printf (" \n fcut=%le Hz", fc);

    double f1,f2, df , f;
    int N,i;
    double argK[128];
    double absK[128];

    printf("\n Input f1 , Hz ");
    scanf("%le",&f1);
    printf("\n Input f2 , Hz ");
    scanf("%le",&f2);
    printf("\n Input delta f (step) , Hz ");
    scanf("%le",&df);
    printf("\n Input N  (< 128) ");
    scanf("%d",&N);
    printf("\n Input number of cascades  (1...100) ");
    int N1;
    scanf("%d",&N1);

    f=f1;
    for(i=0;i<N;i++)
    {
    absK[i]=1/pow(((2*M_PI*f*R*C)*(2*M_PI*f*R*C)+1),0.5);
    absK[i]=pow(absK[i],N1);
    argK[i]=-atan(2*M_PI*f*R*C)*180/M_PI;
    argK[i]=N1*argK[i];
    printf ("\nf= %le Hz; |K(jw)|=%le;fi(w)=%le deg ,order=%d",f,absK[i], argK[i], N1  );
    f=f+df;
    }


    printf ("\n Return - enter 0 \n New sequence- enter 1 :");
    scanf("%s",&key);
    if (key=='1') { goto label1 ; }
    return 0;
    }

    Friday, November 25, 2016 12:34 PM
  • // RCF2.cpp

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

    double f1,f2, df , f;
    int N,i;
    double argK[128];
    double absK[128];
    int N1;
    /*
    X=-j*(1/(w*C))
     Z=sqrt(R^2+X^2)
     UR=U*R/Z;
     Ur=Uc=U/sqrt(2); AFR
     I = C(dU/dt),
     C(dU/dt)= -U/R
     ln|U| = - t/RC + Const
     U = e^-t/RC * e^Const
     U = e-t/RC * Const.

    UC = U(1 - e^-t/RC)


    UC = U(1 - exp(-t/RC)) integr
     t = RC, UC = U(1 - e^-1) = U(1 - 1/e) .
     UCð = Ue^(-t/tau) = U/e^(t/tau)

    K(jw)rcint=1/(jwRC+1);

    |K(jw)|=1/sqrt((wRC)^2+1)

    fi(w)=-arctg(wRC);

    tau=RC

    x(t)=1(t)

    h(t)int=1-exp(-t/tau);
     Uout(t)=Uinp(t)*k(t);
    k(t)=dh(t)/dt
    */


    double R,C,tau,fc;
    char key;

    int main ()
    {

    printf ("\n  'n -RC  circuits based   LPF, HPF  " );


    label1:

    printf ("\n Input R , Ohm " );
    scanf("%le",&R);
    printf ("\n Input C, uF " );
    scanf("%le", &C);
    C=C*10E-6;
    tau=R*C;
    printf (" \n tau=%le s", tau);
    fc=1/(2*M_PI*tau);

    printf (" \n fcut=%le Hz", fc);


    printf ("\n for  RC  LPF  input 1 ");
    printf ("\n for  RC  HPF  input 2 ");
    printf ("\n for  exit     input 0 : ");
    scanf("%s",&key);

    if (key=='1') { goto label_int ; }
    if (key=='2') { goto label_dif ; }
    if (key=='0') { goto label_int ; }

    label_int:


    printf("\n Input f1 , Hz ");
    scanf("%le",&f1);
    printf("\n Input f2 , Hz ");
    scanf("%le",&f2);
    printf("\n Input delta f (step) , Hz ");
    scanf("%le",&df);
    printf("\n Input N  (< 128) ");
    scanf("%d",&N);
    printf("\n Input number of cascades  (1...100) ");

    scanf("%d",&N1);


    f=f1;
    for(i=0;i<N;i++)
    {
    absK[i]=1/pow(((2*M_PI*f*R*C)*(2*M_PI*f*R*C)+1),0.5);
    absK[i]=pow(absK[i],N1);
    argK[i]=-atan(2*M_PI*f*R*C)*180/M_PI;
    argK[i]=N1*argK[i];
    printf ("\nf= %le Hz; |K(jw)|lpf=%le;fi(w)=%le deg ,order=%d",f,absK[i], argK[i], N1  );
    f=f+df;
    }
     goto end_label;


    label_dif:

    //
    //  K(jw)=jwRC/(1+jwRC)
    //  |K|=wRC/sqrt(1+(wRC)^2)
    //  fi=arctg(1/wRC)
    //  fc=1/(2*pi*RC)
    //  f<< fc ->Uout=RC*(dUinp/dt)
    //  fc=sqrt(sum (i; fcut ^2))
    //

     

    printf("\n Input f1 , Hz ");
    scanf("%le",&f1);
    printf("\n Input f2 , Hz ");
    scanf("%le",&f2);
    printf("\n Input delta f (step) , Hz ");
    scanf("%le",&df);
    printf("\n Input N  (< 128) ");
    scanf("%d",&N);
    printf("\n Input number of cascades  (1...100) ");
     
    scanf("%d",&N1);

    printf (" \n fcut of one circuit=%le Hz", fc);


    f=0;
    for (i=0;i<N;i++) {f+=fc*fc; }
    f=pow(f,0.5);
    printf("\n fcut=%le hz \n",f);

    f=f1;
    for(i=0;i<N;i++)
    {
    absK[i]=2*M_PI*f*R*C/pow(((2*M_PI*f*R*C)*(2*M_PI*f*R*C)+1),0.5);
    absK[i]=pow(absK[i],N1);
    argK[i]= atan(1/(2*M_PI*f*R*C))*180/M_PI;
    argK[i]=N1*argK[i];
    printf ("\nf= %le Hz; |K(jw)hpf|=%le;fi(w)=%le deg ,order=%d",f,absK[i], argK[i], N1  );
    f=f+df;
    }


    end_label:


    printf ("\n Return - enter 0 \n New sequence- enter 1 :");
    scanf("%s",&key);
    if (key=='1') { goto label1 ; }
    return 0;
    }

    Tuesday, November 29, 2016 11:43 AM
  • //ichar2.cpp

    #include <stdio.h> // or  cstdio.h
    #include <stdlib.h> //or  cstdlib.h
    #include <math.h>   //or  cmath.h

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


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

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

    */


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

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


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

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


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

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

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

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

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

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

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

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

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

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

    return 0;
    }

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

    }

    */

    Wednesday, November 30, 2016 11:50 AM
  • //afrtohtcoefs.cpp
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>

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


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

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

    double getarg(float Re_H,float Im_H)
    {

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


    int main()
    {


    printf ("\n For input  |H(jw)|, fi(w) from file input '1' ");
    printf ("\n For manual  input |H(jw)|, fi(w) input   '0' ");
    scanf("%s",&key);
    if (key=='0' ) { goto label_02;}
    if (key=='1' ) { goto label_01;}

    label_01:
     printf("\n Open afrcoefs.txt for reading AFR coefficients ");
    if ((fp=fopen("afrcoefs.txt","r"))==NULL)  printf ("\n File opening error");
    fscanf(fp,"%le",&fs);
    fscanf(fp,"%d",&N);
    Tmax=N*1/fs;
    for (k=0 ; k<N;k++)
    {
    // f[k]=k*fs/N;
    //fprintf(fp1,"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
     fscanf (fp,"%d",&n);
     fscanf (fp,"%le",&abs_H[k]);
     fscanf (fp,"%le",&arg_H[k]);
     
      fscanf (fp,"%le",&f[k]);
     
     arg_H[k]=arg_H[k]*M_PI/180;
    }
    fclose(fp);
    goto label_03;

    label_02:

    printf("\n Input fs ,Hz ");
    scanf("%le",&fs);
    printf(" Input  number of points ");
    scanf("%d",&N);
    Tmax=N*1/fs;
    printf("\n Tmax =%lf s", Tmax);
    for (k=0 ; k<N;k++)
    {
    f[k]=k*fs/N;

    //fprintf(fp1,"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
          printf("\n F=%le Hz" ,f[k]);
       printf("\n k=%d;    |H(jw)|[k]=",k );
       scanf ( "%le",&abs_H[k]);
          printf(" k=%d;   arg(H(jw))[k]=",k );
          scanf ( "%le",&arg_H[k]);
    //abs_H[k]=1;
    //arg_H[k]=0;
    arg_H[k]=arg_H[k]*M_PI/180;
    }

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


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

    label1:


    //cidft
    printf("\n Open hcoefs1.txt for writing coeffs h[n] ");
    if ((fp1=fopen("hcoefs1.txt","a"))==NULL)  printf ("\n File opening error");
    fprintf(fp1,"\n\n  %le   %d", fs, N );
    for (n=0 ; n<N; n++)
    {
     Re_x[n]=0;
     Im_x[n]=0;
      for(k=0;k<N;k++)
      {    
       angle=2*M_PI*k*n/N;   
            //Uoutp[w]/Uinp[w]   exp(j( fi out[w]-fi inp[w]) )
       Re_H[k]=abs_H[k]*cos(arg_H[k]);
       Im_H[k]=abs_H[k]*sin(arg_H[k]); 
        Re_x[n] = Re_x[n]+  Re_H[k]*cos(angle)-Im_H[k]*sin(angle);
        Im_x[n] = Im_x[n]+  Re_H[k]*sin(angle) + Im_H[k]*cos(angle);
      }
     Re_x[n]=Re_x[n]/N ;
     Im_x[n]=Im_x[n]/N ;
     // n scaled CIDFT  , not used in this, use if sygnal synth.
     //Re_x[n]=Re_x[n]*N;  
     //Im_x[n]=Im_x[n]*N;  
        //printf("\n n=%d  Re(h(n))=%lf  Im(h(n))=%lf ",n, Re_x[n] ,Im_x[n]  ) ;
        printf("\n n=%d   h[n]=%lf   ",n, Re_x[n]    ) ;
        fprintf(fp1, "\n  %d   %le   ",n, Re_x[n]  ) ; 
    }

     fclose (fp1);
     printf( "\n  Re(h(n)=h[n]), use it to load into FIR  "   ) ;
    label3:
    printf("\n Input 0 for exit ");
    scanf("%s",&key);
    return 0;
    }

    Wednesday, November 30, 2016 11:51 AM
  • //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;
    }

     

    Wednesday, November 30, 2016 7:55 PM