morse2txt: fir_win.c Source File

00001 /*
00002  * $Id: fir_win.c,v 1.16 2009/05/30 00:45:51 kprox Exp $
00003  *
00004  * Copyright (C) 2005:
00005  * Ken Prox <kprox@users.sourceforge.net>
00006  *
00007  * This program is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2 of the License, or
00010  * (at your option) any later version.
00011  *
00012  * This program is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this program; if not, write to the Free Software
00019  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA
00020  *
00021  * --- This object implements a FIR Filter Design Based on the Windows Method ---
00022  *
00023  * Example:
00024  *
00025  * Let's say you want to apply a bandpass filter to `double x[128]'
00026  *
00027  * p_fwin fw;
00028  * fw = fwin_create(8000.0,
00029                     200.0,
00030                     300.0,
00031                     FW_BANDPASS,
00032                     FW_NONE,
00033                     16);
00034  * fwin_apply(fw, x, 128);
00035  * fwin_destroy(fw);
00036  *
00037  * To test this object:
00038  * 
00039  * gcc -D TEST_FIR_WIN -o fir_win fir_win.c -lm 
00040  * with profiling --> add -pg -g
00041  *
00042  * With Discrete time window support, put this in a simple Makefile:
00043  * 
00044  *GCC=gcc
00045  *
00046  *all:
00047  *       $(GCC) -DTEST_FIR_WIN -c fir_win.c -g -Wall -O2
00048  *       $(GCC) -c dt_window.c -g -Wall -O2
00049  *       $(GCC) -lm -g -O2 -o fir_win.exe fir_win.o dt_window.o
00050  */
00051 
00052 #include <stdlib.h>
00053 #ifdef TEST_FIR_WIN
00054     #include <stdio.h>
00055 #endif
00056 #include <string.h>
00057 #include <math.h>
00058 #include <stdint.h>
00059 
00060 #include "fir_win.h"
00061 
00062 #if (FWIN_USE_WINDOW_FUNCTION == 1) 
00063     #include "dt_window.h"
00064 #endif    
00065 
00066 #define FW_2PI              (2.0 * M_PI)
00067 #define DOUBLE_Z_FIR        2
00068 
00072 typedef struct ts_fwin {
00074     double                  sample_freq;
00076     double                  f_low;
00078     double                  f_high;
00080     te_fir_win_filter_type  filter_type;
00082     te_fir_win_window_type  window_type;
00084     double                  *z;
00086     int32_t                 z_state;
00088     double                  *h;
00090     uint32_t                taps;
00092 #if (FWIN_USE_WINDOW_FUNCTION == 1) 
00093     dt_window               h_window;
00094 #endif    
00095 } ts_fwin;
00096 
00097 #if (FWIN_USE_WINDOW_FUNCTION == 1) 
00098 static const te_dt_window_type FW_WINDOW_TABLE[FW_NUMBER_OF_WINDOW_TYPES] = {
00099     0,                      /* FW_NONE */
00100     DT_HAMMING_WINDOW,      /* FW_HAMMING */
00101     DT_HANNING_WINDOW,      /* FW_HANNING */
00102     DT_BLACKMAN_WINDOW      /* FW_BLACKMAN */
00103 };
00105 static void fwin_apply_window(p_fwin fw);
00106 static void fwin_reapply_window(p_fwin fw);
00107 #else
00108     #define fwin_apply_window(a)
00109     #define fwin_reapply_window(a)
00110 #endif
00111 
00112 typedef void ( * fwin_filter_function)(p_fwin fw);
00113 
00114 /* Computes H[] for a bandpass FIR filter using the Window method. */
00115 static void band_pass(p_fwin fw);
00117 static void band_stop(p_fwin fw);
00119 static void low_pass(p_fwin fw);
00121 static void high_pass(p_fwin fw);
00122 /* Computes H[] for a averaging FIR filter. H[] = { 1/taps, ... } */
00123 static void average(p_fwin fw);
00124 
00125 static const fwin_filter_function FW_FILTER_FUNCTIONS[FW_NUMBER_OF_FILTER_TYPES] = {
00126     band_pass,      /* FW_BANDPASS */
00127     band_stop,      /* FW_BANDSTOP */
00128     low_pass,       /* FW_LOWPASS  */
00129     high_pass,      /* FW_HIGHPASS */
00130     average         /* FW_AVERAGE  */
00131 };
00132 
00133 p_fwin fwin_create(double sample_freq,
00134                    double f_low,
00135                    double f_high,
00136                    te_fir_win_filter_type filter_type,
00137                    te_fir_win_window_type window_type,
00138                    uint32_t taps)
00139 {
00140     p_fwin fw;
00141 
00142     if( (fw = malloc(sizeof(ts_fwin))) != NULL) {
00143         #if (FWIN_USE_WINDOW_FUNCTION == 1) 
00144         fw->h_window = NULL;
00145         #endif    
00146         fw->taps = taps;
00147         fw->window_type = window_type;
00148         fw->filter_type = filter_type;
00149         fw->sample_freq = sample_freq;
00150         fw->f_low = f_low;
00151         fw->f_high = f_high;
00152         fw->z_state = 0;
00153         if(
00154             ((fw->h = malloc((taps * sizeof(double)))) != NULL)
00155                 && \
00156             ((fw->z = malloc((DOUBLE_Z_FIR * taps * sizeof(double)))) != NULL)
00157           )  
00158         {
00159             if(filter_type < FW_NUMBER_OF_FILTER_TYPES) {
00160                 FW_FILTER_FUNCTIONS[filter_type](fw);
00161             }
00162             fwin_apply_window(fw);
00163         }
00164         else {
00165             free(fw->h);
00166             free(fw->z);
00167             free(fw);
00168             fw = NULL;
00169         }            
00170     }        
00171     return fw;
00172 }
00173 
00174 void fwin_change_freq(p_fwin fw,
00175                       double f_low,
00176                       double f_high)
00177 {
00178     if(fw != NULL) {
00179         fw->f_low = f_low;
00180         fw->f_high = f_high;
00181         FW_FILTER_FUNCTIONS[fw->filter_type](fw);
00182         fwin_reapply_window(fw);
00183     }        
00184 }
00185 
00186 /* 
00187  * Double delay line FIR filter.
00188  * Applied in double z fashion.
00189  */
00190 void fwin_apply(p_fwin fw, double *data, uint32_t size)
00191 {
00192     uint32_t i, j;
00193     int32_t state;
00194     double accum;
00195     double *p_z, *p_h;
00196         
00197     if(fw != NULL) {
00198         state = fw->z_state;
00199         for(i = 0; i < size; i++) {
00200             /* store input at the beginning of the delay line as well as taps more */
00201             fw->z[state] = fw->z[state + fw->taps] = data[i];
00202             /* calculate the filter */
00203             p_z = fw->z + state; 
00204             p_h = fw->h;
00205             accum = 0.0;
00206             for (j = 0; j < fw->taps; j++) {
00207                 accum += *p_h++ * *p_z++;
00208             }
00209             /* decrement state, wrap if below zero */
00210             state--;
00211             if(state < 0) {
00212                 state += fw->taps;
00213             }
00214             fw->z_state = state;       
00215             data[i] = accum;            
00216         }
00217     }
00218 }
00219 
00220 void fwin_destroy(p_fwin *fw)
00221 {
00222     if(*fw != NULL) {
00223         free((*fw)->h);
00224         free((*fw)->z);
00225 #if (FWIN_USE_WINDOW_FUNCTION == 1) 
00226         free((*fw)->h_window);
00227 #endif    
00228         free(*fw);
00229         *fw = NULL;
00230     }        
00231 }
00232 
00233 static void band_pass(p_fwin fw)
00234 {
00235     uint32_t i;
00236     double w1, w2;
00237     double den, num;
00238     double N = (double)fw->taps;
00239         
00240     w1 = fw->f_low / fw->sample_freq * FW_2PI;  
00241     w2 = fw->f_high / fw->sample_freq * FW_2PI;  
00242     
00243     for(i = 0; i < fw->taps; i++) {
00244         den = M_PI * ((double)i - (N - 1.0) / 2.0);
00245         if(den != 0.0) {
00246             num = sin(w2 * ((double)i - (N - 1.0) / 2.0) ) - \
00247                   sin(w1 * ((double)i - (N - 1.0) / 2.0));
00248             fw->h[i] = num / den;
00249         }            
00250         else {
00251             fw->h[i] = ( w2 - w1 ) / M_PI;
00252         }            
00253     }
00254 }                 
00255 
00256 static void band_stop(p_fwin fw)
00257 {
00258     fw = fw;
00259 }
00260 
00261 static void low_pass(p_fwin fw)
00262 {
00263     fw = fw;
00264 }
00265 
00266 static void high_pass(p_fwin fw)
00267 {
00268     fw = fw;
00269 }
00270 
00271 static void average(p_fwin fw)
00272 {
00273     uint32_t i;
00274     double tmp = 1.0 / (double)fw->taps;
00275         
00276     for(i = 0; i < fw->taps; i++) {
00277         fw->h[i] = tmp;
00278     }
00279 }                 
00280 
00281 #if (FWIN_USE_WINDOW_FUNCTION == 1) 
00282 static void fwin_apply_window(p_fwin fw)
00283 {
00284     if(
00285         (fw->window_type < FW_NUMBER_OF_WINDOW_TYPES) 
00286             && \
00287         (fw->window_type != FW_NONE)
00288       )
00289     {
00290         fw->h_window = dt_window_create((te_dt_window_type)FW_WINDOW_TABLE[fw->window_type], fw->taps);
00291         dt_window_apply(fw->h_window, fw->h);
00292     }       
00293 }
00294 
00295 static void fwin_reapply_window(p_fwin fw)
00296 {
00297     if(fw->window_type != FW_NONE) {
00298         dt_window_apply(fw->h_window, fw->h);
00299     }       
00300 }
00301 #endif
00302 
00303 #ifdef TEST_FIR_WIN
00304 
00305 #define TAPS    12
00306 
00307 /* Desired H with out window function applied. */
00308 static const double BP_TEST_VECTOR[] = {
00309    -2.3592e-02,
00310    -2.3653e-02,
00311    -1.4648e-02,
00312     1.4136e-17,
00313     1.4686e-02,
00314     2.3775e-02,
00315     2.3775e-02,
00316     1.4686e-02,
00317     1.4136e-17,
00318    -1.4648e-02,
00319    -2.3653e-02,
00320    -2.3592e-02
00321 };
00322 
00323 /* Desired H with Blackman window function applied. */
00324 static const double BP_TEST_VECTOR_BM[] = {
00325     3.2740e-19,
00326    -7.7124e-04,
00327    -2.3423e-03,
00328     5.8578e-18,
00329     1.0810e-02,
00330     2.2991e-02,
00331     2.2991e-02,
00332     1.0810e-02,
00333     5.8578e-18,
00334    -2.3423e-03,
00335    -7.7124e-04,
00336     3.2740e-19
00337 };
00338 
00340 static double test[] = {
00341     0.0, 0.0, 0.0, 0.0, 0.0, 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00342     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00343     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00344     100.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00345     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00346     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 100.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00347     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
00348 };    
00349 #define SAMPLES (sizeof(test) / sizeof(double))
00350 
00351 
00352 int main(int argc, char *arg[])
00353 {
00354     int i;
00355     p_fwin fw;
00356 
00357     fw = fwin_create(8000.0,
00358                      750.0,
00359                      850.0,
00360                      FW_BANDPASS,
00361 #if (FWIN_USE_WINDOW_FUNCTION == 1) 
00362                      FW_BLACKMAN,
00363 #else
00364                      FW_NONE,
00365 #endif                     
00366                      TAPS);
00367 
00368     printf("H[] computed, H[] test vector, H[] w/Blackman test vector\n");
00369     for(i = 0; i < TAPS; i++) {
00370         printf("%12.4E, %12.4E, %12.4E\n",
00371                fw->h[i],
00372                BP_TEST_VECTOR[i],
00373                BP_TEST_VECTOR_BM[i]);
00374     }        
00375     fwin_apply(fw, test, SAMPLES);
00376     fwin_change_freq(fw, 100.0, 200.0); 
00377     fwin_apply(fw, test, SAMPLES);
00378     fwin_destroy(&fw);
00379     return 0;
00380 }
00381 #endif /* TEST_FIR_WIN */