# 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];

}

//fft

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

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

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

h[i]=0;
}

if((n>0)&&((2*floor(n/2))==(n))) // n is even
{
printf("\n n is even ");
h[0]=1;
// h[n/2]=1;
h[n/2]=1;

for(k=1;k<n/2;k++)
{
h[k]=2;
}
}

else //n is odd
{
printf("\n n is odd ");
if(n>0)
{
h[0]=1;

for(k=1;k<(n+1)/2;k++)
{
h[k]=2;
}
}

}

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

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

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

return ;
}

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

for (j=0;j<num;j++)
{

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

hilbert(m, num);

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

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

hilbert(m, num);

return 0;
}

Friday, October 23, 2015 3:58 AM

### All replies

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

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

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

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

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

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

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

//dftlib.h

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

#define NMAX 1000

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

typedef struct  {
double data[NMAX];

} DArray;

CVector  dft(CVector xn  ,int N )
{

int j;
CVector Xk ;
int k,n  ;
double angle;

for (k=0;k<N;k++)
{
Xk.re[k]=0;
Xk.im[k]=0;
for(n=0;n<N;n++)
{
angle=-2*M_PI*k*n/N;
//dft_Xk->re[k]=dft_Xk->re[k]+dft_xn->re[n]*cos(angle);
//dft_Xk->im[k]=dft_Xk->im[k]+dft_xn->re[n]*sin(angle);

//  if angle with '-'
//a=xnre
//b=xnim
//c=cos
//d=sin
//(a+bi)(c+di)=  (ac -bd) + (ad + bc)i
Xk.re[k]=Xk.re[k]+xn.re[n]*cos(angle)- xn.im[n]*sin(angle);
Xk.im[k]=Xk.im[k]+xn.re[n]*sin(angle)+ xn.im[n]*cos(angle);
}

}

return Xk;
}

CVector DoScaleDFT(CVector xk,int  num)
{
int i;
for (i=0; i<num ; i++)
{
xk.re[i]=xk.re[i]/num;
xk.im[i]=xk.im[i]/num;
}
return xk;
}

CVector  idft(CVector idft_Xk, int N  )
{

int k,n ;
double angle;
CVector  idft_xn;

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

idft_xn.re[n]=idft_xn.re[n]+idft_Xk.re[k]*cos(angle)-idft_Xk.im[k]*sin(angle);
idft_xn.im[n]=idft_xn.im[n]+idft_Xk.re[k]*sin(angle)+idft_Xk.im[k]*cos(angle);
}
idft_xn.im[n]=idft_xn.im[n]/N;
idft_xn.re[n]=idft_xn.re[n]/N;

}

return idft_xn ;
}

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

return X;
}

CVector  get_r_fi(CVector x, int n, int mode   )
{
int i;
CVector res;
double fi;
for (i=0;i<n;i++)
{
res.re[i]=pow((x.re[i]*x.re[i]+x.im[i]*x.im[i]),0.5);
//double treshh=1e9;
// x.im[i]=x.im[i]*treshh;
// x.im[i]=floor(x.im[i])/treshh;
// x.re[i]=x.re[i]*treshh;
// x.re[i]=floor(x.re[i])/treshh;
res.im[i]=Arg(x.re[i],x.im[i]);
if (mode==1)
{ res.im[i]=res.im[i]*180/M_PI;  }

}
return  res;
}

CVector ScaleIQidft(CVector Xk, int num)
{
int i;
for (i=0;  i<num ; i++)
{
Xk.re[i]=Xk.re[i]*num;
Xk.im[i]=Xk.im[i]*num;
}

return Xk;

}

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

return res;
}

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

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

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

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

void printarray(DArray xn , int n)
{
int i;
printf("\n   " );
for (i=0;i<n;i++)
{
printf("\n %d   %le  ",i,xn.data[i] );
}
printf("\n " );
return;
}
void writetofile(CVector xn , int n,  char *fname, char *message,int  mode      )
{
FILE *fp;
int i;

fp=fopen(fname,"a");

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

void writetofile_darray(DArray xn , int n,  char *fname, char *message     )
{
FILE *fp;
int i;

fp=fopen(fname,"a");
fprintf(fp,"\n ; %s ;  ; " , message );
fprintf(fp,"\n ; ;  " );
fprintf(fp,"\n n ; Array ;   " );

for (i=0;i<n;i++)
{
fprintf(fp,"\n    %d ;  %le ;    ",i,xn.data[i]  );
}
fprintf(fp,"\n ; ;    " );
fclose(fp);
printf("\n data saved to file   %s    ", fname);
return;
}

CVector hilbert(CVector xr , int n)
{
int  k,i,no2;
double h[NMAX];
CVector dft__xn,dft__Xk,idft__XK,idft__XN;

if (n==0)
{
printf("\n Error");
idft__XN.re[0]=0;
idft__XN.im[0]=0;
return idft__XN;
}
//xnreal=real(xr);
no2=floor (n/2);
for (i=0;i<n;i++)
{
dft__xn.re[i]=xr.re[i];
dft__xn.im[i]=0;
}

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

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

if((n>0)&&((2*floor(n/2))==(n))) // n is even
{

h[0]=1;
h[n/2]=1;
for(k=1;k<n/2;k++)
{
h[k]=2;
}
}
else //n is odd
{
if(n>0)
{
h[0]=1;
for(k=1;k<(n+1)/2;k++)
{
h[k]=2;
}
}
}

for(i=0;i<n;i++)
{
idft__XK.re[i]= dft__Xk.re[i]*h[i];
idft__XK.im[i]= dft__Xk.im[i]*h[i];
}

idft__XN= idft(idft__XK, n);

return idft__XN;
}

CVector getanalyticsignal(CVector xn, int num)
{
CVector x_n,Hilb;
int i;
//x=x+jy,y=0
//hilb(x)=a+jb
//z=x+jHilb(x)=x+j(a+jb)=x-b+j*a
//I=x(t), Q=j*Hilb(xre(t))

Hilb=hilbert(xn,num);

for (i=0; i<num; i++)
{
x_n.re[i]=xn.re[i];
x_n.im[i]=Hilb.im[i];
}

return x_n;
}

DArray getenvelope(CVector x, int num)
{
int i;
CVector  y;
DArray env;
y=hilbert(x,num);
for (i=0;i<num;i++)
{

env.data[i]=pow(((y.re[i]*y.re[i])+(y.im[i]*y.im[i])),0.5);

}

return env;
}

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

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

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

return Xk;
}

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

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

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

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

