1 //
2 // iirhilb_filter_example.c
3 //
4 // Hilbert transform example. This example demonstrates the
5 // functionality of iirhilbf (infinite impulse response Hilbert transform)
6 // as a filter to remove the negative half of the spectrum.
7 //
8 // SEE ALSO: iirhilb_interp_example.c
9 // iirhilb_example.c
10 //
11
12 #include <stdio.h>
13 #include <complex.h>
14 #include <math.h>
15
16 #include "liquid.h"
17
18 #define OUTPUT_FILENAME "iirhilb_filter_example.m"
19
main()20 int main() {
21 unsigned int order = 7; // Hilbert filter order
22 unsigned int n = 128; // number of input samples
23
24 // derived values
25 unsigned int num_samples = n + 50;
26
27 // create Hilbert transform objects
28 iirhilbf q0 = iirhilbf_create_default(order);
29 iirhilbf q1 = iirhilbf_create_default(order);
30 iirhilbf_print(q0);
31
32 // data arrays
33 float complex x[num_samples]; // complex input
34 float y[num_samples]; // real output
35 float complex z[num_samples]; // complex output
36
37 // run transform
38 unsigned int i;
39 for (i=0; i<num_samples; i++) {
40 // generate input
41 x[i] = 1.0f*cexpf( 0.12f*_Complex_I*2*M_PI*i) + // primary tone
42 0.1f*cexpf( 0.17f*_Complex_I*2*M_PI*i) + // secondary tone
43 0.2f*cexpf(-0.40f*_Complex_I*2*M_PI*i); // tone in negative spectrum
44 x[i] *= (i < n) ? 1.855f*hamming(i,n) : 0.0f;
45
46 // convert to real
47 iirhilbf_c2r_execute(q0, x[i], &y[i]);
48
49 // convert back to complex
50 iirhilbf_r2c_execute(q1, y[i], &z[i]);
51 }
52
53 // destroy Hilbert transform object
54 iirhilbf_destroy(q0);
55 iirhilbf_destroy(q1);
56
57 //
58 // export results to file
59 //
60 FILE*fid = fopen(OUTPUT_FILENAME,"w");
61 fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
62 fprintf(fid,"clear all;\n");
63 fprintf(fid,"close all;\n");
64 fprintf(fid,"n=%u;\n", n);
65 fprintf(fid,"num_samples=%u;\n", num_samples);
66 fprintf(fid,"t = 0:(num_samples-1);\n");
67
68 for (i=0; i<num_samples; i++) {
69 // print results
70 fprintf(fid,"x(%3u) = %12.4e + %12.4ej;\n", i+1, crealf(x[i]), cimagf(x[i]));
71 fprintf(fid,"y(%3u) = %12.4e;\n", i+1, y[i]);
72 fprintf(fid,"z(%3u) = %12.4e + %12.4ej;\n", i+1, crealf(z[i]), cimagf(z[i]));
73 }
74
75 fprintf(fid,"figure;\n");
76 fprintf(fid,"subplot(3,1,1);\n");
77 fprintf(fid," plot(t,real(x),'Color',[0.00 0.25 0.50],'LineWidth',1.3,...\n");
78 fprintf(fid," t,imag(x),'Color',[0.00 0.50 0.25],'LineWidth',1.3);\n");
79 fprintf(fid," legend('real','imag','location','northeast');\n");
80 fprintf(fid," ylabel('transformed/complex');\n");
81 fprintf(fid," axis([0 num_samples -2 2]);\n");
82 fprintf(fid," grid on;\n");
83 fprintf(fid,"subplot(3,1,2);\n");
84 fprintf(fid," plot(t,y,'Color',[0.00 0.25 0.50],'LineWidth',1.3);\n");
85 fprintf(fid," ylabel('original/real');\n");
86 fprintf(fid," axis([0 num_samples -2 2]);\n");
87 fprintf(fid," grid on;\n");
88 fprintf(fid,"subplot(3,1,3);\n");
89 fprintf(fid," plot(t,real(z),'Color',[0.00 0.25 0.50],'LineWidth',1.3,...\n");
90 fprintf(fid," t,imag(z),'Color',[0.00 0.50 0.25],'LineWidth',1.3);\n");
91 fprintf(fid," legend('real','imag','location','northeast');\n");
92 fprintf(fid," ylabel('transformed/complex');\n");
93 fprintf(fid," axis([0 num_samples -2 2]);\n");
94 fprintf(fid," grid on;\n");
95
96 // plot results
97 fprintf(fid,"nfft=4096;\n");
98 fprintf(fid,"%% compute normalized windowing functions\n");
99 fprintf(fid,"X=20*log10(abs(fftshift(fft(x/n,nfft))));\n");
100 fprintf(fid,"Y=20*log10(abs(fftshift(fft(y/n,nfft))));\n");
101 fprintf(fid,"Z=20*log10(abs(fftshift(fft(z/n,nfft))));\n");
102 fprintf(fid,"f =[0:(nfft-1)]/nfft-0.5;\n");
103 fprintf(fid,"figure; plot(f,X,'LineWidth',1,'Color',[0.50 0.50 0.50],...\n");
104 fprintf(fid," f,Y,'LineWidth',2,'Color',[0.00 0.50 0.25],...\n");
105 fprintf(fid," f,Z,'LineWidth',1,'Color',[0.00 0.25 0.50]);\n");
106 fprintf(fid,"grid on;\n");
107 fprintf(fid,"axis([-0.5 0.5 -80 20]);\n");
108 fprintf(fid,"xlabel('normalized frequency');\n");
109 fprintf(fid,"ylabel('PSD [dB]');\n");
110 fprintf(fid,"legend('original/cplx','transformed/real','regenerated/cplx','location','northeast');");
111
112 fclose(fid);
113 printf("results written to %s\n", OUTPUT_FILENAME);
114
115 printf("done.\n");
116 return 0;
117 }
118