Asked by:
FIR and spectrum tools for VC2005 Express

General discussion
-
afrtohtcoefs -program for converting AFR coeffs to h(t) coeffs for FIR
// afrtohtcoefs.cpp : Defines the entry point for the console application.
//#include "stdafx.h"
#include <stdio.h>
#include <tchar.h>
#include <iostream>
#include <stdlib.h>
#include <cmath>#define M_PI 3.1415926535897932384626433832795
//int _tmain(int argc, _TCHAR* argv[])
//{
// return 0;
//}//afrtohtcoefs.cpp
#define NMAX 512
FILE *fp, *fp1;
int k,n;
int w,t,N;
double h[NMAX];
double h_im[NMAX];
double Re_H[NMAX];
double Im_H[NMAX];
double Re_x[NMAX];
double Im_x[NMAX];
double abs_H[NMAX];
double arg_H[NMAX];
double f[NMAX];
float angle;
double K, fs,T,finp,Tmax;
char str;
char key;/*
float heavi( double t)
{
if (t<0 ) return 0;
if (t>=0) return 1;
}float heavi1( int n)
{
if (n<0 ) return 0;
if (n>=0) return 1;
}*/double getarg(float Re_H,float Im_H)
{return atan2(Im_H,Re_H) * 180 / M_PI;
}
int _tmain(int argc, _TCHAR* argv[])
{
printf ("\n For input |H(jw)|, fi(w) from file input '1' ");
printf ("\n For manual input |H(jw)|, fi(w) input '0' ");//std::cin>>key;
scanf( "%s", &key);
if (key=='0' ) { goto label_02;}
if (key=='1' ) { goto label_01;}label_01:
printf("\n Open afrcoefs.txt for reading AFR coefficients ");
if ((fp=fopen("afrcoefs.txt","r"))==NULL) printf ("\n File opening error");
fscanf(fp,"%le",&fs);
fscanf(fp,"%d",&N);
Tmax=N*1/fs;
for (k=0 ; k<N;k++)
{
// f[k]=k*fs/N;
//fprintf(fp1,"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
fscanf (fp,"%d",&n);
fscanf (fp,"%le",&abs_H[k]);
fscanf (fp,"%le",&arg_H[k]);
fscanf (fp,"%le",&f[k]);
arg_H[k]=arg_H[k]*M_PI/180;
}
fclose(fp);
goto label_03;label_02:
printf("\n Input fs ,Hz ");
scanf("%le",&fs);
printf(" Input number of points ");
scanf("%d",&N);
Tmax=N*1/fs;
printf("\n Tmax =%lf s", Tmax);
for (k=0 ; k<N;k++)
{
f[k]=k*fs/N;//fprintf(fp1,"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
printf("\n F=%le Hz" ,f[k]);
printf("\n k=%d; |H(jw)|[k]=",k );
scanf ( "%le",&abs_H[k]);
printf(" k=%d; arg(H(jw))[k]=",k );
scanf ( "%le",&arg_H[k]);
//abs_H[k]=1;
//arg_H[k]=0;
arg_H[k]=arg_H[k]*M_PI/180;
}label_03:
for (k=0 ; k<N;k++)
{
f[k]=k*fs/N;
printf("\n n=%d;|H(jw)|= %le ;arg(H(jw))= %le deg;f=%lf Hz",k, abs_H[k] , arg_H[k]*180/M_PI , f[k]) ;
}
printf("\n Input 0 for continue ");
scanf("%s",&key);
printf("\n OK. Parsing ...");
// H(jw)=integr(0;inf;h(tau)*exp(-j*w*tau) ; d tau)
//K=2*3.1415926/fs;
/*
DFT
X[k]=sum(n=0;N-1; x[n]*exp(-2*pi*j*k*n/N) ), k=0,...,N-1
IDFT
x[n]=(1/N)*sum(k=0;N-1;X[k]*exp(2*pi*k*n/N) ), n=0,...,N-1
100*10^-6*500
*/label1:;
//cidft
printf("\n Open hcoefs1.txt for writing coeffs h[n] ");
if ((fp1=fopen("hcoefs1.txt","a"))==NULL) printf ("\n File opening error");
fprintf(fp1,"\n\n %le %d", fs, N );
for (n=0 ; n<N; n++)
{
Re_x[n]=0;
Im_x[n]=0;
for(k=0;k<N;k++)
{
angle=2*M_PI*k*n/N;
//Uoutp[w]/Uinp[w] exp(j( fi out[w]-fi inp[w]) )
Re_H[k]=abs_H[k]*cos(arg_H[k]);
Im_H[k]=abs_H[k]*sin(arg_H[k]);
Re_x[n] = Re_x[n]+ Re_H[k]*cos(angle)-Im_H[k]*sin(angle);
Im_x[n] = Im_x[n]+ Re_H[k]*sin(angle) + Im_H[k]*cos(angle);
}
Re_x[n]=Re_x[n]/N ;
Im_x[n]=Im_x[n]/N ;
// n scaled CIDFT , not used in this, use if sygnal synth.
//Re_x[n]=Re_x[n]*N;
//Im_x[n]=Im_x[n]*N;
//printf("\n n=%d Re(h(n))=%lf Im(h(n))=%lf ",n, Re_x[n] ,Im_x[n] ) ;
printf("\n n=%d h[n]=%lf ",n, Re_x[n] ) ;
fprintf(fp1, "\n %d %le ",n, Re_x[n] ) ;
}fclose (fp1);
printf( "\n Re(h(n)=h[n]), use it to load into FIR " ) ;
label3:;
printf("\n Input 0 for exit ");
scanf("%s",&key);
return 0;
}// stdafx.cpp : source file that includes just the standard includes
// afrtohtcoefs.pch will be the pre-compiled header
// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"
// TODO: reference any additional headers you need in STDAFX.H
// and not in this file// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//#pragma once
#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers#include <stdio.h>
#include <tchar.h>
#include <iostream>
#include <stdlib.h>
#include <cmath>// TODO: reference additional headers your program requires here
afrcoefs.txt
4.410000e+004 16
0 1.000000 0.000000 0.000000
1 1.000000 -168.750000 2756.250000
2 1.000000 22.500000 5512.500000
3 1.000000 -146.249999 8268.750000
4 0.000000 45.000005 11025.000000
5 0.000000 -123.750001 13781.250000
6 0.000000 67.500002 16537.500000
7 0.000000 78.750002 19293.750000
8 0.000000 90.000002 22050.000000
9 0.000000 -78.749998 24806.250000
10 0.000000 -67.499999 27562.500000
11 0.000000 123.750002 30318.750000
12 0.000000 135.000007 33075.000000
13 1.000000 146.250002 35831.250000
14 1.000000 -22.499997 38587.500000
15 1.000000 168.750003 41343.750000
4.410000e+004 16
0 1.000001e+000 0.000000e+000 0.000000e+000
1 1.000000e+000 -1.687500e+002 2.756250e+003
2 9.999998e-001 2.250004e+001 5.512500e+003
3 1.000000e+000 -1.462500e+002 8.268750e+003
4 3.961544e-007 -2.516211e+001 1.102500e+004
5 5.013260e-007 5.027288e+001 1.378125e+004
6 3.873424e-007 1.400188e+002 1.653750e+004
7 4.839081e-007 1.122935e+002 1.929375e+004
8 3.130084e-007 1.795795e+002 2.205000e+004
9 4.799148e-007 -1.126101e+002 2.480625e+004
10 3.857124e-007 -1.408047e+002 2.756250e+004
11 4.938736e-007 -5.018234e+001 3.031875e+004
12 3.914724e-007 2.715836e+001 3.307500e+004
13 1.000000e+000 1.462500e+002 3.583125e+004
14 9.999998e-001 -2.250004e+001 3.858750e+004
15 1.000000e+000 1.687500e+002 4.134375e+004
4.410000e+004 16
0 1.000001e+000 0.000000e+000 0.000000e+000
1 1.000000e+000 -1.687500e+002 2.756250e+003
2 9.999998e-001 2.250004e+001 5.512500e+003
3 1.000000e+000 -1.462500e+002 8.268750e+003
4 3.949392e-007 -2.565670e+001 1.102500e+004
5 4.989970e-007 5.024488e+001 1.378125e+004
6 3.867226e-007 1.403127e+002 1.653750e+004
7 4.821592e-007 1.124314e+002 1.929375e+004
8 3.130000e-007 1.800000e+002 2.205000e+004
9 4.821592e-007 -1.124314e+002 2.480625e+004
10 3.867226e-007 -1.403127e+002 2.756250e+004
11 4.989970e-007 -5.024488e+001 3.031875e+004
12 3.949392e-007 2.565670e+001 3.307500e+004
13 1.000000e+000 1.462500e+002 3.583125e+004
14 9.999998e-001 -2.250004e+001 3.858750e+004
15 1.000000e+000 1.687500e+002 4.134375e+004hcoefs1.txt
4.410000e+004 16
0 -4.854692e-002
1 3.078799e-002
2 6.781642e-002
3 -7.925028e-003
4 -9.804488e-002
5 -3.848729e-002
6 1.898830e-001
7 4.045169e-001
8 4.045168e-001
9 1.898829e-001
10 -3.848736e-002
11 -9.804485e-002
12 -7.924644e-003
13 6.781681e-002
14 3.078794e-002
15 -4.854676e-002- Edited by USERPC01 Tuesday, January 23, 2018 10:03 PM adding afrcoefs.txt
Tuesday, January 23, 2018 9:58 PM
All replies
-
Program for converting h(t) coefficients of the FIR to AFR data
// irtoafr.cpp : Defines the entry point for the console application.
//#include "stdafx.h"
#include <cstdio> // or cstdio.h
#include <cstdlib> //or cstdlib.h
#include <cmath> //or cmath.h#define M_PI 3.1415926535897932384626433832795
#define PI M_PI // 3.14159265
#define NMAX 512
FILE *fp, *fp1;
int w,t,N;
double h[NMAX];
double h_im[NMAX];
double Re_H[NMAX];
double Im_H[NMAX];
double Re_x[NMAX];
double Im_x[NMAX];
double abs_H[NMAX];
double arg_h[NMAX];
double f[NMAX];
double angle;
double K, fs,T,Tmax;
char str;
char key;
int k,n;
/*
float x[N];
float y[N];double my_FIR(double sample_data, int Nf)
{
double result = 0;
for ( int i = Nf - 2 ; i >= 0 ; i-- )
{
x[i + 1] = x[i];
y[i + 1] = y[i];
}
x[0] = (double)sample_data;
for (int k = 0; k < N; k++)
{
result = result + x[k]*h[k];
}
y[0] = result;
return ((double)result);
}*/
float heavi( double t)
{
if (t<0 ) return 0;
if (t>=0) return 1;
}float heavi1( int n)
{
if (n<0 ) return 0;
if (n>=0) return 1;
}
float dirac( int n)
{
if (n!=0 ) return 0;
if (n=0) return 1;
}double getarg(float Re_H,float Im_H)
{
return atan2(Im_H,Re_H) * 180/ PI;
}
int _tmain(int argc, _TCHAR* argv[])
{
//double tau=0.02;
//printf(" Input tau , s ");
//scanf("%le",&tau);
printf ("\n For input h(t) from file input '1' ");
printf ("\n For manual input h(t) input '0' ");
scanf("%s",&key);
if (key=='0' ) { goto label_02;}
if (key=='1' ) { goto label_01;}label_01:
printf("\n Open htkoefs.txt for reading h[n] coefficients ");
if ((fp=fopen("htkoefs.txt","r"))==NULL) printf ("\n File opening error");
fscanf(fp,"%le",&fs);
fscanf(fp,"%d",&N);
Tmax=N*1/fs;
for (n=0 ; n<N;n++)
{
fscanf (fp,"%d",&k);
fscanf (fp,"%le",&h[n]);
}
fclose(fp);
goto label_03;
label_02:
printf("\n Input fs ,Hz ");
scanf("%le",&fs);
printf(" Input number of points ");
scanf("%d",&N);
Tmax=N*1/fs;
printf("\n Tmax =%lf s \n", Tmax);
for (n=0 ; n<N;n++)
{
T=n*1/(N*fs);
printf("\nT=%le s; n =%d ; input h(n): ",T,n);
scanf ("%le",&h[n]);/*
h[0]=-0.0485469;
h[1]=+0.0307880;
h[2]=+0.0678165;
h[3]=-0.0079250;
h[4]=-0.0980449;
h[5]=-0.0384873;
h[6]=+0.1898828;
h[7]=+0.4045168;
h[8]=+0.4045168;
h[9]=+0.1898828;
h[10]=-0.0384873;
h[11]=-0.0980449;
h[12]=-0.0079250;
h[13]=+0.0678165;
h[14]=+0.0307880;
h[15]=-0.0485469;
*/
//tau=0.1*Tmax;
//h[t]=1-exp(-T/tau);
}label_03:
printf("\n Fs=%lf Hz",fs);
printf("\n N=%d points",N);
for (n=0 ; n<N;n++)
{
T=n*1/(N*fs);
printf("\n t=%le s; h[t]=%le ",T,h[n]);
}printf("\n Input 0 for next ");
scanf("%s",&key);
printf("\n OK. Parsing ... ");
// H(jw)=integr(0;inf;h(tau)*exp(-j*w*tau) ; d tau)
//
/*
DFT
X[k]=sum(n=0;N-1; x[n]*exp(-2*pi*j*k*n/N) ), k=0,...,N-1
IDFT
x[n]=(1/N)*sum(k=0;N-1;X[k]*exp(2*pi*k*n/N) ), n=0,...,N-1
*//*********************************************************/
label1:
printf("\n Open afrcoefs.txt for reading AFR coefficients ");
if ((fp1=fopen("afrcoefs.txt","a"))==NULL) printf ("\n File opening error");
fprintf(fp1,"\n %le %d ",fs,N);
/*
fprintf(fp1,"\n Re(H(jw)) ; Im(H(jw)) ");
for (k=0 ; k<N; k++)
{
Re_H[k]=0;
Im_H[k]=0;
for(n=0;n<N;n++)
{
angle=-2*PI*k*n/N;
// Re_H[k] = Re_H[k]+ h[n]*cos(angle) - h_im[n]*sin(angle);
// Im_H[k] = Im_H[k]+ h[n]*sin(angle) + h_im[n]*cos(angle);
Re_H[k] = Re_H[k]+ h[n]*cos(angle) ;//- h_im[n]*sin(angle);
Im_H[k] = Im_H[k]+ h[n]*sin(angle) ;//+ h_im[n]*cos(angle);
}// 1/n scaled DFT(h(t)*1(t))
//Re_H[k]=Re_H[k]/N ;
//Im_H[k]=Im_H[k]/N ;
abs_H[k]=pow(((Re_H[k]*Re_H[k])+(Im_H[k]*Im_H[k])), 0.5);
arg_h[k]=getarg(Re_H[k],Im_H[k]);
f[k]=k*fs/N;
printf("\n n=%d Re(H(jw))= %lf ; Im(H(jw))= %lf ; f=%lf Hz",k, Re_H[k] ,Im_H[k], f[k]) ;
fprintf(fp1,"\n n=%d Re(H(jw))= %lf ; Im(H(jw))= %lf ; f=%lf Hz",k, Re_H[k] ,Im_H[k], f[k]) ;
}printf("\n Input 0 for continue ");
scanf("%s",&key);
*/
//fprintf(fp1,"\n |H(jw)| , arg( H(jw) ) ");
for (k=0 ; k<N; k++)
{
Re_H[k]=0;
Im_H[k]=0;
for(n=0;n<N;n++)
{
angle=-2*PI*k*n/N;
// Re_H[k] = Re_H[k]+ h[n]*cos(angle) - h_im[n]*sin(angle);
// Im_H[k] = Im_H[k]+ h[n]*sin(angle) + h_im[n]*cos(angle);
Re_H[k] = Re_H[k]+ h[n]*cos(angle) ;//- h_im[n]*sin(angle);
Im_H[k] = Im_H[k]+ h[n]*sin(angle) ;//+ h_im[n]*cos(angle);
}// 1/n scaled DFT(h(t)*1(t))
//Re_H[k]=Re_H[k]/N ;
//Im_H[k]=Im_H[k]/N ;
abs_H[k]=pow(((Re_H[k]*Re_H[k])+(Im_H[k]*Im_H[k])), 0.5);
arg_h[k]=getarg(Re_H[k],Im_H[k]);
f[k]=k*fs/N;
printf("\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
//fprintf(fp1,"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
fprintf(fp1,"\n %d %le %le %le ",k, abs_H[k] ,arg_h[k], f[k]) ;
}
fclose(fp1);
printf("\n Input 0 for exit ");
scanf("%s",&key);return 0;
}//cidft
/*
for (n=0 ; n<N; n++)
{
Re_x[n]=0;
Im_x[n]=0;
for(k=0;k<N;n++)
{
// angle = 2 * PI * k * N / numidft
// idft_xnre(N + 1) = idft_xnre(N + 1) + idft_Xkre(k + 1) * Cos(angle) - idft_Xkim(k + 1) * Sin(angle)
// idft_xnim(N + 1) = idft_xnim(N + 1) + idft_Xkre(k + 1) * Sin(angle) + idft_Xkim(k + 1) * Cos(angle)
angle=2*PI*k*n/N;
Re_H[k]=abs_H[k]*cos(angle);
Im_H[k]=abs_H[k]*sin(angle);
Re_x[n] = Re_x[n]+ Re_H[k]*cos(angle)- Im_H[k]*sin(angle);
Im_x[n] = Im_x[n]+ Re_H[k]*sin(angle)+ Im_H[k]*cos(angle);
}
Re_x[n]=Re_x[n]/N ;
Im_x[n]=Im_x[n]/N ;
// n scaled CIDFT( )
//Re_x[n]=Re_x[n]*N;
//Im_x[n]=Im_x[n]*N;
// abs_H[k]=pow(((Re_H[k]*Re_H[k])+(Im_H[k]*Im_H[k])), 0.5);
// arg_h[k]=getarg(Re_H[k],Im_H[k]);
// T[k]=k*fs/N;
printf("\n n=%d Re(h(n))=%d ",n, Re_x[n] ,Im_x[n] ) ;}
*/
// stdafx.cpp : source file that includes just the standard includes
// irtoafr.pch will be the pre-compiled header
// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"
// TODO: reference any additional headers you need in STDAFX.H
// and not in this file// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//#pragma once
#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
#include <stdio.h>
#include <tchar.h>// TODO: reference additional headers your program requires here
Tuesday, January 23, 2018 10:01 PM -
Implementation of the FIR
// fir.cpp : Defines the entry point for the console application.
//#include "stdafx.h"
#include <cstdio> // or cstdio.h
#include <cstdlib> //or cstdlib.h
#include <cmath> //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 < Nf; 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 _tmain(int argc, _TCHAR* argv[])
{
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;
}// stdafx.cpp : source file that includes just the standard includes
// fir.pch will be the pre-compiled header
// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"
// TODO: reference any additional headers you need in STDAFX.H
// and not in this file// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//#pragma once
#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
#include <stdio.h>
#include <tchar.h>// TODO: reference additional headers your program requires here
input.txt
4.410000e+004 160 1.000000e+001
1 7.071068e+000
2 6.123032e-016
3 -7.071068e+000
4 -1.000000e+001
5 -7.071068e+000
6 -1.836910e-015
7 7.071068e+000
8 1.000000e+001
9 7.071068e+000
10 3.061516e-015
11 -7.071068e+000
12 -1.000000e+001
13 -7.071068e+000
14 -4.286122e-015
15 7.071068e+000htkoefs.txt
4.410000e+004 16
0 0
1 1
2 0
3 0
4 0
5 0
6 0
7 0
8 0
9 0
10 0
11 0
12 0
13 0
14 0
15 0output.txt
4.410000e+004 160 0.000000e+000
1 1.000000e+001
2 7.071068e+000
3 6.123032e-016
4 -7.071068e+000
5 -1.000000e+001
6 -7.071068e+000
7 -1.836910e-015
8 7.071068e+000
9 1.000000e+001
10 7.071068e+000
11 3.061516e-015
12 -7.071068e+000
13 -1.000000e+001
14 -7.071068e+000
15 -4.286122e-015- Edited by USERPC01 Tuesday, January 23, 2018 10:08 PM fixed bugs
Tuesday, January 23, 2018 10:07 PM -
Program for spectrum calculating, using CDFT,CIDFT, Mathlab's algorithm for Hilbert transform
You can use complex.h, vector.h for data arrays and pointered arrays for arrays.
How to use radix more then 20000 in #define NMAX , how to fix memory mapping error in arrays if NMAX>20000...32768?
// getspec.cpp : Defines the entry point for the console application.
//#include "stdafx.h"
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <limits>
#include <iostream>
using namespace std;#define NMAX 2048
#define M_PI 3.1415926535897932384626433832795
typedef struct {
double re[NMAX];
double im[NMAX];
} CData;typedef struct {
double data[NMAX];
} Array;
CData cdft(CData xn ,int N )
{
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 ((double)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((double)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 _tmain(int argc, _TCHAR* argv[])
{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] : ");
if ((fp2=fopen("spectrum.txt","a"))==NULL) printf ("\n File opening error");
fprintf(fp2,"\n %le %d \n",fs,N1);
//fprintf(fp2, "\n n;\t f, Hz;\t |Xf(jw)| ;\t arg(Xf), deg ");
for (i=0; i<N1;i++)
{
f=fs*i/N1;
//cout <<"\n f="<< f<<" Hz;|Xf(jw)|="<<Xf.re[i]<<";arg(Xf)="<<Xf.im[i]<<" deg ";
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 ");
char key1;
cin>>key1;return 0;
}
// stdafx.cpp : source file that includes just the standard includes
// getspec.pch will be the pre-compiled header
// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"
// TODO: reference any additional headers you need in STDAFX.H
// and not in this file// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//#pragma once
#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
#include <stdio.h> //edit
#include <tchar.h> //edit// TODO: reference additional headers your program requires here
xinput.txt
4.410000e+004 16
0 1.000000e+001
1 7.071068e+000
2 6.123032e-016
3 -7.071068e+000
4 -1.000000e+001
5 -7.071068e+000
6 -1.836910e-015
7 7.071068e+000
8 1.000000e+001
9 7.071068e+000
10 3.061516e-015
11 -7.071068e+000
12 -1.000000e+001
13 -7.071068e+000
14 -4.286122e-015
15 7.071068e+000spectrum.txt
4.410000e+004 16
0 7.216450e-016 -1.126199e+002
1 7.237769e-016 -8.560129e+001
2 1.000000e+001 9.223608e-015
3 9.436896e-016 1.800000e+002
4 3.373574e-016 3.464047e+001
5 5.950643e-015 -5.660642e+001
6 1.330312e-007 1.800000e+002
7 7.196406e-015 -2.652734e+000
8 3.760637e-015 -7.363206e+001
9 3.624946e-015 -1.272209e+002
10 2.820163e-015 1.438068e+002
11 5.525033e-015 -1.015922e+002
12 4.814310e-015 7.423349e+000
13 1.110917e-014 6.612472e+001
14 6.445273e-015 -9.246812e+001
15 8.760391e-015 -1.326889e+002
- Edited by USERPC01 Tuesday, January 23, 2018 10:29 PM adding xinput.txt
Tuesday, January 23, 2018 10:14 PM -
Program for signal synth. from freq. data
// synthsig.cpp : Defines the entry point for the console application.
//#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>#define NMAX 1000
#define M_PI 3.1415926535897932384626433832795typedef struct {
double re[NMAX];
double im[NMAX];
} CData;typedef struct {
double data[NMAX];
} Array;
CData cdft(CData xn ,int N )
{
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 ((double)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((double) 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 ()
int _tmain(int argc, _TCHAR* argv[])
{double xin[NMAX];
CData Xf;
CData xt;
char key;
int k,n,i,N1;
double fs=44100, 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 ");
std::cin>>key;
return 0;}
// stdafx.cpp : source file that includes just the standard includes
// synthsig.pch will be the pre-compiled header
// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"
// TODO: reference any additional headers you need in STDAFX.H
// and not in this file// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//#pragma once
#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
#include <stdio.h>
#include <tchar.h>// TODO: reference additional headers your program requires here
spectrum.txt
44100 16
0 0 0
1 0 0
2 10 0
3 0 0
4 0 0
5 0 0
6 0 0
7 0 0
8 0 0
9 0 0
10 0 0
11 0 0
12 0 0
13 0 0
14 0 0
15 0 0Friday, January 26, 2018 1:27 PM -
// orthosynth.cpp : Defines the entry point for the console application.
//#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>#define NMAX 1000
#define M_PI 3.1415926535897932384626433832795typedef struct {
double re[NMAX];
double im[NMAX];
} CData;typedef struct {
double data[NMAX];
} Array;
CData cdft(CData xn ,int N )
{
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 ((double)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((double)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 _tmain(int argc, _TCHAR* argv[])
int main ()
{double xin[NMAX];
CData Xf;
CData xt;int k,n,i,N1;
double fs=44100,T,x,f;
int key=0;printf ("\n For input x(t) from file input '1' ");
printf ("\n For manual input x(t) input '0' ");
std::cin>>key;
switch(key)
{
case 1 :
goto label_05;
break;
case 0 :
goto label_06;
break;
}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 ");
char key1;
std::cin>>key1;
std::cout << key1;
return 0;}
// stdafx.cpp : source file that includes just the standard includes
// orthosynth.pch will be the pre-compiled header
// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"
// TODO: reference any additional headers you need in STDAFX.H
// and not in this file// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//#pragma once
#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
#include <stdio.h>
#include <tchar.h>// TODO: reference additional headers your program requires here
Friday, January 26, 2018 1:52 PM -
// spectrumratio.cpp : Defines the entry point for the console application.
//#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>#define NMAX 1000
#define M_PI 3.1415926535897932384626433832795typedef struct {
double re[NMAX];
double im[NMAX];
} CData;typedef struct {
double data[NMAX];
} Array;
CData cdft(CData xn ,int N )
{
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 ((double)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((double)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 _tmain(int argc, _TCHAR* argv[])
//int main ()
{
double delta =1e-15;
double xin[NMAX];
double yout[NMAX];
CData Xf1,Xf2,Xf3;
double Kmax=0, Kmin=0, Kmid=0;
int 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' ");
std::cin>>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 ");
char key1;
std::cin>>key1;
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' ");
std::cin>>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 ");
std::cin>>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 ");
char key2;
std::cin>>key2;return 0;
}
// stdafx.cpp : source file that includes just the standard includes
// spectrumratio.pch will be the pre-compiled header
// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"
// TODO: reference any additional headers you need in STDAFX.H
// and not in this file// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//#pragma once
#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
#include <stdio.h>
#include <tchar.h>// TODO: reference additional headers your program requires here
Friday, January 26, 2018 2:04 PM -
make.bat
set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"
javac THDObtainer.java
jar cvfe THDObtainer.jar THDObtainer *.class
jar uf THDObtainer.jar THDObtainer.java
pause 0set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"
rem java THDObtainer
java -jar THDObtainer.jarpause 0
import java.util.Scanner;
import java.io.*;
class THDParser
{
static int NMAX = 1024;public static double[] Getfassoc(double fs , int Numb )
{
double[] fvect =new double[Numb];
for(int k = 0;k<Numb;k++)
{
fvect[k] = fs * k / Numb;
}
return fvect;
}
public static double[] Gettassoc(double fs , int Numb )
{
double[] tvect=new double[Numb];
for(int N = 0 ; N<Numb; N++)
{
tvect[N] = N/fs ;
}
return tvect;
}
public static double Arg(double x_re , double x_im )
{
double Arg=0;
if (x_re > 0) { Arg = Math.atan(x_im / x_re); }
if ((x_re == 0) && (x_im > 0)) { Arg = Math.PI / 2;}
if ((x_re == 0) && (x_im < 0)) { Arg = Math.PI / 2 ; }
if ((x_re < 0) && (x_im >= 0)) { Arg = Math.atan(x_im / x_re) + Math.PI; }
if ((x_re < 0) && (x_im < 0)) { Arg = Math.atan(x_im / x_re) - Math.PI; }
if ((x_re == 0) && (x_im == 0)) {
Arg = 0;
/* MsgBox (" Xre=0 and Xim=0, Invalid result"); */
}
return Arg;
}
/*
'36. - 4. + 9.6568542i - 4. + 4.i - 4. + 1.6568542i - 4. - 4. - 1.6568542i
' column 7 to 8
' - 4. - 4.i - 4. - 9.6568542i
*/public static double[][] GetIQfromR_fi( int Num ,double [][] Xk, int mode )
{
//long int i ;
double arg1 ;
double[][] Xk1=new double [Num][2];
arg1 = 1;
if( mode == 1) { arg1 = Math.PI / 180;}
for(int i = 0 ;i< Num;i++)
{
Xk1[i][0] = Xk[i][0]*Math.cos(arg1*Xk[i][1]);
Xk1[i][1] = Xk[i][0]*Math.sin(arg1*Xk[i][1]);
}
return Xk;
}
public static double[][] CDFT(int numdft ,double [][] dft_xn )
{
int i, k, N;
double angle ;
double[][] dft_Xk=new double [numdft ][2];
for( k = 0 ;k<numdft ;k++)
{
dft_Xk[k][0] = 0;
dft_Xk[k][1] = 0;
for( N = 0;N<numdft;N++ )
{
angle = -2*Math.PI*k*N/numdft;
dft_Xk[k][0]+=dft_xn[N][0]*Math.cos(angle)-dft_xn[N][1]*Math.sin(angle);
dft_Xk[k][1]+=dft_xn[N][0]*Math.sin(angle)+dft_xn[N][1]*Math.cos(angle);
}
}
return dft_Xk;
}
public static double[][] ScaleDFT(int Num , double [][] Xk )
{
for(int i = 0 ;i<Num;i++)
{
Xk[i][0] = Xk[i][0]/Num;
Xk[i][1] = Xk[i][1]/Num;
}
return Xk;
}public static double[][] ScaleIQidft(int Num , double[][] ScXk)
{
for (int i = 0;i<Num;i++)
{
ScXk[i][0] = ScXk[i][0]*Num;
ScXk[i][1] = ScXk[i][1]*Num;
}return ScXk;
}public static double[][] CIDFT(int numidft , double[][] idft_Xk)
{
double angle;
double[][] idft_xn=new double[numidft][2] ;
for(int n = 0 ; n<numidft; n++)
{
idft_xn[n][0] = 0;
idft_xn[n][1] = 0;
for (int k=0; k<numidft; k++)
{
angle =2*Math.PI*k*n/numidft;
idft_xn[n][0]+=idft_Xk[k][0]*Math.cos(angle)-idft_Xk[k][1]*Math.sin(angle);
idft_xn[n][1]+=idft_Xk[k][0]*Math.sin(angle)+idft_Xk[k][1]*Math.cos(angle);
}
idft_xn[n][0]= idft_xn[n][0]/numidft;
idft_xn[n][1]= idft_xn[n][1]/numidft;
}
return idft_xn;
}public static double[][] Hilbert(double[] xre , int numb )
{
/* MATLAB algorithm */
int k, i, no2;
double[] h= new double [NMAX+1];
double[][] dft_xn =new double [NMAX][2];
double[][] dft_Xk =new double [NMAX][2];
double[][] idft_Xk =new double [NMAX][2];
double[][] idft_xn =new double [NMAX][2];if (numb == 0)
{
idft_xn[0][0]=0;
idft_xn[0][1]=0;
/*MsgBox (" Hilbert() : Parsing error ")*/
}
no2 = (int) Math.floor(numb / 2); //
for(i = 0;i<numb; i++)
{
dft_xn[i][0]=xre[i];
dft_xn[i][1]=0;
}
dft_Xk=CDFT(numb, dft_xn );
i=0;
for (i=0;i <= numb ;i++){ h[i]=0; }/*if N is even*/
if( (numb> 0)&&((2*no2)==numb) )
{
h[0] = 1; // 0
h[no2] = 1;// f Niquist
for( k = 1; k<no2;k++)
{
h[k]= 2;
}
}/* if N is odd */
else {
if(numb > 0)
{h[0] = 1;
for( k = 1; k<=((0.5*(numb + 1))-1); k++)
{
h[k] = 2;
}
}
}
for( i=0;i<numb;i++)
{
idft_Xk[i][0]=dft_Xk[i][0]*h[i];
idft_Xk[i][1]=dft_Xk[i][1]*h[i];
}
idft_xn=CIDFT(numb, idft_Xk);
return idft_xn;
}
public static double[][] GetAnalyticSignal(double [] xnre , int Numb )
{
// 'I=x(t)
// 'Q=jHilb(t)
// 'Hilb=Hilbert()
double[][] Hilb= new double[Numb][2];
double[][] x_iq= new double[Numb][2];Hilb= Hilbert(xnre, Numb );
for(int i = 0; i<Numb;i++)
{
x_iq[i][0] = xnre[i] ;
x_iq[i][1] = Hilb[i][1];
}
return x_iq;
}public static double[] getenvelope(double [] xnre , int Numb)
{
double[][] Hilb= new double[Numb][2];
double[] env= new double[Numb] ;
Hilb= Hilbert(xnre, Numb );for(int i = 0 ; i<Numb; i++)
{
env[i]=Math.sqrt((Hilb[i][0]*Hilb[i][0])+(Hilb[i][1]*Hilb[i][1]));
}
return env;}
public static double[][] GetSpectrumIQ(double[] xn , int Num )
{
double[][] x_iq= new double[Num ][2];
double[][] Xf_iq= new double[Num ][2];
x_iq=GetAnalyticSignal(xn , Num );
Xf_iq=CDFT(Num, x_iq);
Xf_iq=ScaleDFT(Num,Xf_iq);
return Xf_iq;
}
public static double[][] GetR_fi(int Numb , double[][] Xf_iq , int mode )
{double[][] Xf_rfi=new double [Numb][2];
for (int i = 0 ;i<Numb;i++)
{
Xf_rfi[i][0] =Math.sqrt((Xf_iq[i][0]*Xf_iq[i][0])+(Xf_iq[i][1]*Xf_iq[i][1]) );
Xf_rfi[i][1] = Arg(Xf_iq[i][0], Xf_iq[i][1]);
if( mode == 1) { Xf_rfi[i][1] = Xf_rfi[i][1]*180/Math.PI; }
}return Xf_rfi;
}
public static double[][] GetSpectrum(double[] xn , int Num , int mode )
{
double[][] x_iq= new double[Num][2];
double[][] Xf_iq= new double[Num][2];
double[][] Xf_rfi=new double [Num][2];x_iq=GetAnalyticSignal(xn, Num );
Xf_iq=CDFT(Num,x_iq);
Xf_iq=ScaleDFT(Num, Xf_iq);
Xf_rfi=GetR_fi(Num, Xf_iq, mode);
return Xf_rfi;
}
public static double[][] GetSpectrumfromIQ(double[][] x_iq,int Num, int mode)
{
double[][] Xf_iq= new double[Num][2];
double[][] Xf_rfi=new double [Num][2];Xf_iq=CDFT(Num,x_iq);
Xf_iq=ScaleDFT(Num, Xf_iq);
Xf_rfi=GetR_fi(Num,Xf_rfi, mode);
return Xf_rfi;
}
public static double[][] SynthSignalfromIQ(int Num , double[][] Xk_iq)
{
double[][] xn_iq= new double[Num][2];
xn_iq= CIDFT(Num, Xk_iq);
xn_iq= ScaleIQidft(Num, xn_iq);
return xn_iq;
}
public static double[][] SynthSignalfromRFi( double[][] Xk_rfi, int Num, int mode)
{
double[][] Xk_iq= new double [Num][2];
Xk_iq=GetIQfromR_fi(Num, Xk_rfi, mode );
return SynthSignalfromIQ(Num, Xk_iq);
}public static void printcomplex(double[][] xn , int num, int col1,int col2)
{
for (int i=0 ;i<num;i++)
{
System.out.print("\n i= ");
System.out.print(i);
System.out.print(" ");
System.out.print(String.format("%.08f",xn[i][0]));
System.out.print(" ");
System.out.print(String.format("%.08f",xn[i][1]));
}
System.out.print("\n ");
return;
}public static void printarray(double[] xn , int Num , int col1 )
{
for (int i = 0; i<Num;i++)
{
System.out.print("\n ");
System.out.print(String.format("%.08f",xn[i]));
}
System.out.print("\n ");
return;
}public static double deltat(double[] argF , double f0 , double fs , int N )
{
double[] fv =new double[N+1];
int N0;
double deltat;
fv=Getfassoc(fs, N);
N0 = 0;
for(int i = 0;i<N-1;i++)
{
if ((fv[i]>= f0)&&(fv[i+1]<=f0)) { N0= i; }
}
deltat =-(argF[N0] -argF[0])/(2*Math.PI*fv[N0]-2*Math.PI*fv[0]);
return deltat;
}
public static void writetofile(double[][] xnc , String FilePath , int N ) throws IOException
{
File outFile1 = new File (FilePath);
FileWriter fWriter1 = new FileWriter (outFile1);
PrintWriter pWriter1 = new PrintWriter (fWriter1);for(int i = 0 ;i< N;i++)
{
pWriter1.println("\n ");
pWriter1.println(" ");
pWriter1.print(String.format("%d",i) );
pWriter1.print(" ");
pWriter1.print(String.format("%e",xnc[i][0] ));
pWriter1.print(" ");
pWriter1.print(String.format("%e",xnc[i][1] ));
pWriter1.print(" ");
}
pWriter1.close();
return;
}
public static void TestHilbert( )
{int i ;
int N = 8;
double fs ;
double[] fv=new double[N ];
double[] tv=new double[N ];
double[][] Xf=new double[N ][2];double[][] Yf=new double[N ][2];
double[][] Hilb=new double[N][2];
double[][] Idft=new double[N ][2];
double[][] SpX=new double[N][2];
double[] env=new double[N];
double[][] Xn=new double[N][2];
double[] Xn1=new double[N];
double[] x_re=new double[N];
double[][] Yn=new double[N][2];
Xf[0][0] = 1;
Xf[1][0] = 2;
Xf[2][0] = 3;
Xf[3][0] = 4;
Xf[4][0] = 5;
Xf[5][0] = 6;
Xf[6][0] = 7;
Xf[7][0] = 8;Xf[0][1] = 0;
Xf[1][1] = 0;
Xf[2][1] = 0;
Xf[3][1] = 0;
Xf[4][1] = 0;
Xf[5][1] = 0;
Xf[6][1] = 0;
Xf[7][1] = 0;
fs = 1000;printcomplex(Xf,N,1,2);
fv=Getfassoc(fs, N);
tv=Gettassoc(fs, N);
System.out.print("\n f= ");
printarray(fv, N, 0);
System.out.print("\n t= ");
printarray(tv, N, 0);
Yn=CDFT(N,Xf);
System.out.print("\n CDFT(XF)= ");
printcomplex(Yn, N, 5, 6);
Idft=CIDFT(N, Yn );
System.out.print("\n CIDFT(Yn)= ");
printcomplex(Idft, N, 7, 8);
for(i=0;i<N;i++)
{
Xn1[i]=Idft[i][0];
}Hilb=Hilbert(Xn1, N );
System.out.print("\n Hilbert(Xn1)= ");
printcomplex(Hilb, N, 9, 10);
SpX=GetSpectrum(Xn1, N , 1);
System.out.print("\n Get spectrum= ");
printcomplex(SpX, N, 11, 12);
System.out.print("\n Get envelope (Xn1) ");
env=getenvelope(Xn1, N);
printarray(env, N, 13);
System.out.print("\nPoke new aray ");// 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
N = 16;
double [][] Xn2=new double[N][2] ;
double [][] Xf2=new double[N][2] ;
for (i= 0;i<N; i++)
{
Xf2[i][0] = 0;
Xf2[i][1] = 0;
}
Xf2[1][0] = 1;
Xf2[1][1] = -45;
Xf2[5][0] = 0;
Xf2[5][1] = 0;
Xn2 =SynthSignalfromRFi(Xf2, N, 1);
System.out.print("\n Get Xn from Xf");
printcomplex(Xn2, N, 14, 15);
fv=Getfassoc(fs, N);
System.out.print("\n f[]=");
printarray(fv, N, 16);
System.out.print("\n SpX[]=");
SpX=GetSpectrumfromIQ(Xn2, N, 1);
printcomplex(SpX, N, 17, 18);
double[] Xn4=new double[N];
for (i = 0 ; i<N; i++)
{
Xn4[i] = 0.5*Math.cos(2*Math.PI*250*i/fs+0);
Xn4[i] +=1.5*Math.cos(2*Math.PI*125*i/fs+Math.PI/2);
Xn4[i] +=2*Math.cos(2*Math.PI*375 *i/fs-Math.PI/2);
}double [][] Xn3=new double[N][2] ;
double [][] Xf3=new double[N][2] ;
double [][] Yn3=new double[N][2] ;
System.out.print("\n Spectrum (Xn4)=");
Xf3=GetSpectrum(Xn4, N , 1);
printcomplex(Xf3, N, 19, 20);
System.out.print("\n Spectrum IQ(Xn4)=");
Yn3=GetSpectrumIQ(Xn4, N );
printcomplex(Yn3, N, 21, 22);
return;}
public int GetNumF1(double f1 ,int Numb ,double fs )
{
return (int) (Math.floor(f1*Numb/fs));
}public static double GetTHD_F(double[] xn , int Numb, double fs, double f1)
{int numbf1;
double[][] Xf =new double[Numb][2];
double Urms1h;
double Uhharm;
int i ;//
//f 0 5 10 15 20 25 30 35 Hz N=8
//i 0 1 2 3 4 5 6 7
// f1 2f1 3f1
// numb1=2
//
//f 0 5 10 15 20 25 30 35 40 45 50 55 60 Hz N=13 imax=12
//i 0 1 2 3 4 5 6 7 8 9 10 11 12
// f1 2f1 3f1 4f1
// numb1=3
//numbf1=(int)Math.floor(f1*Numb/fs);// 2 , get associated vector of frequency
Xf=GetSpectrum(xn, Numb, 0);
Urms1h= Xf[numbf1][0]/Math.sqrt(2);
i = numbf1*2;
Uhharm=0;
while (i<Numb)
{
if (i<Numb) { Uhharm+=(Xf[i][0]*Xf[i][0])/2 ; }
i=i+numbf1;
}Uhharm =Math.sqrt(Uhharm) ;
return Uhharm*100/(Urms1h);
}
public static double GetTHD_R(double[] xn , int Numb, double fs, double f1)
{
int numbf1;
double[][] Xf =new double[Numb][2];
double Urmssigah;
double Uhharm;
int i ;
numbf1=(int)Math.floor(f1*Numb/fs); // get associated vector of frequency, ORIGIN=0
Xf=GetSpectrum(xn, Numb, 0);
/* parser of Urms of all harmonics */// THD-R=Urmshharm*100/urmsig =sqrt( sum(Xrms[i] ^2))/Vrmssig;
Urmssigah=0;
i=numbf1;while (i<Numb)
{
if (i<Numb) { Urmssigah+=(Xf[i][0]*Xf[i][0])/2; }
i=i+1;
}Urmssigah=Math.sqrt (Urmssigah );
Uhharm=0;
/* parser of Urms of higher harmonics */i=numbf1*2; // ' number of bin of 2-nd harmonic
while (i<Numb)
{
if (i<Numb) { Uhharm+=(Xf[i][0]*Xf[i][0])/2; }
i=i+numbf1;
}Uhharm=Math.sqrt (Uhharm) ;
return (Uhharm*100/Urmssigah);
}public static double GetRMS(double[] xn , int Numb )
{
int i ;
double Urms;
Urms = 0;
for (i=0;i<Numb;i++)
{
Urms+=xn[i]*xn[i];
}
Urms=Math.sqrt(Urms/Numb);return Urms;
}
public static double GetTHD_N(double[] xn , int Numb, double fs , double f1 )
{
int numbf1;
double[][] Xf=new double [Numb][2];
double Urmssig;
double Uhharmn;
int i ;
/*
Urmssig = 0;
for(i=0 ;i<Numb;i++)
{
Urmssig+=xn[i]*xn[i];
}
Urmssig=Math.sqrt(Urmssig/Numb);*/
numbf1=(int)Math.floor (f1*Numb/fs);
Xf=GetSpectrum(xn, Numb, 0);i=0;
Urmssig=0;
while (i<Numb)
{
{ Urmssig+=Xf[i][0]*Xf[i][0]/2; }
i=i+1;
}Urmssig= Math.sqrt(Urmssig);
Uhharmn=0;
i=0;while (i<Numb)
{
if (i!= numbf1) { Uhharmn+=Xf[i][0]*Xf[i][0]/2; }
i=i+1;
}Uhharmn= Math.sqrt(Uhharmn);
return (Uhharmn*100/Urmssig);
}public static void PrintValueDouble(String str, double s,String str1)
{
System.out.print(str);
System.out.print(String.format("%.08f",s));
System.out.print(str1);
return;
}public static void PrintValueInt(String str, int s,String str1)
{
System.out.print(str);
System.out.print(String.format("%d",s));
System.out.print(str1);
return;
}public static double GetValueDouble(String s,Scanner scan1)
{
System.out.print(s);double result;
result=(double) scan1.nextDouble();
return result;
}
public static int GetValueInt(String s,Scanner scan1)
{
System.out.print(s);
int result=(int) scan1.nextInt();
return result;
}public static void TestTHDN( )
{
int N;
Scanner sc2 = new Scanner (System.in);
PrintValueInt("\n N=",16," ");
N =16;// GetValueInt("\n Input N ",sc2);double [][] Xf =new double[N][2];
double [][] X_iq=new double[N ][2];
double [][] Xn =new double[N ][2] ;
double [] Xn1=new double[N ] ;
double [] fvect =new double[N ] ;
double fs;
double f1;
int i;fs =44100;// GetValueDouble("\n Input fs,Hz ",sc2);
f1 =11025;//GetValueDouble("\n Input f [1-st harm] ,Hz ",sc2);
fvect= Getfassoc(fs, N);for(i = 0;i<N;i++)
{
//PrintValueDouble("\n f=",fvect[i]);
Xf[i][0] =0;//GetValueDouble("\n input Um[f] " ,sc2);
Xf[i][1] =0;//GetValueDouble("\n input fi[f] ",sc2);}
PrintValueDouble("\n fs=",fs," ");
PrintValueDouble("\n f1=",f1," ");PrintValueDouble("\n f=",fvect[4]," Hz");
Xf[4][0]=1;
PrintValueDouble(" |X[f]|=",Xf[4][0]," ");PrintValueDouble("\n f=",fvect[8]," Hz");
Xf[8][0]=0.3;
PrintValueDouble(" |X[f]|=",Xf[8][0]," ");PrintValueDouble("\n f=",fvect[9]," Hz");
Xf[9][0]=0.1;
PrintValueDouble(" |X[f]|=",Xf[9][0]," ");for(i = 0;i<N;i++)
{
PrintValueDouble("\n f=",fvect[i]," ");
PrintValueDouble(" |X[f]|=",Xf[i][0], " ");
PrintValueDouble(" fi=",Xf[i][1]," ");}
Xn=SynthSignalfromRFi(Xf, N, 1);
printcomplex(Xn, N, 5, 6);
for (i=0;i<N;i++)
{
Xn1[i]=Xn[i][0];
}
double THD_N= GetTHD_N(Xn1, N, fs, f1);
PrintValueDouble("\n THD+N= ", THD_N," %");
double THD_R= GetTHD_R(Xn1, N, fs, f1);
PrintValueDouble("\n THD-R= ",THD_R," %");
double THD_F= GetTHD_F(Xn1, N, fs, f1 );
PrintValueDouble("\n THD-F= ",THD_F," %");
double RMS=GetRMS(Xn1, N);
PrintValueDouble("\n U RMS=",RMS," units of ADC");
sc2.close();
return ;
}}
public class THDObtainer
{
public static void main(String[] args) throws IOException
{
THDParser.TestHilbert();
THDParser.TestTHDN();
}}
Friday, January 26, 2018 2:08 PM -
spectrumstat - utility for obtaining max. levels of signals on frequencies of the spectrum using short-time Fourier transform and founder of the max. vaues
xinputns.txt
4,410000e+04 16 8
0 1,000000e+01
1 -7,071068e+00
2 3,061617e-15
3 7,071068e+00
4 -1,000000e+01
5 7,071068e+00
6 -2,694842e-14
7 -7,071068e+00
8 1,000000e+01
9 -7,071068e+00
10 3,061617e-15
11 7,071068e+00
12 -1,000000e+01
13 7,071068e+00
14 -2,694842e-14
15 -7,071068e+004,410000e+04 16
maxspectr.txt
0 3,020133e-15 0,000000e+00
1 1,330312e-07 5,512500e+03
2 6,447902e-15 1,102500e+04
3 1,000000e+01 1,653750e+04
4 4,642237e-15 2,205000e+04
5 5,151483e-15 2,756250e+04
6 2,334110e-15 3,307500e+04
7 7,772354e-15 3,858750e+04make.bat
set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"
javac Spectrumstat.java
jar cvfe Spectrumstat.jar Spectrumstat *.class
jar uf Spectrumstat.jar Spectrumstat.java xinputns.txt
pause 0run.bat
set PATH="C:\Program Files\Java\jdk1.7.0_17\bin"
java -jar Spectrumstat.jar
pause 0//javac Spectrum.java
//java Spectrumimport java.util.Scanner;
import java.io.*;
public class Spectrumstat
{public static double[][] cdft(double[][] xn ,int N )
{
int j;
double[][] Xk=new double[N][2];
int k,n ;
double angle;
for (k=0;k<N;k++)
{
Xk[k][0]=0;
Xk[k][1]=0;
for(n=0;n<N;n++)
{
angle=-2*Math.PI*k*n/N;
Xk[k][0]+= xn[n][0]*Math.cos(angle)-xn[n][1]*Math.sin(angle);
Xk[k][1]+= xn[n][0]*Math.sin(angle)+xn[n][1]*Math.cos(angle);
}
}return Xk;
}
public static double[][] ScaleDFT(double[][] xk,int num)
{
for (int i=0; i<num ; i++)
{
xk[i][0] =xk[i][0] /num;
xk[i][1] =xk[i][1] /num;
}
return xk;
}public static double[][] hilbert(double[] xr )
{
double angle;
double[] h=new double[xr.length+1] ;
double[] dftxn_re=new double[xr.length] ;
double[] dftxn_im=new double[xr.length] ;
double[] dftXk_re=new double[xr.length] ;
double[] dftXk_im=new double[xr.length] ;
double[] idftXK_re=new double[xr.length] ;
double[] idftXK_im=new double[xr.length] ;
double[][] idftXN=new double[xr.length][2] ;
if (xr.length==0)
{
System.out.println("\n Error");
idftXN[0][0] =0;
idftXN[0][1]=0;
return idftXN;
}
//xnreal=real(xr);
for (int i=0;i<xr.length;i++)
{
dftxn_re[i]=xr[i];
dftxn_im[i]=0;
}//fft
// dftXk=dft(dftxn,n);for (int k=0;k<xr.length;k++)
{
dftXk_re[k]=0;
dftXk_im[k]=0;
for(int n=0;n<xr.length;n++)
{
angle=-2*Math.PI*k*n/xr.length;
dftXk_re[k]+= dftxn_re[n]*Math.cos(angle)- dftxn_im[n]*Math.sin(angle);
dftXk_im[k]+= dftxn_re[n]*Math.sin(angle)+ dftxn_im[n]*Math.cos(angle);
}}
//X={Xk.re,Xk.im};
for (int i=0;i<=xr.length;i++)
{
h[i]=0;
}
//int no2=(int) Math.floor (xr.length/2);if((xr.length>0)&&((2* Math.floor(xr.length/2))==(xr.length))) // n is even
{h[0]=1;
h[xr.length/2]=1;
for(int k=1;k<xr.length/2;k++)
{
h[k]=2;
}
}
else //n is odd
{
if(xr.length>0)
{
h[0]=1;
for(int k=1;k<(xr.length+1)/2;k++)
{
h[k]=2;
}
}
}for(int i=0;i<xr.length;i++)
{
idftXK_re[i]= dftXk_re[i]*h[i];
idftXK_im[i]= dftXk_im[i]*h[i];
}
//idftXN= idft(idftXK, n);
//N=sizeof(xk)/sizof(double);
//ifft
for (int n=0;n<xr.length;n++)
{
idftXN[n][0] =0;
idftXN[n][1] =0;
for(int k=0;k<xr.length;k++)
{
angle=2*Math.PI*k*n/xr.length;
idftXN[n][0] +=idftXK_re[k]*Math.cos(angle)-idftXK_im[k]*Math.sin(angle);
idftXN[n][1] +=idftXK_re[k]*Math.sin(angle)+idftXK_im[k]*Math.cos(angle);
}
idftXN[n][0] =idftXN[n][0] /xr.length;
idftXN[n][1] =idftXN[n][1] /xr.length;}
return idftXN;
}public static double[][] getanalyticsignal(double [] xn )
{
int i;
double[][] Hilb =new double[xn.length][2];
double[][] x_n =new double[xn.length][2];Hilb=hilbert(xn );
for (i=0; i<xn.length; i++)
{
x_n[i][0] =xn[i];
x_n[i][1] =Hilb[i][1] ;
}
return x_n;
}
public static double Arg(double re, double im)
{
return Math.atan2(im,re) ;
}public static double[][] get_mod_fi(double[][] x, int n )
{
double[][] res=new double[ n][2] ;
double fi;
for (int i=0;i<n;i++)
{
res[i][0] =Math.sqrt(x[i][0] *x[i][0] +x[i][1] *x[i][1] );
res[i][1] =Arg(x[i][0] ,x[i][1] );
{ res[i][1] =res[i][1] *180/Math.PI; }
}
return res;
}
public static double[][] GetSpectrum( double [] xn )
{
double[][] xn_iq=new double[xn.length][2];
double[][] Xf_IQ =new double[xn.length][2] ;
double[][] Xf_modfi =new double[xn.length][2];
xn_iq =getanalyticsignal( xn );
Xf_IQ =cdft(xn_iq, xn.length);
Xf_IQ =ScaleDFT(Xf_IQ,xn.length);
Xf_modfi =get_mod_fi(Xf_IQ, xn.length );
return Xf_modfi ;
}
public static double[][][] GetSpectrumarray( double [] xn , int num,int nseq)
{
//int np=(int) num/nseq;
int npages=(int) Math.floor (num/nseq);double [] xnseq =new double [nseq];
double [][] Xf = new double [nseq][2];
double [][][] Xfarray= new double [npages][nseq][2];
for (int i=0;i<npages; i++) //0,1,32
{
for (int j=0;j<nseq;j++)
{
xnseq[j]=xn[i+j]; //if nseq=16, xnseq[j]=xn[k],k= 0...15,16...31
// if ((i+j)<num) &&( (nseq*npages)>(num/nseq) ) xnseq[j]=0;
}Xf=GetSpectrum(xnseq );
for (int j=0;j<nseq; j++) //0,1,32
{
Xfarray[i][j][0]= Xf[j][0];
Xfarray[i][j][1]=Xf[j][1];
}}
return Xfarray;
}public static double get_max(double a,double b )
{
if (a>b) {return a;}return b;
}
public static double[] GetPeakSpectrumComponents(double [] xn , int num, int nseq)
{
//int np=(int)num/nseq;
int npages=(int) Math.floor (num/nseq);
double [][][] Xfarray =new double [npages][nseq][2];
double [] Xfpeak=new double[nseq];
Xfarray=GetSpectrumarray( xn , num, nseq);for (int j=0;j<nseq;j++) // for each Xn [frequency](page)
{
Xfpeak[j]=0;
for (int i=0;i<npages;i++) //for pages of spectrum
{
Xfpeak[j]=get_max(Xfpeak[j], Xfarray[i][j][0] );
}}
return Xfpeak;
}
public static void main(String[] args) throws IOException
{Scanner sc = new Scanner(System.in);
int NMAX=32767;
int N=0,n,i, ns=0;
double[] xn =new double[NMAX] ;
double fs=0,T;
System.out.print ("\n For input x[n] from file input '1' ");
System.out.print ("\n For manual input x[n] input '0' ");
int key=sc.nextInt();
switch( key)
{case 1:
System.out.print ("\n Open xinput.txt \n");
String fileNameDefined = "xinputns.txt";
File file1 = new File(fileNameDefined);
try{
// -read from filePooped with Scanner class
Scanner inputStream = new Scanner(file1);
fs = (double) inputStream.nextDouble();
N =inputStream.nextInt();
ns =inputStream.nextInt();
System.out.print ("\n ns=");
System.out.print (ns);
for (n=0 ; n<N ;n++)
{
inputStream.nextInt();
xn[n]= inputStream.nextDouble();}
inputStream.close();
} catch (FileNotFoundException e) {e.printStackTrace();
}
break;case 0:
System.out.print ("\n Input fs input,Hz ");
fs = (double) sc.nextDouble();
System.out.print ("\n Input number of samples of xin[n] ");
N =sc.nextInt();
System.out.print ("\n Input number of samples in CDFT sequence ");
ns=sc.nextInt();
xn =new double[N] ;
for (n=0 ; n<N ;n++)
{
T=n*1/(N *fs );
System.out.print ("\n n=");
System.out.print (n);
System.out.print (" t=");
System.out.print (T);
System.out.print (" s; xn[n]= " );
xn[n]= sc.nextDouble();
}
break;
default:
break;
}
System.out.print ("\n fs =");
System.out.print (fs);
System.out.print (" Hz");
System.out.print ("\n Nsamples1=");
System.out.print (N);for (n=0 ; n<N ;n++)
{
T=n*1/(N *fs );
System.out.print ("\n n=");
System.out.print (n);
System.out.print (" x[n]=");
System.out.print (xn[n]);
System.out.print (" t= ");
System.out.print (T);
System.out.print (" s");
}
System.out.println ("\n Input 0 for continue ");
sc.nextLine();
System.out.print ("\n OK");System.out.println("Array of Xn is: ");
for ( i=0;i<N;i++)
{
System.out.print("\n");
System.out.print("N= ");
System.out.print(i);
System.out.print (" ; xn[n]= ");
System.out.print(xn[i]);
}double Xfmax[] =new double[ns] ;
System.out.println("\n Parsing spectrum ");
Xfmax=GetPeakSpectrumComponents(xn, N, ns );
//Xf=GetSpectrum(xn, N);
System.out.println("\n Array of Xfmax is: ");
for (i=0;i<ns ;i++)
{
System.out.print("\n");
System.out.print("ns= ");
System.out.print(i);System.out.print (" ; max |Xf[k]|= ");
System.out.print(Xfmax[i] );
System.out.print (" ; f=");
System.out.print ( fs*i/ns );
System.out.print(" Hz");}
File outFile = new File ("maxspectr.txt");
FileWriter fWriter = new FileWriter (outFile);
PrintWriter pWriter = new PrintWriter (fWriter);
pWriter.println("\n");
pWriter.print(String.format("%e",fs ));
pWriter.print(" ");
pWriter.print(String.format("%d",N ) );
pWriter.println("\n");
for (i=0; i<ns ;i++)
{
pWriter.println("\n ");
pWriter.print(String.format("%d",i) );
pWriter.print(" ");
pWriter.print(String.format("%e",Xfmax[i] ));
pWriter.print(" ");
pWriter.print (String.format("%e", fs*i/ns) );
pWriter.print(" ");
}
pWriter.close();System.out.println("\n End of sequence ");
sc.close();
}}
Friday, January 26, 2018 3:01 PM -
You may convert Java programs into VC++ console applications (if the problem of work with big 2-dimension arrays is solved ).Friday, January 26, 2018 3:04 PM
-
Program for obtaining envelope of sinal sequence using Hilbert transform .
// getenv.cpp : Defines the entry point for the console application.
//#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <iostream>
#define NMAX 1000
#define M_PI 3.1415926535897932384626433832795typedef struct {
double re[NMAX];
double im[NMAX];
} CData;typedef struct {
double data[NMAX];
} Array;
CData cdft(CData xn ,int N )
{
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 ((double)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((double)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 _tmain(int argc, _TCHAR* argv[])
//int main ()
{
double delta =1e-15;
double xin[NMAX];
double yout[NMAX];
Array Xenv1 ;
double Kmax=0, Kmin=0, Kmid=0;
int 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' ");
std::cin>>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 ");
std::cin>>key;
printf("\n OK. ");
/*****************************/label1:
printf("\n Obtaining envelope of x[n]");
Xenv1=getenvelope(xin, N1);printf("\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++)
{
printf ("\n %d |Hilbert(x[n])|= %le ", i, Xenv1.data[i] );
fprintf (fp2,"\n %d %le ", i, Xenv1.data[i]);}
fclose(fp2);printf ("\n input '0' for continue ");
std::cin>>key;Kmax=getmax(Xenv1,N1);
Kmin=getmin(Xenv1,N1);
Kmid=getmid(Xenv1,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 ");
std::cin>>key;return 0;
}
// stdafx.cpp : source file that includes just the standard includes
// getenv.pch will be the pre-compiled header
// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"
// TODO: reference any additional headers you need in STDAFX.H
// and not in this file// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//#pragma once
#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
#include <stdio.h>
#include <tchar.h>// TODO: reference additional headers your program requires here
Friday, January 26, 2018 4:19 PM -
Program for obtaining of the Bessel's functions for obtainer of the PWM and FM spectrums of signals.
// bessel.cpp : Defines the entry point for the console application.
//#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include "bessel.h"
int main(){
int n;
double x;
double result1, result2,result3,result4;std::cout<<"Input order n : "<<"\n";
std::cin>>n;
std::cout<<"Input x : "<<"\n";
std::cin>>x;
result1=jn(n,x);
std::cout<<"jn(n,x)= "<<result1<<"\n";
result2=bessj( n, x );
std::cout<<"jn(n,x)= "<<result2<<"\n";
result3=yn(n,x);
std::cout<<"yn(n,x)= "<<result3<<"\n";
result4=bessy( n, x );
std::cout<<"yn(n,x)= "<<result4<<"\n";std::cout<<"Input '0' for exit";
char key;
std::cin>>key;
return 0;
}//bessel.h (add into project using Add Item-> .h file,then put the code into file)
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define ACC 40.0
#define BIGNO 1.0e10
#define BIGNI 1.0e-10/*
http://www.astro.rug.nl/~gipsy/sub/bessel.c
*/static double bessj0( double x )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate Bessel function of first kind and order */
/* 0 at input x */
/*------------------------------------------------------------*/
{
double ax,z;
double xx,y,ans,ans1,ans2;if ((ax=fabs(x)) < 8.0) {
y=x*x;
ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
+y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
+y*(59272.64853+y*(267.8532712+y*1.0))));
ans=ans1/ans2;
} else {
z=8.0/ax;
y=z*z;
xx=ax-0.785398164;
ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
+y*(-0.2073370639e-5+y*0.2093887211e-6)));
ans2 = -0.1562499995e-1+y*(0.1430488765e-3
+y*(-0.6911147651e-5+y*(0.7621095161e-6
-y*0.934935152e-7)));
ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
}
return ans;
}static double bessj1( double x )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate Bessel function of first kind and order */
/* 1 at input x */
/*------------------------------------------------------------*/
{
double ax,z;
double xx,y,ans,ans1,ans2;if ((ax=fabs(x)) < 8.0) {
y=x*x;
ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
+y*(99447.43394+y*(376.9991397+y*1.0))));
ans=ans1/ans2;
} else {
z=8.0/ax;
y=z*z;
xx=ax-2.356194491;
ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+y*(0.2457520174e-5+y*(-0.240337019e-6))));
ans2=0.04687499995+y*(-0.2002690873e-3
+y*(0.8449199096e-5+y*(-0.88228987e-6
+y*0.105787412e-6)));
ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
if (x < 0.0) ans = -ans;
}
return ans;
}double bessj( int n, double x )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate Bessel function of first kind and order */
/* n at input x */
/* The function can also be called for n = 0 and n = 1. */
/*------------------------------------------------------------*/
{
int j, jsum, m;
double ax, bj, bjm, bjp, sum, tox, ans;
if (n < 0)
{
double dblank;
printf("Error");
// setdblank_c( &dblank );
return( 0 );
}
ax=fabs(x);
if (n == 0)
return( bessj0(ax) );
if (n == 1)
return( bessj1(ax) );if (ax == 0.0)
return 0.0;
else if (ax > (double) n) {
tox=2.0/ax;
bjm=bessj0(ax);
bj=bessj1(ax);
for (j=1;j<n;j++) {
bjp=j*tox*bj-bjm;
bjm=bj;
bj=bjp;
}
ans=bj;
} else {
tox=2.0/ax;
m=2*((n+(int) sqrt(ACC*n))/2);
jsum=0;
bjp=ans=sum=0.0;
bj=1.0;
for (j=m;j>0;j--) {
bjm=j*tox*bj-bjp;
bjp=bj;
bj=bjm;
if (fabs(bj) > BIGNO) {
bj *= BIGNI;
bjp *= BIGNI;
ans *= BIGNI;
sum *= BIGNI;
}
if (jsum) sum += bj;
jsum=!jsum;
if (j == n) ans=bjp;
}
sum=2.0*sum-bj;
ans /= sum;
}
return x < 0.0 && n%2 == 1 ? -ans : ans;
}
static double bessy0( double x )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate Bessel function of second kind and order */
/* 0 at input x. */
/*------------------------------------------------------------*/
{
double z;
double xx,y,ans,ans1,ans2;if (x < 8.0) {
y=x*x;
ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
+y*(10879881.29+y*(-86327.92757+y*228.4622733))));
ans2=40076544269.0+y*(745249964.8+y*(7189466.438
+y*(47447.26470+y*(226.1030244+y*1.0))));
ans=(ans1/ans2)+0.636619772*bessj0(x)*log(x);
} else {
z=8.0/x;
y=z*z;
xx=x-0.785398164;
ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
+y*(-0.2073370639e-5+y*0.2093887211e-6)));
ans2 = -0.1562499995e-1+y*(0.1430488765e-3
+y*(-0.6911147651e-5+y*(0.7621095161e-6
+y*(-0.934945152e-7))));
ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
}
return ans;
}static double bessy1( double x )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate Bessel function of second kind and order */
/* 1 at input x. */
/*------------------------------------------------------------*/
{
double z;
double xx,y,ans,ans1,ans2;if (x < 8.0) {
y=x*x;
ans1=x*(-0.4900604943e13+y*(0.1275274390e13
+y*(-0.5153438139e11+y*(0.7349264551e9
+y*(-0.4237922726e7+y*0.8511937935e4)))));
ans2=0.2499580570e14+y*(0.4244419664e12
+y*(0.3733650367e10+y*(0.2245904002e8
+y*(0.1020426050e6+y*(0.3549632885e3+y)))));
ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x);
} else {
z=8.0/x;
y=z*z;
xx=x-2.356194491;
ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+y*(0.2457520174e-5+y*(-0.240337019e-6))));
ans2=0.04687499995+y*(-0.2002690873e-3
+y*(0.8449199096e-5+y*(-0.88228987e-6
+y*0.105787412e-6)));
ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
}
return ans;
}
double bessy( int n, double x )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate Bessel function of second kind and order */
/* n for input x. (n >= 0) */
/* Note that for x == 0 the functions bessy and bessk are not */
/* defined and a blank is returned. */
/*------------------------------------------------------------*/
{
int j;
double by,bym,byp,tox;
if (n < 0 || x == 0.0)
{
double dblank;
printf("Error"); // setdblank_c( &dblank );
return(0 );
}
if (n == 0)
return( bessy0(x) );
if (n == 1)
return( bessy1(x) );tox=2.0/x;
by=bessy1(x);
bym=bessy0(x);
for (j=1;j<n;j++) {
byp=j*tox*by-bym;
bym=by;
by=byp;
}
return by;
}static double bessi0( double x )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate modified Bessel function In(x) and n=0. */
/*------------------------------------------------------------*/
{
double ax,ans;
double y;
if ((ax=fabs(x)) < 3.75) {
y=x/3.75,y=y*y;
ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
+y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
} else {
y=3.75/ax;
ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
+y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
+y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
+y*0.392377e-2))))))));
}
return ans;
}
static double bessi1( double x)
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate modified Bessel function In(x) and n=1. */
/*------------------------------------------------------------*/
{
double ax,ans;
double y;
if ((ax=fabs(x)) < 3.75) {
y=x/3.75,y=y*y;
ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
+y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
} else {
y=3.75/ax;
ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
-y*0.420059e-2));
ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
+y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
ans *= (exp(ax)/sqrt(ax));
}
return x < 0.0 ? -ans : ans;
}
double bessi( int n, double x)
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate modified Bessel function In(x) for n >= 0*/
/*------------------------------------------------------------*/
{
int j;
double bi,bim,bip,tox,ans;
if (n < 0)
{
double dblank;
printf("Error");// setdblank_c( &dblank );
return( 0);
}
if (n == 0)
return( bessi0(x) );
if (n == 1)
return( bessi1(x) );
if (x == 0.0)
return 0.0;
else {
tox=2.0/fabs(x);
bip=ans=0.0;
bi=1.0;
for (j=2*(n+(int) sqrt(ACC*n));j>0;j--) {
bim=bip+j*tox*bi;
bip=bi;
bi=bim;
if (fabs(bi) > BIGNO) {
ans *= BIGNI;
bi *= BIGNI;
bip *= BIGNI;
}
if (j == n) ans=bip;
}
ans *= bessi0(x)/bi;
return x < 0.0 && n%2 == 1 ? -ans : ans;
}
}
static double bessk0( double x )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate modified Bessel function Kn(x) and n=0. */
/*------------------------------------------------------------*/
{
double y,ans;if (x <= 2.0) {
y=x*x/4.0;
ans=(-log(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420
+y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
+y*(0.10750e-3+y*0.74e-5))))));
} else {
y=2.0/x;
ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
+y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
+y*(-0.251540e-2+y*0.53208e-3))))));
}
return ans;
}
static double bessk1( double x )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate modified Bessel function Kn(x) and n=1. */
/*------------------------------------------------------------*/
{
double y,ans;if (x <= 2.0) {
y=x*x/4.0;
ans=(log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144
+y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1
+y*(-0.110404e-2+y*(-0.4686e-4)))))));
} else {
y=2.0/x;
ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619
+y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2
+y*(0.325614e-2+y*(-0.68245e-3)))))));
}
return ans;
}double bessk( int n, double x )
/*------------------------------------------------------------*/
/* PURPOSE: Evaluate modified Bessel function Kn(x) and n >= 0*/
/* Note that for x == 0 the functions bessy and bessk are not */
/* defined and a blank is returned. */
/*------------------------------------------------------------*/
{
int j;
double bk,bkm,bkp,tox;
if (n < 0 || x == 0.0)
{
double dblank;
printf("Error");// setdblank_c( &dblank );
return( 0 );
}
if (n == 0)
return( bessk0(x) );
if (n == 1)
return( bessk1(x) );tox=2.0/x;
bkm=bessk0(x);
bk=bessk1(x);
for (j=1;j<n;j++) {
bkp=bkm+j*tox*bk;
bkm=bk;
bk=bkp;
}
return bk;
}
#undef ACC
#undef BIGNO
#undef BIGNI// stdafx.cpp : source file that includes just the standard includes
// bessel.pch will be the pre-compiled header
// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"
// TODO: reference any additional headers you need in STDAFX.H
// and not in this file// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//#pragma once
#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
#include <stdio.h>
#include <tchar.h>
//#include "bessel.h"
// TODO: reference additional headers your program requires here- Edited by USERPC01 Friday, January 26, 2018 6:57 PM bug fix
Friday, January 26, 2018 6:54 PM -
project for obtaining IR coefficients of FIR from AFR data with #include <vector>, #include <complex>, #include <cmath>
//afrtoht.cpp
// afrtoht.cpp : main project file.
#include "stdafx.h"
#include <stdio.h>
/*
#include <tchar.h>
#include <iostream>
#include <stdlib.h>
#include <cmath>
*/
#include <iostream>
#include <fstream>
#include <complex>
#include <vector>
#include <set>
#include <cmath>
#include <limits>using namespace std;
#define M_PI 3.1415926535897932384626433832795
//int _tmain(int argc, _TCHAR* argv[])
//{
// return 0;
//}//afrtohtcoefs.cpp
/* Define complex double data type */
typedef complex<double> dcomp;//std::vector<std::vector<dcomp> >
typedef std::vector<dcomp> cplxvector ;#define NMAX 512
//FILE *fp, *fp1;
int k,n;
int w,t,N;
double h[NMAX];
double h_im[NMAX];
double Re_H[NMAX];
double Im_H[NMAX];
double Re_x[NMAX];
double Im_x[NMAX];
double abs_H[NMAX];
double arg_H[NMAX];
double f[NMAX];
float angle;
double K, fs,T,finp,Tmax;
char str;
char key;/*
float heavi( double t)
{
if (t<0 ) return 0;
if (t>=0) return 1;
}float heavi1( int n)
{
if (n<0 ) return 0;
if (n>=0) return 1;
}*/double getarg(float Re_H,float Im_H)
{return atan2(Im_H,Re_H) * 180 / M_PI;
}cplxvector CIDFT(cplxvector idftXk, int N)
{
double angle;
cplxvector idftxn;
dcomp data;//N=sizeof(xk)/sizof(double);
//ifft
for (int n=0;n<N;n++)
{
data=complex<double>(0.0,0.0);
for(int k=0;k<N;k++)
{
angle=2*M_PI*k*n/N;
data+=idftXk.at(k) * exp(complex<double> (0.0,angle));
// data+=idft_Xk.at(k)* cmath.exp( (0.0,angle));
//data+=idft_Xk.at(k)* exp( 1i*angle );
//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);
}data/=N;
idftxn.push_back( data);}
return idftxn;
}
int main( )
{fstream finp;
std::cout << "\n For input |H(jw)|, fi(w) from file input '1' " ;
std::cout << "\n For manual input |H(jw)|, fi(w) input '0' " ;std::cin>>key;
ofstream fout ("hcoefs1.txt" , ios::ate | ios::app);
fout.close();
if (key=='0' ) { goto label_02;}
if (key=='1' ) { goto label_01;}label_01:;
std::cout << "\n Open afrcoefs.txt for reading AFR coefficients ";
// ifstream finp("afrcoefs.txt", ios::in );
finp.open("afrcoefs.txt", ios::in );
finp>>fs ;
finp>>N ;
Tmax=N*1/fs;
for (k=0 ; k<N;k++)
{
// f[k]=k*fs/N;
//fout<<"\n n="<<k<<" |H(jw)|="<<abs_H[k]<< "; arg(H(jw))= "<< arg_h[k]<<"; f="<<f[k]<<" Hz" ;
finp>> n ;
finp>>abs_H[k] ;
finp>>arg_H[k] ;
// finp>>f[k] ;
arg_H[k]*= M_PI/180;
}
finp.close();
fout.close();
goto label_03;
label_02:;std::cout<<"\n Input fs ,Hz " ;
std::cin>>fs;
std::cout<<" Input number of points " ;
std::cin>>N ;
Tmax=N*1/fs;
std::cout<<"\n Tmax ="<<Tmax<<" s" ;
for (k=0 ; k<N;k++)
{
f[k]=k*fs/N;//fout<<"\n n="<<k<<" |H(jw)|=" << abs_H[k]<<" ; arg(H(jw))="<< arg_h[k]<<"; f="<<f[k]<<" Hz" ;
//fout<<"\n n=%d |H(jw)|= %lf ; arg(H(jw))= %lf; f=%lf Hz",k, abs_H[k] ,arg_h[k], f[k]) ;
std::cout<<"\n F="<<f[k]<<" Hz" ;
std::cout<<"\n k="<<k<<"; |H(jw)|[k]=" ;
std::cin>> abs_H[k] ;
std::cout<<" k="<<k<<" ; arg(H(jw))[k]=" ;
std::cin>>arg_H[k] ;
//abs_H[k]=1;
//arg_H[k]=0;
arg_H[k]*= M_PI/180;}
label_03:
for (k=0 ; k<N;k++)
{
f[k]=k*fs/N;
std::cout << "\n n="<<k<<";|H(jw)|= "<<abs_H[k]<< " ;arg(H(jw))="<<arg_H[k]*180/M_PI<<" deg;f="<< f[k]<<" Hz";
}
std::cout<< "\n Input 0 for continue " ;
std::cin.get();
std::cout<<"\n OK. Parsing ..." ;
// H(jw)=integr(0;inf;h(tau)*exp(-j*w*tau) ; d tau)
//K=2*3.1415926/fs;
/*
DFT
X[k]=sum(n=0;N-1; x[n]*exp(-2*pi*j*k*n/N) ), k=0,...,N-1
IDFT
x[n]=(1/N)*sum(k=0;N-1;X[k]*exp(2*pi*k*n/N) ), n=0,...,N-1
100*10^-6*500
*///label1:;
//cidft
//printf("\n Open hcoefs1.txt for writing coeffs h[n] ");
cplxvector idft_xn, idft_Xk;
dcomp data;for(int k=0;k<N;k++)
{
Re_H[k]=abs_H[k]*cos(arg_H[k]);
Im_H[k]=abs_H[k]*sin(arg_H[k]);
data =complex<double> (Re_H[k] ,Im_H[k]) ;
// std::cout<<"\n k="<<k<< " Xk[k]= " <<data ; /*real(idft_xn.at(n) ) */ ;
idft_Xk.push_back(data);
}
std::cout <<"\n IDFT\n";
idft_xn= CIDFT(idft_Xk,N);fstream fout1 ("hcoefs1.txt" , ios::ate | ios::app);
for (int n=0;n<N;n++)
{
std::cout<<"\n n="<<n<< " h[n]= " <<real(idft_xn.at(n) ) ;
fout1 << "\n"<<n <<" "<< real(idft_xn.at(n) ) ;
}
fout1.close();//fclose (fp1);
std::cout<< "\n Re(h(n)=h[n]), use it to load into FIR " ;
//label3:;
std::cout<<"\n Input 0 for exit " ;
std::cin.get();
return 0;
}//stdafx.cpp
// stdafx.cpp : source file that includes just the standard includes
// afrtoht.pch will be the pre-compiled header
// stdafx.obj will contain the pre-compiled type information#include "stdafx.h"
//stdafx.h
// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//#pragma once
#define WIN32_LEAN_AND_MEAN // Exclude rarely-used stuff from Windows headers
// TODO: reference additional headers your program requires here//afrcoefs.txt
44100
4
0 1 0
1 1 0
2 0 0
3 0 0//hcoefs1.txt
0 0.5
1 0.25
2 0
3 0.25- Edited by USERPC01 Tuesday, February 6, 2018 6:22 PM fix bugs
Tuesday, February 6, 2018 6:20 PM -
Example of using <vector>, <complex> for Hilbert transform
(fix bugs)
// cplxarray.cpp : main project file.
#include "stdafx.h"
#include <iostream>
//#include <tgmath>//for the type generate marcros.#include <complex>
#include <vector>
#include <set>
#include <cmath>
#include <limits>
/*
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <limits>
#include <iostream>
*/
using namespace std;/* Define complex double data type */
typedef complex<double> dcomp;//std::vector<std::vector<dcomp> >
typedef std::vector<dcomp> cplxvector ;#define NMAX 2048
#define M_PI 3.1415926535897932384626433832795
/*
typedef struct {double re[NMAX];
double im[NMAX];
} CData;typedef struct {
double data[NMAX];
} Array;*/
cplxvector cdft( cplxvector xn ,int N )
{
cplxvector Xk ;
dcomp data;
int k,n ;
double angle;
for (k=0;k<N;k++)
{
//Xk.re[k]=0;
//Xk.im[k]=0;
data=(0.0,0.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);
//data+=xn.at(n)*cmath.exp( (0.0,angle));data+=xn.at(n)* exp( (0.0,angle));
}}
return Xk;
}cplxvector cidft(cplxvector idft_Xk, int N )
{int k,n ;
double angle;
cplxvector idft_xn;
dcomp data;
//N=sizeof(xk)/sizof(double);
//ifft
for (n=0;n<N;n++)
{
data=(0.0,0.0);
for(k=0;k<N;k++)
{
angle=2*M_PI*k*n/N;
//creal(z1),
//cimag(z1),
// data+=idft_Xk.at(k)* cmath.exp( (0.0,angle));
data+=idft_Xk.at(k)* exp( (0.0,angle));
//data+=idft_Xk.at(k)* exp( 1i*angle );
//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);
}data/=N;
idft_xn.push_back( data);
// idft_xn.im[n]=idft_xn.im[n]/N;
// idft_xn.re[n]=idft_xn.re[n]/N;}
return idft_xn ;
}
cplxvector ScaleDFT( cplxvector xk, int num)
{int i;
for (i=0; i<num ; i++)
{
xk.at(i)/= num;
// xk.im[i]=xk.im[i]/num;
}
return xk;
}cplxvector ScaleIQidft( cplxvector Xk, int num)
{
int i;
for (i=0; i<num ; i++)
{
Xk.at(i)*=num;
//Xk.re[i]=Xk.re[i]*num;
//Xk.im[i]=Xk.im[i]*num;
}return Xk;
}
cplxvector hilbert(double *xr , int N)
{
int k, i, no2 ;
//double angle;
double h[NMAX];
cplxvector dft__xn;
//double dft__xn_re[NMAX];
//double dft__xn_im[NMAX];
cplxvector dft__Xk;
//double dft__Xk_re[NMAX];
//double dft__Xk_im[NMAX];
cplxvector idft__XK;
// double idft__XK_re[NMAX];
// double idft__XK_im[NMAX];
cplxvector idft__XN;
if (N==0)
{
printf("\n Error");
idft__XN.push_back((0,0)) ;
return idft__XN;
}
//xnreal=real(xr);
no2=(int) floor ((double)N/2);
// dcomp tmpdata;for (i=0;i<N;i++)
{
// tmpdata=(xr[i],0) ;
dft__xn.push_back((xr[i],0)) ;
//dft__xn_re[i]=xr[i];
//dft__xn_im[i]=0;
}//fft
// dft__Xk=dft(dft__xn,n);
dft__Xk=cdft( dft__xn , N );//X={Xk.re,Xk.im};
for (i=0;i<=N;i++)
{
h[i]=0;
}
if((N>0)&&(( 2*floor((double)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.at(i)*=h[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
idft__XN=cidft(idft__XK, N );
/*
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;
}
cplxvector getanalyticsignal(double *xn, int num)
{
cplxvector 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);
dcomp data1;
double Ximg;
for (i=0; i<num; i++)
{
data1=Hilb.at(i) ;
Ximg=data1.imag();
data1=(xn[i],Ximg );
x_n.push_back(data1);
//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) ;
}
*/cplxvector get_mod_fi(cplxvector x, int n, int mode )
{
int i;
cplxvector res;
double fi;
double module;
for (i=0;i<n;i++)
{
module=abs(x.at(i));
fi=arg(x.at(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) { fi=fi*180/M_PI; }
res.push_back((module,fi));
std::cout<< "\n f="<< " Hz;|Xf(jw)|=" << module << ";arg(Xf)="<< fi<<" deg " ;
}
return res;
}
cplxvector GetSpectrum(double *xn, int num)
{
cplxvector xn_iq, Xf_IQ, Xf;
xn_iq=getanalyticsignal( xn, num);
// Xf_IQ=cdft(xn_iq ,num );
//Xf_IQ=ScaleDFT(Xf_IQ,num);
// Xf =get_mod_fi(Xf_IQ, num,1 );
return Xf;
}
/*
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];
cplxvector Xf;char key;
int n,i,N1;
//int k;
double fs,T, f;
//double x;
//printf ("\n For input x(t) from file input '1' ");
printf ("\n For manual input x(t) input '0' ");
std::cin>>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) ;
std::cin>>fs; //scanf( "%le",&fs);
printf("\n Input number of samples of xin[n] ") ;
std::cin>>N1; //scanf( "%d",&N1);for (n=0 ; n<N1;n++)
{
T=n*1/(N1*fs);
std::cout<<"\n n="<<n<<"; t="<<T<<" s; xin[n]= " ;
std::cin>>xin[n];
}// label_07:;
std::cout<< "\n fs=" <<fs << " Hz " ;
std::cout<< "\n Nsamples= "<<N1 ;for (n=0 ; n<N1;n++)
{
T=n*1/(N1*fs);
std::cout<< "\n n="<< n << "; x[n]=" <<xin[n] << ", t="<<T<<" s" ;
}std::cout<<"\n Input 0 for continue " ;
std::cin.get();
std::cout<<"\n OK. " ;
/*****************************///label1:;
//Xf.re= |Xf|, Xf.im -arg(Xf) (only in this program)
std::cout<<"\n Obtaining spectrum of x[n]" ;
Xf=GetSpectrum(xin, N1);std::cout<< "\n Spectrum of x[n] : " ;
//if ((fp2=fopen("spectrum.txt","a"))==NULL) printf ("\n File opening error");
//fprintf(fp2,"\n %le %d \n",fs,N1);
//fprintf(fp2, "\n n;\t f, Hz;\t |Xf(jw)| ;\t arg(Xf), deg ");
for (i=0; i<N1;i++)
{
f=fs*i/N1;
//cout <<"\n f="<< f<<" Hz;|Xf(jw)|="<<Xf.re[i]<<";arg(Xf)="<<Xf.im[i]<<" deg ";
// std::cout<< "\n f="<<f<<" Hz;|Xf(jw)|=" << real(Xf.at(i))<< ";arg(Xf)="<< imag(Xf.at(i))<<" deg " ;
//fprintf(fp2,"\n %d %le %le",i ,Xf.re[i],Xf.im[i]);}
std::cout<< "\n input '0' for exit " ;
std::cin.get();return 0;
}
Tuesday, February 6, 2018 6:25 PM -
// RectAIMSpectr.cpp :
#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>using namespace std;
#define M_PI 3.1415926535897932384626433832795
double tau0, deltatau, A0;
double f0, Fm,f1;
double w0,T0,m,m1,Wm;
int k,n;#define kmax 60 //change
#define nmax 60 //changedouble U0;
double UFm[nmax], Ffm[nmax];
double fi[kmax], Ufi[kmax];
double fk[2*nmax+1][kmax+1], Ufk[2*nmax+1][kmax+1];
double Q0;
double Q;
int kmax1;void InputData3()
{
//Q0=T0/tau0;
//T0=1/f0
//tau0=T0/Q0;std::cout<<"\nInput A0,V : (1) ";
std::cin>>A0;
std::cout<<"\nInput f1,Hz : (16384) ";
std::cin>>f1;std::cout<<"\nInput Q : (10) ";
std::cin>>Q;
std::cout<<"\nInput kmax : (60) ";
std::cin>>kmax1;
return;
}
/*
________ _______
| | |
| |___________|*/
void getFreqRectQ(void)
{double U0=A0/Q;
double timp=1/(f1*Q);
double w1=2*M_PI*f1;
double Ti=1/f1;
std::cout <<"\n f0="<<0<<"Hz ; U0=" <<U0<<" V ";for(k=1;k<=kmax1; k++)
{
//Ufi~ A0*|sinc(k*pi/q)| =A0*sin (k*pi/q) /(k*pi/q)
//Ufi[k]=2*A0*sin(k*M_PI/Q)/(k*M_PI);
//fi[k]=k*f1;//std::cout <<"\n k="<<k<<"; f[k]="<<fi[k]<<"Hz ; Um[k]=" <<Ufi[k]<<" V ";
double Um;
//Um=fabs( ( 2*A0 /(k*M_PI))*sin(k*M_PI/Q));
//std::cout <<"\n k="<<k<<"; f[k]="<<k*f1<<"Hz ; Um[k]=" <<Um<<" V ";
Um=(2*timp/Ti)*A0*(sin(k*2*M_PI*f1*timp/2)/(k*2*M_PI*f1*timp/2));
Um=fabs(Um);
std::cout <<"\n k="<<k<<"; f[k]="<<k*f1<<"Hz ; Um[k]=" <<Um<<" V ";}
std::cin.get() ;
return;
}/********************************************/
/*S(t)=A0+sum(from k=1;to k=m;{ Ak*cos(k*wi*t)})=
=(taui*Ui/Ti)+(2*taui/Ti*Ui)*sum(from k=1;to k=m;{
(sin(k*wi*taui/2)/(k*wi*taui/2))*cos(kwit))uaim(t)=(taui*Ui/Ti)+(taui/Ti)*Ui*maim*cos(W*t)+
+(2*taui/Ti)*Ui*sum(from k=1;to k=m;{
(sin(k*wi*taui/2)/(k*wi*taui/2))*(1+maim*cos(W*t))*cos(k*wi*t) })Uk*cos(k*wi*t)+0.5*maim*Uk*cos(k*wi-W)*+0.5*maim*Uk*cos(k*wi+W)*t
Uk=(2*taui/Ti)*Ui*(sin(k*wi*taui/2)/(k*wi*taui/2))f=0 Um[0]= taui*Ui/Ti
f=W/2pi Um[W]=taui*miam*Ui/Tif=k*wi-W Um=0.5*maim*Uk
f=k*wi Uk=Uk
f=k*wi+W Um=0.5*maim*Uk
...
Uenv_level[2*pi/taui]=0fsr=Whi
*/
double maim;
void InputData4()
{
//Q0=T0/tau0;
//T0=1/f0
//tau0=T0/Q0;std::cout<<"\nInput A0,V : (1) ";
std::cin>>A0;
std::cout<<"\nInput f1,Hz : (16384) ";
std::cin>>f1;std::cout<<"\nInput Q : (10) ";
std::cin>>Q;
std::cout<<"\nInput Fmod ,Hz : (1000) ";
std::cin>>Fm;
std::cout<<"\nInput maim,% : (10) ";
std::cin>>maim;
maim/=100;
std::cout<<"\nInput kmax : (60) ";
std::cin>>kmax1;
return;
}
void getFreqRectAIM(void)
{
double timp=1/(f1*Q);
double w1=2*M_PI*f1;
double Ti=1/f1;
double U0=timp*A0/Ti;
std::cout <<"\n ; f[0]="<<0<<"Hz ; Um[k]=" << U0<<" V ";
double Fmod=Fm;
double Ufm=timp*maim*A0/Ti;
std::cout <<"\n ; f[Fm]="<<Fmod<<"Hz ; Um[k]=" << Ufm<<" V ";
double Uk,Ulsbk,Uusbk;
double fk,flsbk,fusbk;
std::cout <<"\n ; f[2*pi/taui /( 2*pi)]="<< 1/timp<<"Hz ; Um[1/tau]=" << 0<<" V ";
for (int k=1;k<kmax1;k++){
Uk=(2*timp/Ti)*A0*(sin(k*2*M_PI*f1*timp/2)/(k*2*M_PI*f1*timp/2));
Uk=fabs(Uk);
Ulsbk=0.5*maim*Uk;
Uusbk=0.5*maim*Uk;
fk=k*f1;
flsbk=fk-Fm;
fusbk=fk+Fm;
std::cout <<"\n ; f[lsbk]="<<flsbk<<"Hz ; Um[flsbk]=" <<Ulsbk<<" V ";
std::cout <<"\n ; f[k]= "<<fk<<" Hz ; Um[fk]=" <<Uk<<" V ";
std::cout <<"\n ; f[usbk]="<<fusbk<<"Hz ; Um[fusbk]=" <<Uusbk<<" V ";
}std::cin.get() ;
return;
}int main( )
{
int key;
std::cout<<"\n Input 1 for rect. signal\n Input 2 for AIM ";
std::cin>>key;
switch(key)
{
case 1:
std::cout<<"\n Rect. signal with Q ";
InputData3();
getFreqRectQ();
break;
case 2:
InputData4();
getFreqRectAIM( );
break;
}std::cout<<"\n Input 0 for exit ";
std::cin.get() ;
return 0;
}Tuesday, February 6, 2018 9:12 PM