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

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-1IDFT
x[n]=(1/N)*sum(k=0;N-1;X[k]*exp(2*pi*k*n/N) ), n=0,...,N-1100*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