//hilb (sin(x)) =-cos(x)
//hilb(cos(x)) = sin(x)
//hilb (exp(it)=-i*exp(it)
//hilb(exp(-it))=i*exp(-it)
//x = 2.5cos(2*pi*20.3*t)+sin(2*pi*72.1*t)+cos(2*pi*100.1*t);
//y = hilbert(x);
//printf("\n xn=  ");
//xn=pokecomplex(m,n,num);
//double inph[NMAX];

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

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

double deg=M_PI/180;

CVector res;

int i;
*fs=1000 ;
*num=*fs  ;
int Nmax=*num;
//inph[i]=m[i];

double ampl1=1;
double freq1=10;
double phase1=90;
double ampl2=0.5;
double freq2=72;
double phase2=35;

printf("\n Input A1 (0.05) ");
scanf("%le",&ampl1);
printf("\n Input f1 (10) ");
scanf("%le",&freq1);
printf("\n Input phase1 ,deg (90) ");
scanf("%le",&phase1);
printf("\n Input A2 (0.1)");
scanf("%le",&ampl2);
printf("\n Input f1 (72) ");
scanf("%le",&freq2);
printf("\n Input phase1 ,deg (35) ");
scanf("%le",&phase2);

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

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

res.im[i]=0;
}

getch();

return res;
}

• Edited by 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];

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

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

}

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

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);
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] ;

}

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 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;
for (i=0;i<numb;i++)
{
//10   ; 1.000000e+000 ;  9.000000e+001 ;
}

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 Thursday, October 29, 2015 10:28 PM
Thursday, October 29, 2015 4:28 PM
• //testspectrum.cpp

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

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

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

getch();

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

getch();

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

return 0;

}

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

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

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

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

getch();

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

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

Xk=get_r_fi(Xk, num, 1  );

xn1=SynthSignalfromRFi(Xk, num,1);

writetofile( xn1 , num, "signalidft.csv", "xre(t) signal ",0 );
return 0;

}

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

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

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

int i,j;
double *xkr,*xki;
CVector xn,xn1,hilb,hexp,fass;
DArray  X;
double fs;
getch();
printf ("\n init  xn   ");
xn=do_initsignal(&fs, &num);
printf("\n xn=  ");
// printcomplex(xn,num);
writetofile( xn , num, "signal.csv", "xre(t) signal ",0 );

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

//analytic signal

xn1=getanalyticsignal(  xn,   num);

writetofile( xn1 , num, "analytic.csv", "analytic signal I(t),Q(t)  ",1 );
getch();

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

return 0;
}

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

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

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

double deg=M_PI/180;

CVector res;

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

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

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

ampl1=1;
ampl2=1;
freq1=5;
freq2=100;

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

res.im[i]=0;
}

getch();

return res;
}

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

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

VBA   script for Excel (ALT+F11)

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

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

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

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

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

Sub GetIQfromR_fi(ByVal Num As Integer, ByRef modXk() As Double, ByRef argXk() As Double, _
ByRef Xkre() As Double, ByRef Xkim() As Double, ByVal mode As Integer)
Dim i As Integer

If mode = 1 Then
End If

For i = 1 To Num
Xkre(i) = modXk(i) * Cos(rad * argXk(i))
Xkim(i) = modXk(i) * Sin(rad * argXk(i))
Next i

End Sub

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

Dim i, k, n As Integer

Dim angle As Double

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

End Sub

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

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

