Asked by:
Obtainer of Hilbert transform (C++)
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];}
//fftgetfft(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+000Friday, 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",&l1);
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",&l2);
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 = 32767Sub 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 IfEnd Function
'36. - 4. + 9.6568542i - 4. + 4.i - 4. + 1.6568542i - 4. - 4. - 1.6568542i
' - 4. - 4.i - 4. - 9.6568542iSub 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 DoubleFor 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 kEnd 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 SubSub 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) / numidftNext n
'idft_xnre(n), idft_xnim(n), numidft
End SubSub 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 DoubleIf (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 iCall DFT(Numb, xre, dft_xnim, dft_Xkre, dft_Xkim)
i = 0
Do While i <= Numbh(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 SubSub 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 SubSub 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 DoubleCall 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 IfNext 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 SubSub 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 SubSub 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 SubSub 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 IntegerOpen 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 IfFor 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 SubSub 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 DoubleDim 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+308Sub 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 IfEnd 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 DoubleFor 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 kEnd 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 SubSub 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) / numidftNext N
'idft_xnre(n), idft_xnim(n), numidft
End SubSub 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 DoubleIf (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 iCall DFT(Numb, xre, dft_xnim, dft_Xkre, dft_Xkim)
i = 0
Do While i <= Numbh(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 SubSub 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 SubSub 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 DoubleCall 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 IfNext 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 SubSub 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 SubSub 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 SubSub 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 SubFunction 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 LongOpen 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 IfFor 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 SubSub 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 DoubleDim 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 SubTuesday, 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 LongDim 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 Stringstr1 = ""
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 Stringstr1 = ""
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 LongOpen 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 IfFor 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 LongOpen 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 SubSub 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 SubThursday, 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 = 131072Sub 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 LongDim 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 Stringstr1 = ""
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 Stringstr1 = ""
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 LongOpen 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 IfFor 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 LongOpen 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 IfEnd 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 DoubleFor 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 kEnd 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 SubSub 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) / numidftNext N
'idft_xnre(n), idft_xnim(n), numidft
End SubSub 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 DoubleIf (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 iCall CDFT(Numb, xre, dft_xnim, dft_Xkre, dft_Xkim)
i = 0
Do While i <= Numbh(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 SubSub 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 SubSub 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 DoubleCall 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 IfNext 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 SubSub 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 iEnd 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 0End 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 SubSub 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
