1 //
2 // firdespm_callback_example.c
3 //
4 // This example demonstrates finite impulse response filter design
5 // using the Parks-McClellan algorithm with callback function for
6 // arbitrary response and weighting function.
7 //
8 // SEE ALSO: firdes_kaiser_example.c
9
10 #include <stdio.h>
11 #include <stdlib.h>
12 #include <math.h>
13 #include "liquid.h"
14
15 #define OUTPUT_FILENAME "firdespm_callback_example.m"
16
17 // user-defined callback function defining response and weights
callback(double _frequency,void * _userdata,double * _desired,double * _weight)18 int callback(double _frequency,
19 void * _userdata,
20 double * _desired,
21 double * _weight)
22 {
23 // de-reference pointer as floating-point value
24 unsigned int n = *((unsigned int*)_userdata);
25 double v = sincf(n*_frequency);
26 double fc = 1.0f / (float)n;
27
28 // inverse sinc
29 if (_frequency < fc) {
30 *_desired = 1.0f / v; // inverse of sinc
31 *_weight = 4.0f;
32 } else {
33 *_desired = 0.0f; // stop-band
34 *_weight = 10*fabs(v) * exp(4.0*_frequency);
35 }
36 return 0;
37 }
38
main(int argc,char * argv[])39 int main(int argc, char*argv[])
40 {
41 // filter design parameters
42 unsigned int n = 8; // sinc filter length
43 unsigned int h_len = 81; // inverse sinc filter length
44 liquid_firdespm_btype btype = LIQUID_FIRDESPM_BANDPASS;
45 unsigned int num_bands = 2;
46 float bands[4] = {0.00f, 0.75f/(float)n, // pass-band
47 1.05f/(float)n, 0.5f}; // stop-band
48
49 // design filter
50 float h[h_len];
51 firdespm q = firdespm_create_callback(h_len,num_bands,bands,btype,callback,&n);
52 firdespm_execute(q,h);
53 firdespm_destroy(q);
54
55 // print coefficients
56 unsigned int i;
57 for (i=0; i<h_len; i++)
58 printf("h(%4u) = %16.12f;\n", i+1, h[i]);
59
60 // open output file
61 FILE*fid = fopen(OUTPUT_FILENAME,"w");
62 fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
63 fprintf(fid,"clear all;\n");
64 fprintf(fid,"close all;\n\n");
65 fprintf(fid,"n=%u;\n", n);
66 fprintf(fid,"h_len=%u;\n", h_len);
67
68 for (i=0; i<h_len; i++)
69 fprintf(fid,"h(%4u) = %20.8e;\n", i+1, h[i]);
70
71 fprintf(fid,"nfft=1024;\n");
72 fprintf(fid,"H0=20*log10(abs(fftshift(fft(ones(1,n)/n,nfft))));\n");
73 fprintf(fid,"H1=20*log10(abs(fftshift(fft(h, nfft))));\n");
74 fprintf(fid,"Hc=H0+H1;\n");
75 fprintf(fid,"f=[0:(nfft-1)]/nfft-0.5;\n");
76 fprintf(fid,"figure;\n");
77 fprintf(fid,"hold on;\n");
78 fprintf(fid,"plot(f,H0,'Color',[0.5 0.5 0.5],'LineWidth',1);\n");
79 fprintf(fid,"plot(f,H1,'Color',[0.0 0.2 0.5],'LineWidth',1);\n");
80 fprintf(fid,"plot(f,Hc,'Color',[0.0 0.5 0.2],'LineWidth',2);\n");
81 fprintf(fid,"hold off;\n");
82 fprintf(fid,"grid on;\n");
83 fprintf(fid,"xlabel('normalized frequency');\n");
84 fprintf(fid,"ylabel('PSD [dB]');\n");
85 fprintf(fid,"legend('sinc','inverse sinc','composite');\n");
86 fprintf(fid,"title('Filter design (firdespm), inverse sinc');\n");
87 fprintf(fid,"axis([-0.5 0.5 -80 20]);\n");
88
89 fclose(fid);
90 printf("results written to %s.\n", OUTPUT_FILENAME);
91
92 printf("done.\n");
93 return 0;
94 }
95
96