Main Archive Specials Wiki | FAQ Links Submit Forum


Windowed Sinc FIR Generator

Type : LP, HP, BP, BS
References : Posted by Bob Maling
Linked file : wsfir.h

Notes :
This code generates FIR coefficients for lowpass, highpass, bandpass, and bandstop filters by windowing a sinc function.

The purpose of this code is to show how windowed sinc filter coefficients are generated. Also shown is how highpass, bandpass, and bandstop filters can be made from lowpass filters.

Included windows are Blackman, Hanning, and Hamming. Other windows can be added by following the structure outlined in the opening comments of the header file.





Comments


Added on : 14/04/05 by Christian[ AT ]savioursofsoul[ DOT ]de
Comment :
// Object Pascal Port...

unit SincFIR;

(* Windowed Sinc FIR Generator
   Bob Maling (BobM.DSP@gmail.com)
   Contributed to musicdsp.org Source Code Archive
   Last Updated: April 12, 2005
   Translated to Object Pascal by Christian-W. Budde

   Usage:
   Lowpass:wsfirLP(H, WindowType, CutOff)
   Highpass:wsfirHP(H, WindowType, CutOff)
   Bandpass:wsfirBP(H, WindowType, LowCutOff, HighCutOff)
   Bandstop:wsfirBS(H, WindowType, LowCutOff, HighCutOff)

   where:
   H (TDoubleArray): empty filter coefficient table (SetLength(H,DesiredLength)!)
   WindowType (TWindowType): wtBlackman, wtHanning, wtHamming
   CutOff (double): cutoff (0 < CutOff < 0.5, CutOff = f/fs)
   --> for fs=48kHz and cutoff f=12kHz, CutOff = 12k/48k => 0.25

   LowCutOff (double):low cutoff (0 < CutOff < 0.5, CutOff = f/fs)
   HighCutOff (double):high cutoff (0 < CutOff < 0.5, CutOff = f/fs)

   Windows included here are Blackman, Blackman-Harris, Gaussian, Hanning
   and Hamming.*)

interface

uses Math;

type TDoubleArray = array of Double;
     TWindowType = (wtBlackman, wtHanning, wtHamming, wtBlackmanHarris,
                    wtGaussian); // Window type contstants

// Function prototypes
procedure wsfirLP(var H : TDoubleArray; const WindowType : TWindowType; const CutOff : Double);
procedure wsfirHP(var H : TDoubleArray; const WindowType : TWindowType; const CutOff : Double);
procedure wsfirBS(var H : TDoubleArray; const WindowType : TWindowType; const LowCutOff, HighCutOff : Double);
procedure wsfirBP(var H : TDoubleArray; const WindowType : TWindowType; const LowCutOff, HighCutOff : Double);
procedure genSinc(var Sinc : TDoubleArray; const CutOff : Double);
procedure wGaussian(var W : TDoubleArray);
procedure wBlackmanHarris(var W : TDoubleArray);
procedure wBlackman(var W : TDoubleArray);
procedure wHanning(var W : TDoubleArray);
procedure wHamming(var W : TDoubleArray);

implementation

// Generate lowpass filter
// This is done by generating a sinc function and then windowing it
procedure wsfirLP(var H : TDoubleArray; const WindowType : TWindowType; const CutOff : Double);
begin
genSinc(H, CutOff);     // 1. Generate Sinc function
case WindowType of  // 2. Generate Window function -> lowpass filter!
  wtBlackman: wBlackman(H);
  wtHanning: wHanning(H);
  wtHamming: wHamming(H);
  wtGaussian: wGaussian(H);
  wtBlackmanHarris: wBlackmanHarris(H);
end;
end;

// Generate highpass filter
// This is done by generating a lowpass filter and then spectrally inverting it
procedure wsfirHP(var H : TDoubleArray; const WindowType : TWindowType; const CutOff : Double);
var i : Integer;
begin
wsfirLP(H, WindowType, CutOff); // 1. Generate lowpass filter

// 2. Spectrally invert (negate all samples and add 1 to center sample) lowpass filter
// = delta[n-((N-1)/2)] - h[n], in the time domain
for i:=0 to Length(H)-1
  do H[i]:=H[i]*-1.0;
H[(Length(H)-1) div 2]:=H[(Length(H)-1) div 2]+1.0;
end;

// Generate bandstop filter
// This is done by generating a lowpass and highpass filter and adding them
procedure wsfirBS(var H : TDoubleArray; const WindowType : TWindowType; const LowCutOff, HighCutOff : Double);
var i  : Integer;
    H2 : TDoubleArray;
begin
SetLength(H2,Length(H));

