ソースコードは、テストのためだけです。
// test.c #include <stdio.h> #include <math.h> #include <complex.h> #include <stdlib.h> #define PI 3.1415926535 int main () { const int nhilbert = 51; const double normalization = 1.0028050845; double hilbert[51] = { -1.244274239685367e-02 , // 0 +8.812508503481302e-16 , // 1 -9.007745694645722e-03 , // 2 -4.892921491956048e-16 , // 3 -1.224827000343725e-02 , // 4 -2.457355030119336e-16 , // 5 -1.626770260654829e-02 , // 5 -1.076580799233997e-15 , // 7 -2.126259131785366e-02 , // 8 +1.551646837953304e-15 , // 9 -2.754013546165147e-02 , // 10 -8.271548913440781e-16 , // 11 -3.555144114893793e-02 , // 12 +1.538194937222543e-16 , // 13 -4.616069183085086e-02 , // 14 -2.517547480270933e-16 , // 15 -6.087912031697079e-02 , // 16 +1.529522013663762e-16 , // 17 -8.310832605502005e-02 , // 18 -7.153673299041007e-16 , // 19 -1.216347994939241e-01 , // 20 -3.571474290448791e-17 , // 21 -2.087518418318809e-01 , // 22 -3.409056833722323e-16 , // 23 -6.354638176591967e-01 }; // 24 double wave0[1024], wave1[1024], wave2[1024]; double complex cwave1[1024], cwave2[1024], cbfo[1024]; double fsample = 32000.0; double fsignal = 1000.0; double fbfo = 1010.0; double amp = 1.0; double phase = (45.0/360.0)*2.0*PI; hilbert[(nhilbert-1)/2] = 0.0; // 25 for (int i=((nhilbert-1)/2)+1;i<nhilbert;i++) { hilbert[i] = - hilbert[(nhilbert-1)-i]; // 26 <- -24 } for (int i=0;i<nhilbert;i++) { printf("a+ %20.10f \n", hilbert[i]); } for (int i=0;i<1024;i++) { if( (i>=100) && (i<150) ) { amp = (i-100)/50.0; } else if ( (i>=120) && (i<800) ) { amp = 1.0; } else if ( (i>=800) && (i<850) ) { amp = 1.0 - (i-800)/50.0; } else { amp = 0.0; } cbfo [i] = cos(i*2.0*PI/(fsample/fbfo)+phase) - I * sin(i*2.0*PI/(fsample/fbfo)+phase); wave0[i] = amp * sin(i*2.0*PI/(fsample/fsignal)); wave1[i] = 0.0; wave2[i] = 0.0; printf("b+ %20.10f \n", wave0[i]); } for (int i=nhilbert;i<1024-nhilbert;i++) { wave1[i] = wave0[i-((nhilbert-1)/2)]; printf("c+ %20.10f \n", wave1[i]); } for (int i=nhilbert;i<1024-nhilbert;i++) { for(int j=0;j<nhilbert;j++) { wave2[i] += wave0[i-j] * hilbert[j]; } wave2[i] /= normalization; printf("d+ %20.10f \n", wave2[i]); } for (int i=0;i<1024;i++) { cwave1[i] = wave1[i] + I * wave2[i]; } for (int i=0;i<1024;i++) { cwave2[i] = cwave1[i] * cbfo[i]; printf("e+ %20.10f \n", creal(cwave2[i])); } for (int i=0;i<1024;i++) { printf("f+ %20.10f \n", cimag(cwave2[i])); } for (int i=0;i<1024;i++) { printf("g+ %20.10f \n", cabs(cwave2[i])); } return 0; }
bash.sh
grep a+ ttt | colrm 1 2 > Hilbert grep b+ ttt | colrm 1 2 > Waveform grep c+ ttt | colrm 1 2 > Re grep d+ ttt | colrm 1 2 > Im grep e+ ttt | colrm 1 2 > Sig-Re grep f+ ttt | colrm 1 2 > Sig-Im grep g+ ttt | colrm 1 2 > Sig-Abs
% gcc test.c -o test -lm && ./test > ttt && ./bash.sh % gnuplot gnuplot> plot "Re" w lp, "Sig-Re" w lp, "Sig-Im" w lp, "Sig-abs" w lp