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

All replies

• 

  

  0

  Sign in to vote

    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
• 

  

  0

  Sign in to vote

  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
• 

  

  0

  Sign in to vote

  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
• 

  

  0

  Sign in to vote

  // 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
• 

  

  0

  Sign in to vote

  //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
• 

  

  0

  Sign in to vote

  //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
• 

  

  0

  Sign in to vote

  //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