// 1. Generate lowpass filter at first (low) cutoff frequency
wsfirLP(H, WindowType, LowCutOff);

// 2. Generate highpass filter at second (high) cutoff frequency
wsfirHP(H2, WindowType, HighCutOff);

// 3. Add the 2 filters together
for i:=0 to Length(H)-1
  do H[i]:=H[i]+H2[i];

SetLength(H2,0);
end;

// Generate bandpass filter
// This is done by generating a bandstop filter and spectrally inverting it
procedure wsfirBP(var H : TDoubleArray; const WindowType : TWindowType; const LowCutOff, HighCutOff : Double);
var i : Integer;
begin
wsfirBS(H, WindowType, LowCutOff, HighCutOff); // 1. Generate bandstop filter

// 2. Spectrally invert (negate all samples and add 1 to center sample) lowpass filter
// = delta[n-((N-1)/2)] - h[n], in the time domain
for i:=0 to Length(H)-1
  do H[i]:=H[i]*-1.0;
H[(Length(H)-1) div 2]:=H[(Length(H)-1) div 2]+1.0;
end;

// Generate sinc function - used for making lowpass filter
procedure genSinc(var Sinc : TDoubleArray; const Cutoff : Double);
var i,j  : Integer;
    n,k  : Double;
begin
j:=Length(Sinc)-1;
k:=1/j;
// Generate sinc delayed by (N-1)/2
for i:=0 to j do
  if (i=j div 2)
   then Sinc[i]:=2.0*Cutoff
   else
    begin
     n:=i-j/2.0;
     Sinc[i]:=sin(2.0*PI*Cutoff*n)/(PI*n);
    end;
end;

// Generate window function (Gaussian)
procedure wGaussian(var W : TDoubleArray);
var i,j : Integer;
    k   : Double;
begin
j:=Length(W)-1;
k:=1/j;
for i:=0 to j
  do W[i]:=W[i]*(exp(-5.0/(sqr(j))*(2*i-j)*(2*i-j)));
end;

// Generate window function (Blackman-Harris)
procedure wBlackmanHarris(var W : TDoubleArray);
var i,j : Integer;
    k   : Double;
begin
j:=Length(W)-1;
k:=1/j;
for i:=0 to j
  do W[i]:=W[i]*(0.35875-0.48829*cos(2*PI*(i+0.5)*k)+0.14128*cos(4*PI*(i+0.5)*k)-0.01168*cos(6*PI*(i+0.5)*k));
end;

// Generate window function (Blackman)
procedure wBlackman(var W : TDoubleArray);
var i,j : Integer;
    k   : Double;
begin
j:=Length(W)-1;
k:=1/j;
for i:=0 to j
  do W[i]:=W[i]*(0.42-(0.5*cos(2*PI*i*k))+(0.08*cos(4*PI*i*k)));
end;

// Generate window function (Hanning)
procedure wHanning(var W : TDoubleArray);
var i,j : Integer;
    k   : Double;
begin
j:=Length(W)-1;
k:=1/j;
for i:=0 to j
  do W[i]:=W[i]*(0.5*(1.0-cos(2*PI*i*k)));
end;

// Generate window function (Hamming)
procedure wHamming(var W : TDoubleArray);
var i,j : Integer;
    k   : Double;
begin
j:=Length(W)-1;
k:=1/j;
for i:=0 to j
  do W[i]:=W[i]*(0.54-(0.46*cos(2*PI*i*k)));
end;

end.              




Added on : 05/01/07 by scoofy[ AT ]inf[ DOT ]elte[ DOT ]hu
Comment :
The Hanning window is often incorrecty referred to as 'Hanning', since it was named after a guy called Julius von Hann. So it's more appropriate to call it 'Hann' window.




Added on : 02/09/07 by Dave in sinc land
Comment :
I've seen a BASIC version of the same genSinc code.
Shouldn't a 'sin(2.0*Cutoff)' be used when the divide by zero check is 0?
...
  if (i=j div 2)
   then Sinc[i]:=2.0*Cutoff
   else
    begin
     n:=i-j/2.0;
     Sinc[i]:=sin(2.0*PI*Cutoff*n)/(PI*n);
    end;
...




Added on : 02/09/07 by Dave in sinc land
Comment :
Scrap that, I've just FFT'd the response, and it appears to be correct as it was. Hey it just looked wrong o.k.  ;)



Added on : 07/09/08 by rousicoello[ AT ]gmail[ DOT ]com
Comment :
What do I have to do to apply this windowed sinc filter to a signal "Y" for example... what is the code for this?
Imagine signal "Y" is a stereo song which means that I have Y1 from channel 1 and Y2 from channel 2.. help me please as soon as posible because at least i want to know how to aplly the filter to any signal...




