Analytic Signals (2)

The source code is just for a test.

// 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