Sub IDFT(ByVal numidft As Integer, ByRef idft_Xkre() As Double, _
ByRef idft_Xkim() As Double, ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
Dim k, n    As Integer
For n = 0 To numidft - 1
'idft_item(n) = numidft - 1
idft_xnre(n + 1) = 0
idft_xnim(n + 1) = 0
For k = 0 To numidft - 1
angle = 2 * M_PI * k * n / numidft
idft_xnre(n + 1) = idft_xnre(n + 1) + idft_Xkre(k + 1) * Cos(angle) - idft_Xkim(k + 1) * Sin(angle)
idft_xnim(n + 1) = idft_xnim(n + 1) + idft_Xkre(k + 1) * Sin(angle) + idft_Xkim(k + 1) * Cos(angle)
Next k

idft_xnre(n + 1) = idft_xnre(n + 1) / numidft
idft_xnim(n + 1) = idft_xnim(n + 1) / numidft

Next n

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

Sub Hilbert(ByRef xre() As Double, ByVal Numb As Integer, ByRef Hilbertre() As Double, ByRef Hilbertim() As Double)
'MATLAB algorithm
Dim k, i, no2 As Integer
Dim h(1 To NMAX) As Double

' Dim dft_xnre(1 To NMAX)     As Double
Dim dft_xnim(1 To NMAX)     As Double
Dim dft_Xkre(1 To NMAX)    As Double
Dim dft_Xkim(1 To NMAX)    As Double
Dim idft_Xkre(1 To NMAX)   As Double
Dim idft_Xkim(1 To NMAX)   As Double
Dim idft_xnre(1 To NMAX)   As Double
Dim idft_xnim(1 To NMAX)   As Double

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

no2 = Int(Numb / 2)

For i = 1 To Numb

dft_xnim(i) = 0
Next i

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

i = 0
Do While i <= Numb

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

'if  N is even

If (Numb > 0) And (2 * no2) = Numb Then
h(1) = 1
h(no2 + 1) = 1
For k = 2 To no2
h(k) = 2
Next k

'if N is odd

ElseIf (Numb > 0) Then

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

End If

For i = 1 To Numb
idft_Xkre(i) = dft_Xkre(i) * h(i)
idft_Xkim(i) = dft_Xkim(i) * h(i)
Next i

Call IDFT(Numb, idft_Xkre, idft_Xkim, idft_xnre, idft_xnim)

For i = 1 To Numb
Hilbertre(i) = idft_xnre(i)
Hilbertim(i) = idft_xnim(i)
Next i

End Sub

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

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

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

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

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

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

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

Next i

End Sub

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

End Sub

Sub GetSpectrumfromIQ(ByRef x_re() As Double, ByRef x_im() As Double, _
ByVal Num As Integer, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Integer)

Dim Xkre(1 To NMAX) As Double
Dim Xkim(1 To NMAX) As Double

Call DFT(Num, x_re, x_im, Xkre, Xkim)
Call ScaleDFT(Num, Xkre, Xkim)
Call GetR_fi(Num, Xkre, Xkim, modXk, argXk, mode)

End Sub

Sub SynthSignalfromIQ(ByVal Num As Integer, ByRef Xkre() As Double, ByRef Xkim() As Double, _
ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)

'Call ScaleIQidft(Num, Xkre, Xkim)
Call IDFT(Num, Xkre, Xkim, idft_xnre, idft_xnim)
Call ScaleIQidft(Num, idft_xnre, idft_xnim)
End Sub

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

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

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

Next i
End Sub

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

Open FilePath For Output As #1

'Write #1,message
'Write #1,"\n ; ; ; ;   "
' If  mode=0  Then
'Write #1,"\n n ; Re(x) ; +j* Im(x)      ;   "
'  End If

'If  mode=1 Then
' Write #1,"\n n ; |X|   ; exp(j* Arg(x) );   "
'  End If

For i = 1 To n
Write #1, i
Write #1, xnre(i)
Write #1, xnim(i)
Next i

Write #1, "\n ; ; ; ;  "
Close #1
MsgBox (" data saved to file")

End Sub

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

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

Xq(1) = 0
Xq(2) = 0
Xq(3) = 0
Xq(4) = 0
Xq(5) = 0
Xq(6) = 0
Xq(7) = 0
Xq(8) = 0
fs = 1000
Call printcomplex(Xi, Xq, n, 1, 2)
Call Getfassoc(fs, fv, n)
Call Gettassoc(fs, tv, n)
Call printcomplex(tv, fv, n, 3, 4)

Call DFT(n, Xi, Xq, Yi, Yq)

Call printcomplex(Yi, Yq, n, 5, 6)

Call IDFT(n, Yi, Yq, Ii, Iq)
Call printcomplex(Ii, Iq, n, 7, 8)
Call Hilbert(Xi, n, Hi, Hq)
Call printcomplex(Hi, Hq, n, 9, 10)
Call GetSpectrum(Xi, n, modX, argX, 1)
Call printcomplex(modX, argX, n, 11, 12)

Call getenvelope(Xi, n, env)
Call printarray(env, n, 13)

n = 16

For i = 0 To n - 1
argX(i + 1) = 0
modX(i + 1) = 0
Next i

modX(2) = 1
argX(2) = 0
modX(6) = 0
argX(6) = 0

'  cos (90° – x) = sin x
'  -sin ( x-90) = cos x
'

' Hilbert(sin(x)) = -cos(t)
' Hilbert(Cos(x)) = sin(x)

' For i = 0 To n - 1
'   Xi1(i + 1) = 1 * Sin(2 * M_PI * 15.625 * i / fs + 0) + 0.5 * Cos(2 * M_PI * 93.75 * i / fs + 0)

'  Next i

Call SynthSignalfromRFi(modX, argX, Xi, Xq, n, 1)
Call printcomplex(Xi, Xq, n, 14, 15)

Call Getfassoc(fs, fv, n)
Call printarray(fv, n, 16)

Call GetSpectrumfromIQ(Xi, Xq, n, modX1, argX1, 1)
Call printcomplex(modX1, argX1, n, 17, 18)

For i = 0 To n - 1
Xi1(i + 1) = 0.5 * Cos(2 * M_PI * 250 * i / fs + 0) + 1 * Sin(2 * M_PI * 125 * i / fs + 0)

Next i
'Xi1(1) = 0

Call GetSpectrum(Xi1, n, modX1, argX1, 1)
Call printcomplex(modX1, argX1, n, 19, 20)

Call GetSpectrumIQ(Xi1, n, Yre, Yim)
Call printcomplex(Yre, Yim, n, 21, 22)

End Sub

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

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

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

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

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

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

'         column 7 to 8

' - 4. - 4.i  - 4. - 9.6568542i
Sub GetIQfromR_fi(ByVal Num As Long, ByRef modXk() As Double, ByRef argXk() As Double, _
ByRef Xkre() As Double, ByRef Xkim() As Double, ByVal mode As Long)
Dim i As Long

If mode = 1 Then
End If

For i = 1 To Num
Xkre(i) = modXk(i) * Cos(rad * argXk(i))
Xkim(i) = modXk(i) * Sin(rad * argXk(i))
Next i

End Sub

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

Dim i, k, N As Long

Dim angle As Double

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

End Sub

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

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

Sub IDFT(ByVal numidft As Long, ByRef idft_Xkre() As Double, _
ByRef idft_Xkim() As Double, ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
Dim k, N    As Long
For N = 0 To numidft - 1
'idft_item(n) = numidft - 1
idft_xnre(N + 1) = 0
idft_xnim(N + 1) = 0
For k = 0 To numidft - 1
angle = 2 * M_PI * k * N / numidft
idft_xnre(N + 1) = idft_xnre(N + 1) + idft_Xkre(k + 1) * Cos(angle) - idft_Xkim(k + 1) * Sin(angle)
idft_xnim(N + 1) = idft_xnim(N + 1) + idft_Xkre(k + 1) * Sin(angle) + idft_Xkim(k + 1) * Cos(angle)
Next k

idft_xnre(N + 1) = idft_xnre(N + 1) / numidft
idft_xnim(N + 1) = idft_xnim(N + 1) / numidft

Next N

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

Sub Hilbert(ByRef xre() As Double, ByVal Numb As Long, ByRef Hilbertre() As Double, ByRef Hilbertim() As Double)
'MATLAB algorithm
Dim k, i, no2 As Long
Dim h(1 To NMAX) As Double

' Dim dft_xnre(1 To NMAX)     As Double
Dim dft_xnim(1 To NMAX)     As Double
Dim dft_Xkre(1 To NMAX)    As Double
Dim dft_Xkim(1 To NMAX)    As Double
Dim idft_Xkre(1 To NMAX)   As Double
Dim idft_Xkim(1 To NMAX)   As Double
Dim idft_xnre(1 To NMAX)   As Double
Dim idft_xnim(1 To NMAX)   As Double

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

no2 = Int(Numb / 2)

For i = 1 To Numb

dft_xnim(i) = 0
Next i

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

i = 0
Do While i <= Numb

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

'if  N is even

If (Numb > 0) And (2 * no2) = Numb Then
h(1) = 1
h(no2 + 1) = 1
For k = 2 To no2
h(k) = 2
Next k

'if N is odd

ElseIf (Numb > 0) Then

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

End If

For i = 1 To Numb
idft_Xkre(i) = dft_Xkre(i) * h(i)
idft_Xkim(i) = dft_Xkim(i) * h(i)
Next i

Call IDFT(Numb, idft_Xkre, idft_Xkim, idft_xnre, idft_xnim)

For i = 1 To Numb
Hilbertre(i) = idft_xnre(i)
Hilbertim(i) = idft_xnim(i)
Next i

End Sub

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

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

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

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

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

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

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

Next i

End Sub

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

End Sub

Sub GetSpectrumfromIQ(ByRef x_re() As Double, ByRef x_im() As Double, _
ByVal Num As Long, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Long)

Dim Xkre(1 To NMAX) As Double
Dim Xkim(1 To NMAX) As Double

Call DFT(Num, x_re, x_im, Xkre, Xkim)
Call ScaleDFT(Num, Xkre, Xkim)
Call GetR_fi(Num, Xkre, Xkim, modXk, argXk, mode)

End Sub

Sub SynthSignalfromIQ(ByVal Num As Long, ByRef Xkre() As Double, ByRef Xkim() As Double, _
ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)

'Call ScaleIQidft(Num, Xkre, Xkim)
Call IDFT(Num, Xkre, Xkim, idft_xnre, idft_xnim)
Call ScaleIQidft(Num, idft_xnre, idft_xnim)
End Sub

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

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

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

Next i
End Sub

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

deltat = -(argF(N0) - argF(0)) / (2 * M_PI * fv(N0) - 2 * M_PI * fv(0))
End Function

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

Open FilePath For Output As #1

'Write #1,message
'Write #1,"\n ; ; ; ;   "
' If  mode=0  Then
'Write #1,"\n n ; Re(x) ; +j* Im(x)      ;   "
'  End If

'If  mode=1 Then
' Write #1,"\n n ; |X|   ; exp(j* Arg(x) );   "
'  End If

For i = 1 To N
Write #1, i
Write #1, xnre(i)
Write #1, xnim(i)
Next i

Write #1, "\n ; ; ; ;  "
Close #1
MsgBox (" data saved to file")

End Sub

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

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

Xq(1) = 0
Xq(2) = 0
Xq(3) = 0
Xq(4) = 0
Xq(5) = 0
Xq(6) = 0
Xq(7) = 0
Xq(8) = 0
fs = 1000
Call printcomplex(Xi, Xq, N, 1, 2)
Call Getfassoc(fs, fv, N)
Call Gettassoc(fs, tv, N)
Call printcomplex(tv, fv, N, 3, 4)

Call DFT(N, Xi, Xq, Yi, Yq)

Call printcomplex(Yi, Yq, N, 5, 6)

Call IDFT(N, Yi, Yq, Ii, Iq)
Call printcomplex(Ii, Iq, N, 7, 8)
Call Hilbert(Xi, N, Hi, Hq)
Call printcomplex(Hi, Hq, N, 9, 10)
Call GetSpectrum(Xi, N, modX, argX, 1)
Call printcomplex(modX, argX, N, 11, 12)

Call getenvelope(Xi, N, env)
Call printarray(env, N, 13)

N = 16

For i = 0 To N - 1
argX(i + 1) = 0
modX(i + 1) = 0
Next i

modX(2) = 1
argX(2) = -45
modX(6) = 0
argX(6) = 0

'  cos (90° – x) = sin x
'  -sin ( x-90) = cos x
'

' Hilbert(sin(x)) = -cos(t)
' Hilbert(Cos(x)) = sin(x)

' For i = 0 To n - 1
'   Xi1(i + 1) = 1 * Sin(2 * M_PI * 15.625 * i / fs + 0) + 0.5 * Cos(2 * M_PI * 93.75 * i / fs + 0)

'  Next i

Call SynthSignalfromRFi(modX, argX, Xi, Xq, N, 1)
Call printcomplex(Xi, Xq, N, 14, 15)

Call Getfassoc(fs, fv, N)
Call printarray(fv, N, 16)

Call GetSpectrumfromIQ(Xi, Xq, N, modX1, argX1, 1)
Call printcomplex(modX1, argX1, N, 17, 18)

For i = 0 To N - 1
Xi1(i + 1) = 0.5 * Cos(2 * M_PI * 250 * i / fs + 0) + 1 * Sin(2 * M_PI * 125 * i / fs + 0)

Next i
'Xi1(1) = 0

Call GetSpectrum(Xi1, N, modX1, argX1, 1)
Call printcomplex(modX1, argX1, N, 19, 20)

Call GetSpectrumIQ(Xi1, N, Yre, Yim)
Call printcomplex(Yre, Yim, N, 21, 22)

End Sub

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

Option Explicit
Const NMAX As Long = 65536

'  ByRef xnre() As Double, ByRef xnim() As Double, _

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

'Dim FilePath As String
Dim i As Long

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

str1 = ""
Open FilePath For Input As #1

Line Input #1, str1 ' Preamble          Input #1,str1
Line Input #1, str1 ' "Function name"

If f = 0 Then
Line Input #1, str1 '   ; ; ; 8 ;
word = Split(str1, ";")
N = CLng(word(3))
End If

If f = 1 Then
Input #1, str1, str1, str1, str1, str2  '   , , ,   8 ,

N = CLng(str1)
End If

Line Input #1, str1 ' "  n ; Re(x) ; +j* Im(x)      ;   " or "\n n ; |X|   ; exp(j* Arg(x) );   "

For i = 1 To N

If f = 0 Then

Line Input #1, str
word = Split(str, ";")
'  n1 = UBound(word)
str1 = word(0)
str2 = word(1)
str3 = word(2)
End If

If f = 1 Then

Input #1, str1, str2, str3

End If

ind(i) = CLng(Replace(str1, Comma, ","))
If f = 0 Then
ind(i) = ind(i) + 1
End If

xnre(i) = CDbl(Replace(str2, Comma, ","))
xnim(i) = CDbl(Replace(str3, Comma, ","))

'  Cells(i + 1, 1).Value = ind(i)
'   Cells(i + 1, 2).Value = xnre(i)
'  Cells(i + 1, 3).Value = xnim(i)

Next i

Line Input #1, str1   '  "\n ; ; ; ;  "
Close #1

End Sub

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

Dim i As Long

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

str1 = ""

Open FilePath For Input As #1

Line Input #1, str1 ' Preamble          Input #1,str1
Line Input #1, str1 ' "Function name"

If f = 0 Then   'default  format , which  uses  ';' as delim. symbol
Line Input #1, str1 '   ; ; ; 8 ;
word = Split(str1, ";")
N = CLng(word(3))
End If

If f = 1 Then '  uses  ' "string1","string2", "string3"' as format , internal
Input #1, str1, str1, str1, str1, str2  '   , , ,   8 ,

N = CLng(str1)
End If

Line Input #1, str1 ' "  n ; Re(x) ; +j* Im(x)      ;   " or "\n n ; |X|   ; exp(j* Arg(x) );   "

For i = 1 To N

If f = 0 Then

Line Input #1, str
word = Split(str, ";")
'  n1 = UBound(word)
str1 = word(0)
str2 = word(1)
' str3 = word(2)
End If

If f = 1 Then

Input #1, str1, str2 ', str3

End If

ind(i) = CLng(Replace(str1, Comma, ","))
If f = 0 Then
ind(i) = ind(i) + 1
End If

xnre(i) = CDbl(Replace(str2, Comma, ","))
'   xnim(i) = CDbl(Replace(str3, Comma, ","))

'   Cells(i + 1, 1).Value = ind(i)
'  Cells(i + 1, 2).Value = xnre(i)
'   Cells(i + 1, 3).Value = xnim(i)

Next i

Line Input #1, str1   '  "\n ; ; ; ;  "
Close #1

End Sub

Sub WriteToFileCplx(ByRef xnre() As Double, ByRef xnim() As Double, _
ByRef FilePath As String, ByVal N As Long, ByVal mode As Integer)
' Dim FilePath As String
Dim i As Long

Open FilePath For Output As #1

Write #1, , , , , " Data sequence     "
Write #1, , , , , " Begin of sequence "
Write #1, , , , N
If mode = 0 Then
Write #1, " n ", " Re(x) ", " +j* Im(x) "
End If

If mode = 1 Then
Write #1, " n ", " |X|  ", " Exp(j * Arg(x)) "

End If

For i = 1 To N

Write #1, i, xnre(i), xnim(i)

Next i

Write #1, , , , , "  End of sequence "
Close #1
MsgBox ("Saved , OK")

End Sub

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

Open FilePath For Output As #1

Write #1, , , , , " Data sequence     "
Write #1, , , , , " Begin of sequence "
Write #1, , , , N

Write #1, , , , , " Data  "

For i = 1 To N
Write #1, i, xn(i)
Next i

Write #1, , , , , "  End of sequence "
Close #1
MsgBox ("  Saved , OK")

End Sub

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

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

Sub PrintArray(ByRef xnre() As Double, ByVal num As Long, ByVal Col1 As Long)
Dim i As Long
For i = 1 To num
Cells(i + 1, Col1).Value = xnre(i)

Next i
End Sub

Sub GetArray(ByRef xnre() As Double, ByVal num As Long, ByVal Col1 As Long)
Dim i As Long
For i = 1 To num
xnre(i) = Cells(i + 1, Col1).Value

Next i
End Sub

Sub doInput()
Dim ind(1 To NMAX) As Long
Dim xre(1 To 100) As Double
Dim xim(1 To 100) As Double
Dim num As Long

Call ReadFromFileCplx("file21.csv", num, ind, xre, xim, 1, ".")
Call PrintCplx(xre, xim, num, 2, 3)

Call WriteToFileCplx(xre, xim, "file21.csv", 8, 0)

End Sub

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

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

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

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

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

|

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

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

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

'Dim FilePath As String
Dim i As Long

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

str1 = ""
Open FilePath For Input As #1

Line Input #1, str1 ' Preamble          Input #1,str1
Line Input #1, str1 ' "Function name"

If f = 0 Then
Line Input #1, str1 '   ; ; ; 8 ;
word = Split(str1, ";")
N = CLng(word(3))
End If

If f = 1 Then
Input #1, str1, str1, str1, str1, str2  '   , , ,   8 ,

N = CLng(str1)
End If

Line Input #1, str1 ' "  n ; Re(x) ; +j* Im(x)      ;   " or "\n n ; |X|   ; exp(j* Arg(x) );   "

For i = 1 To N

If f = 0 Then

Line Input #1, str
word = Split(str, ";")
'  n1 = UBound(word)
str1 = word(0)
str2 = word(1)
str3 = word(2)
End If

If f = 1 Then

Input #1, str1, str2, str3

End If

ind(i) = CLng(Replace(str1, Comma, ","))
If f = 0 Then
ind(i) = ind(i) + 1
End If

xnre(i) = CDbl(Replace(str2, Comma, ","))
xnim(i) = CDbl(Replace(str3, Comma, ","))

'  Cells(i + 1, 1).Value = ind(i)
'   Cells(i + 1, 2).Value = xnre(i)
'  Cells(i + 1, 3).Value = xnim(i)

Next i

Line Input #1, str1   '  "\n ; ; ; ;  "
Close #1

End Sub

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

Dim i As Long

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

str1 = ""

Open FilePath For Input As #1

Line Input #1, str1 ' Preamble          Input #1,str1
Line Input #1, str1 ' "Function name"

If f = 0 Then   'default  format , which  uses  ';' as delim. symbol
Line Input #1, str1 '   ; ; ; 8 ;
word = Split(str1, ";")
N = CLng(word(3))
End If

If f = 1 Then '  uses  ' "string1","string2", "string3"' as format , internal
Input #1, str1, str1, str1, str1, str2  '   , , ,   8 ,

N = CLng(str1)
End If

Line Input #1, str1 ' "  n ; Re(x) ; +j* Im(x)      ;   " or "\n n ; |X|   ; exp(j* Arg(x) );   "

For i = 1 To N

If f = 0 Then

Line Input #1, str
word = Split(str, ";")
'  n1 = UBound(word)
str1 = word(0)
str2 = word(1)
' str3 = word(2)
End If

If f = 1 Then

Input #1, str1, str2 ', str3

End If

ind(i) = CLng(Replace(str1, Comma, ","))
If f = 0 Then
ind(i) = ind(i) + 1
End If

xnre(i) = CDbl(Replace(str2, Comma, ","))
'   xnim(i) = CDbl(Replace(str3, Comma, ","))

'   Cells(i + 1, 1).Value = ind(i)
'  Cells(i + 1, 2).Value = xnre(i)
'   Cells(i + 1, 3).Value = xnim(i)

Next i

Line Input #1, str1   '  "\n ; ; ; ;  "
Close #1

End Sub

Sub WriteToFileCplx(ByRef xnre() As Double, ByRef xnim() As Double, _
ByRef FilePath As String, ByVal N As Long, ByVal mode As Integer)
' Dim FilePath As String
Dim i As Long

Open FilePath For Output As #1

Write #1, , , , , " Data sequence     "
Write #1, , , , , " Begin of sequence "
Write #1, , , , N
If mode = 0 Then
Write #1, " n ", " Re(x) ", " +j* Im(x) "
End If

If mode = 1 Then
Write #1, " n ", " |X|  ", " Exp(j * Arg(x)) "

End If

For i = 1 To N

Write #1, i, xnre(i), xnim(i)

Next i

Write #1, , , , , "  End of sequence "
Close #1
MsgBox ("Saved , OK")

End Sub

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

Open FilePath For Output As #1

Write #1, , , , , " Data sequence     "
Write #1, , , , , " Begin of sequence "
Write #1, , , , N

Write #1, , , , , " Data  "

For i = 1 To N
Write #1, i, xn(i)
Next i

Write #1, , , , , "  End of sequence "
Close #1
MsgBox ("  Saved , OK")

End Sub

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

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

Sub PrintArray(ByRef xnre() As Double, ByVal num As Long, ByVal Col1 As Long)
Dim i As Long
For i = 1 To num
Cells(i + 1, Col1).Value = xnre(i)

Next i
End Sub

Sub GetArray(ByRef xn() As Double, ByVal num As Long, ByVal Col1 As Long)
Dim i As Long
For i = 1 To num
xn(i) = Cells(i + 1, Col1).Value

Next i
End Sub

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

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

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

End Function

Sub GetIQfromR_fi(ByVal num As Long, ByRef modXk() As Double, ByRef argXk() As Double, _
ByRef Xkre() As Double, ByRef Xkim() As Double, ByVal mode As Long)
Dim i As Long

If mode = 1 Then
End If

For i = 1 To num
Xkre(i) = modXk(i) * Cos(rad * argXk(i))
Xkim(i) = modXk(i) * Sin(rad * argXk(i))
Next i

End Sub

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

Dim i, k, N As Long

Dim angle As Double

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

End Sub

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

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

Sub CIDFT(ByVal numidft As Long, ByRef idft_Xkre() As Double, _
ByRef idft_Xkim() As Double, ByRef idft_xnre() As Double, ByRef idft_xnim() As Double)
Dim k, N    As Long
Dim angle As Double

For N = 0 To numidft - 1
'idft_item(n) = numidft - 1
idft_xnre(N + 1) = 0
idft_xnim(N + 1) = 0
For k = 0 To numidft - 1
angle = 2 * M_PI * k * N / numidft
idft_xnre(N + 1) = idft_xnre(N + 1) + idft_Xkre(k + 1) * Cos(angle) - idft_Xkim(k + 1) * Sin(angle)
idft_xnim(N + 1) = idft_xnim(N + 1) + idft_Xkre(k + 1) * Sin(angle) + idft_Xkim(k + 1) * Cos(angle)
Next k

idft_xnre(N + 1) = idft_xnre(N + 1) / numidft
idft_xnim(N + 1) = idft_xnim(N + 1) / numidft

Next N

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

Sub Hilbert(ByRef xre() As Double, ByVal Numb As Long, ByRef Hilbertre() As Double, ByRef Hilbertim() As Double)
'MATLAB algorithm
Dim k, i, no2 As Long
Dim h(1 To NMAX) As Double

' Dim dft_xnre(1 To NMAX)     As Double
Dim dft_xnim(1 To NMAX)     As Double
Dim dft_Xkre(1 To NMAX)    As Double
Dim dft_Xkim(1 To NMAX)    As Double
Dim idft_Xkre(1 To NMAX)   As Double
Dim idft_Xkim(1 To NMAX)   As Double
Dim idft_xnre(1 To NMAX)   As Double
Dim idft_xnim(1 To NMAX)   As Double

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

no2 = Int(Numb / 2)

For i = 1 To Numb

dft_xnim(i) = 0
Next i

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

i = 0
Do While i <= Numb

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

'if  N is even

If (Numb > 0) And (2 * no2) = Numb Then
h(1) = 1
h(no2 + 1) = 1
For k = 2 To no2
h(k) = 2
Next k

'if N is odd

ElseIf (Numb > 0) Then

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

End If

For i = 1 To Numb
idft_Xkre(i) = dft_Xkre(i) * h(i)
idft_Xkim(i) = dft_Xkim(i) * h(i)
Next i

Call CIDFT(Numb, idft_Xkre, idft_Xkim, idft_xnre, idft_xnim)

For i = 1 To Numb
Hilbertre(i) = idft_xnre(i)
Hilbertim(i) = idft_xnim(i)
Next i

End Sub

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

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

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

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

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

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

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

Next i

End Sub

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

End Sub

Sub GetSpectrumfromIQ(ByRef x_re() As Double, ByRef x_im() As Double, _
ByVal num As Long, ByRef modXk() As Double, ByRef argXk() As Double, ByVal mode As Long)

Dim Xkre(1 To NMAX) As Double
Dim Xkim(1 To NMAX) As Double

Call CDFT(num, x_re, x_im, Xkre, Xkim)
Call ScaleCDFT(num, Xkre, Xkim)
Call GetR_fi(num, Xkre, Xkim, modXk, argXk, mode)

End Sub

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

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

Sub GetMFCoesffs(ByRef xn() As Double, ByVal num As Long, ByRef Zasre() As Double, ByRef Zasim() As Double)

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

End Sub

Sub GetMatchedFilter(ByRef xn() As Double, ByVal num As Long, ByRef Zasre() As Double, ByRef Zasim() As Double, _
ByRef xoutre() As Double, ByRef xoutim() As Double)

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

End Sub

Sub GetRejectMatchedFilter(ByRef xn() As Double, ByVal num As Long, ByRef Zasre() As Double, ByRef Zasim() As Double, _
ByRef xoutre() As Double, ByRef xoutim() As Double)

Dim xre(1 To NMAX) As Double
Dim xim(1 To NMAX) As Double
Dim yre(1 To NMAX) As Double
Dim yim(1 To NMAX) As Double
Dim i As Long
Call GetAnalyticSignal(xn, num, xre, xim)
Call GetMatchedFilter(xn, num, Zasre, Zasim, yre, yim)

For i = 1 To num
xoutre(i) = xre(i) - yre(i) ' or another  expression
xoutim(i) = xim(i) - yim(i) ' or another  expression

Next i

End Sub

'*************************************************

Sub doInput()
Dim ind(1 To NMAX) As Long
Dim xre(1 To 100) As Double
Dim xim(1 To 100) As Double
Dim num As Long

Call ReadFromFileCplx("file21.csv", num, ind, xre, xim, 1, ".")
Call PrintCplx(xre, xim, num, 2, 3)

Call WriteToFileCplx(xre, xim, "file21.csv", 8, 0)

End Sub

Sub GetSignalMF()
Dim ModXkre(1 To NMAX) As Double
Dim argXk(1 To NMAX)  As Double
Dim xnre(1 To NMAX) As Double
Dim xnim(1 To NMAX)  As Double
Dim num As Long
num = Cells(1, 1).Value

Call GetArray(ModXkre, num, 2)
Call GetArray(argXk, num, 3)

Call SynthSignalfromRFi(ModXkre, argXk, xnre, xnim, num, 0)
Call PrintCplx(xnre, xnim, num, 4, 5)

End Sub

Sub CoeffsMF()
Dim xnre(1 To NMAX) As Double
Dim xnim(1 To NMAX)  As Double
Dim Zasre(1 To NMAX) As Double
Dim Zasim(1 To NMAX) As Double
Dim num As Long
num = Cells(1, 1).Value
Call GetArray(xnre, num, 4)
'Call GetArray(xnim, num, 5)
Call GetMFCoesffs(xnre, num, Zasre, Zasim)
Call PrintCplx(Zasre, Zasim, num, 7, 8)
End Sub

Sub TestMF()
Dim xnre(1 To NMAX) As Double
Dim xnim(1 To NMAX)  As Double
Dim Zasre(1 To NMAX) As Double
Dim Zasim(1 To NMAX) As Double
Dim xre(1 To NMAX) As Double
Dim xim(1 To NMAX) As Double
Dim env(1 To NMAX) As Double
Dim num As Long
num = Cells(1, 1).Value
Call GetArray(xnre, num, 4)
'Call GetArray(xnim, num, 5)
Call GetArray(Zasre, num, 7)
Call GetArray(Zasim, num, 8)
Call GetMatchedFilter(xnre, num, Zasre, Zasim, xre, xim)

Call PrintCplx(xre, xim, num, 10, 11)
Call GetR_fi(num, xre, xim, xnre, xnim, 1)
'Call GetEnvelope(xnre, num, env)

Call PrintArray(xnre, num, 12)

End Sub

Sub TestRMF()
Dim xnre(1 To NMAX) As Double
Dim xnim(1 To NMAX)  As Double
Dim Zasre(1 To NMAX) As Double
Dim Zasim(1 To NMAX) As Double
Dim xre(1 To NMAX) As Double
Dim xim(1 To NMAX) As Double
Dim env(1 To NMAX) As Double
Dim num As Long
num = Cells(1, 1).Value
Call GetArray(xnre, num, 4)
'Call GetArray(xnim, num, 5)
Call GetArray(Zasre, num, 7)
Call GetArray(Zasim, num, 8)
Call GetRejectMatchedFilter(xnre, num, Zasre, Zasim, xre, xim)

Call PrintCplx(xre, xim, num, 13, 14)
Call GetR_fi(num, xre, xim, xnre, xnim, 1)
Call PrintArray(xnre, num, 15)
End Sub

• Edited by 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