Added on : 24/04/12 by craig101[ DOT ]com101[ AT ]yahoo[ DOT ]com
Comment :
Can anyone, provide a sample code on how to apply this to an array of samples?              



Added on : 10/06/12 by amishman35[ AT ]cox[ DOT ]net
Comment :
Here is a FIR filter that works with up to about 80 coefficients on a Sony PSP running at 280 Megahertz. The multipliers (8191.0 and 32767.0) may need to be incremented by one, or not. You need to call 'fir_filter' with the sample to process in 'sample' and the floating point coefficents in '*coefs' and how long the entire filter is in 'len' and the returned sample is a 16 bit audio (signed short) sample. When it is first run it will normalize the filter.

/* code starts here, remove spaces between the includes */
#include < math.h >
#include < string.h >
#include < stdio.h >
#include < stdlib.h >
              void normalizeCoefs(signed short *firshort, signed short *outfilter, unsigned long len)
{
    // filter gain at uniform frequency intervals
    double *g=NULL;
    double *dtemp=NULL;
    double theta, s, c, sac, sas;
    double gMax = -100.0;
    double sc = 10.0/log(10.0);
    double t = M_PI / len;
    long i=0;
    long n=len-1;
    double normFactor =0;
    g=(double*)malloc((len+1)*sizeof(double));
    dtemp=(double*)malloc((len+1)*sizeof(double));
    for (i=0; i<len; i++)
    {
        dtemp[i]=firshort[i]/32768.0;
    }
    for (i=0; i<len; i++)
    {
        theta = i*t;
        sac = 0.0;
        sas = 0.0;
        long k=0;
        for (k=0; k<len; k++)
        {
            c =cos(k*theta);
            s =sin(k*theta);
            sac += c*(dtemp[k]);
            sas += s*(dtemp[k]);
        }
        g[i] = sc*log(sac*sac + sas*sas);
        if(g[i]>gMax)
        {
            gMax=g[i];
        }      
    }
    // normalise to 0 dB maximum gain
    for (i=0; i<len; i++)
    {
        g[i] -= gMax;
    }
    // normalise coefficients
    normFactor =0;pow(10.0, -0.05*gMax);
    for (i=0; i<len; i++)
    {
        dtemp[i] *= normFactor;
    }
    n=len-1;
    for (i=0; i<len; i++)
    {
        outfilter[i]=dtemp[n--]*32767.0;
    }
    free(dtemp);
    free(g);
  }
  
signed short generalFIR(short input, short *coefs, unsigned long len,unsigned char resetit)
{
    
    static signed short *delay_line=NULL;
    static unsigned int first=1;
    static signed long reacc=0;
    unsigned long filtcnt=0;
    static int consumer=1;
    static int producer=0;
    unsigned long tx=0;
    if(resetit==1)
    {
        first=1;
    }
    if(first==1)
    {
        if(delay_line!=NULL)
            free(delay_line);
        producer=0;
        consumer=1;
        delay_line=(signed short*)malloc(reallen*sizeof(signed short));
        first=0;
    }
    reacc=0;
    delay_line[producer++]=input;
    if(producer>len-1)
    {
        producer-=len;
    }
    int filtptr=consumer;
    for(tx=0;tx<len;tx++)
    {
        reacc+=(delay_line[filtptr++]*coefs[tx]);
        if(filtptr>len-1)
        {
            filtptr-=len;
        }
        
    }
    consumer++;
    if(consumer>len-1)
    {
        consumer-=len;
    }
    reacc = reacc >> 15;
    if(reacc<-32768)
        reacc= -32768;
    if(reacc>32767)
        reacc=32767;
    return (signed short)reacc;
}


signed short fir_filter(signed short sample,double *coefs,unsigned long len)
{
    static unsigned char first=1;
    static signed short *newfir=NULL;
    if(first==1)
    {
        unsigned long i=0;
        if(newfir!=NULL)
        {
            free(newfir);
            newfir=NULL;
        }
        newfir=(signed short*)malloc((filterorder+1) * sizeof(signed short));
        if(newfir==NULL)
        {
            return 0;
        }
        for(k=0;k<len;k++)
        {
            newfir[k]=coefs[k] * 8191.0;
        }
        normalizeCoefs(newfir,newfir,len);
        generalFIR(0,newfir,len,1);
        first=0;
    }
    return generalFIR(sample,newfir,len,0);
}




Add your own comment
Comments are displayed in fixed width, no HTML code allowed!
Email:

Comment:

Are you human?



Site created and maintained by Bram
Graphic design by line.out | Server sponsered by fxpansion