FIR filter programming and testing

About FIR (finite impulse response) filter programming in C I have written in my old post. Now I want to show the object-oriented programming approach to solve this task using Java. 

Also there is some improvements: the cyclic buffer is used to convolve filter kernel with  input data samples. It reduces the number of calculations, because there is no need to shift all input data through array where data stored. You just have to put input sample in write place and shift array index (or "pointer" to element).

The FIR filter equation is:

The output sample y(n) is the sum of  products of kernel element h(k) and input element x(n-k). Inside the sum operator the shift through kernel array and buffer of input data is performed. While the shift inside kernel array is straightforward, the shift inside buffer of input data is tricky: if n = 0, then you have to assess to x(0), x(-1), x(-2), ... x(-(N-1)) elements. How cyclic buffer solves this? Look at fig.1.

Fig.1

We have buffer of length 5 and two pointers: iW (for writing) and iR (for reading). These pointers are moved in opposite directions. We put new sample in place where iW points, then make iR equals iW and read data subsequently by incrementing iR. Then we decrement iW by one and repeat the procedure. So, when S0 comes it becomes "0" indexed, but when the S1 comes it becames "0" indexed while S0 becomdes "-1" indexed and so on. 

With this technique the Java implementation of FIR filter is:

package com.blogspot.shulgadim.filters; public class FirFilter{           private int N;      private double h[];      private double y;      private double x[];      private int n;      private int iWrite = 0;      private int iRead = 0;           public FirFilter(double h[]){           this.h = h;           this.N = h.length;                       this.x = new double[N];                                            }      public double filter(double newSample){                  y = 0;                      x[iWrite] = newSample;                 iRead = iWrite;           for(n=0; n<N; n++){                 y += h[n] * x[iRead];                iRead++;                if(iRead == x.length){                     iRead = 0;                }           }           iWrite--;           if(iWrite < 0){                 iWrite = x.length-1;           }                 return y;      }    }

We have class FirFilter where N is filter length, h  is kernel coefficients, y is output sample, x is input data buffer, iWrite and iRead are indexes for x assess. Method filter() performs filtration: it takes new sample and gives filtered sample back. The code works mainly according to the Fig.1. Inside the for loop the cumulative multiplication is performed. And the boundary check is added to preserve pointer to move outside the buffer.

Let's test out filter. The following code does it:

package com.blogspot.shulgadim.filters; public class FirFilterTest {      public static void main(String[] args){                                 int M = 512;       //signal length          double Fs = 250.0;  //sample frequency          double dt = 1.0/Fs; //sample period in time domain                    double t[] = new double[M];  //time vector          for(int i=0; i<M; i++){          t[i] = dt*i;          }                                        double x[]= new double[M]; //input signal          x[0] = 1.0;                                   new Plot("Input signal", t, x, "t, sec", "x, Volts");                                     double h[] = {0.2, 0.2, 0.2, 0.2, 0.2};          FirFilter firFilter = new FirFilter(h);                                       double y[] = new double[M];  //output signal                     for(int i=0; i<M; i++){          y[i] = firFilter.filter(x[i]);          }                                          new Plot("Output signal", t, y, "t, sec", "x, Volts");                  Spectrum spectrum = new Spectrum(y,Fs);                 new Plot("Spectrum", spectrum.getF(),                                         spectrum.getFreqResponse(), "f, Hz", "S, dB");      } }

We create the delta-impulse input signal with sample frequency 250 Hz and send it to the filter. The filter output is impulse response of filter. So the Fourier transform of it is filter frequency response.
If you run test the following windows appears:

The  Plot class draws the signals in separate windows. It uses open source pure Java JFreeChart library. The  Plot code:

package com.blogspot.shulgadim.filters; import java.awt.*; import javax.swing.*; import org.jfree.chart.*; import org.jfree.chart.plot.*; import org.jfree.data.xy.*; public class Plot {      private String name;     private double x[] = null;     private double y[] = null;        private String xName;     private String yName;                 private JFrame frame;                    public Plot(String name, double x[], double y[],                 String xName, String yName) {         this.name = name;         this.x = x;         this.y = y;                this.xName = xName;         this.yName = yName;                frame = new JFrame("Figure");         frame.getContentPane().setLayout(new BorderLayout());                frame.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);         XYDataset dataset = createDataset();         JFreeChart chart = createChart(dataset);         ChartPanel chartPanel = new ChartPanel(chart);            frame.getContentPane().add(chartPanel, BorderLayout.CENTER);         frame.setSize(400, 400);                frame.setVisible(true);            }             private XYDataset createDataset() {         XYSeries series = new XYSeries("");         for(int i=0; i<y.length;i++){             series.add(x[i],y[i]);           }                XYSeriesCollection dataset = new XYSeriesCollection();         dataset.addSeries(series);                                   return dataset;            }             private JFreeChart createChart(final XYDataset dataset) {              final JFreeChart chart = ChartFactory.createXYLineChart(             name,                              // chart title             xName,                      // x axis label             yName,                      // y axis label             dataset,                      // data             PlotOrientation.VERTICAL,             false,                        // include legend             false,                        // tooltips             false                         // urls         );               chart.setBackgroundPaint(Color.LIGHT_GRAY);            XYPlot plot = chart.getXYPlot();            plot.setBackgroundPaint(Color.WHITE)            plot.setDomainGridlinePaint(Color.LIGHT_GRAY);         plot.setRangeGridlinePaint(Color.LIGHT_GRAY);                return chart;            } }

The  Spectrum class does Fast Fourier transform of signal. It uses open source pure Java library JTransforms. The Spectrum code:

package com.blogspot.shulgadim.filters; import java.util.Arrays; import edu.emory.mathcs.jtransforms.fft.DoubleFFT_1D; public class Spectrum{         private double signal[];      private int N;      private double Fs;      private double df;      private double f[];       private double S[];       public Spectrum(double signal[], double Fs){                  this.N = signal.length;                  this.signal = new double[N];           System.arraycopy( signal, 0, this.signal,                              0, signal.length );                   this.Fs = Fs;                  this.df = this.Fs/N;                f = new double[N/2];           for(int i=0; i<N/2; i++){           f[i] = df*i;           }                 S = new double[N/2];           fft();      }         private void fft(){                                          DoubleFFT_1D fft = new DoubleFFT_1D(N);          fft.realForward(signal);                      for(int i=0; i<N/2; i++){          S[i] = Math.sqrt(signal[2*i]*signal[2*i] +                             signal[2*i+1]*signal[2*i+1]);          }          double S2[] = new double[S.length];          System.arraycopy( S, 0, S2, 0, S.length );                 Arrays.sort(S2);                                             for(int i=0; i<N/2; i++){          S[i] = 20*Math.log10(S[i]/S2[S2.length-1]);          }      }              public double[] getFreqResponse(){           return S;      }      public double[] getF(){           return f;      }    }

Few words about Spectrum class code. The fft.realForward() method modifies the data inside the signal array. So, before call fft.realForward() in signal array is time data, after call there is a fft of these data. To prevent data corruption (because arrays are passed by reference) in Spectrum constructor method input signal array copyes to another array - this.signal. After applying fft to get frequency response we have to add squares of odd end even elements of signal array (in that way the real and imaginary parts are stored). After that we find the maximum element in S by Array.sort(S2). And scale the S to plot values in dB.