DSP Filter Design

From Hamsterworks Wiki!

Jump to: navigation, search

This is only sort of a FPGA_Project. This code might help you design kernels for low pass DSP filters.

You can control the kernel width, cutoff frequency, window type and the number of significant bits used in each term.


Sample Output

C:\Users\Hamster\Documents\Visual Studio 2010\Projects\sinc_window\Debug>.\sinc_window.exe -w 61 -f 0.15 -t b -b 16
Windowed Sinc filter
  Width              : 61
  Critical frequency : 0.150
  Bits               : 16


   0,  0.00000000000000
   1,  0.00021362304688
   2,  0.00036621093750
   3,  0.00018310546875
   4, -0.00054931640625
   5, -0.00140380859375
   6, -0.00119018554688
   7,  0.00088500976563
   8,  0.00372314453125
   9,  0.00427246093750
  10,  0.00000000000000
  11, -0.00738525390625
  12, -0.01113891601563
  13, -0.00460815429688
  14,  0.01101684570313
  15,  0.02331542968750
  16,  0.01693725585938
  17, -0.01092529296875
  18, -0.04104614257813
  19, -0.04244995117188
  20,  0.00000000000000
  21,  0.06243896484375
  22,  0.08926391601563
  23,  0.03549194335938
  24, -0.08346557617188
  25, -0.17901611328125
  26, -0.13690185546875
  27,  0.09899902343750
  28,  0.46719360351563
  29,  0.80545043945313
  30,  0.94247436523438
  31,  0.80545043945313
  32,  0.46719360351563
  33,  0.09899902343750
  34, -0.13690185546875
  35, -0.17901611328125
  36, -0.08346557617188
  37,  0.03549194335938
  38,  0.08926391601563
  39,  0.06243896484375
  40,  0.00000000000000
  41, -0.04244995117188
  42, -0.04104614257813
  43, -0.01092529296875
  44,  0.01693725585938
  45,  0.02331542968750
  46,  0.01101684570313
  47, -0.00460815429688
  48, -0.01113891601563
  49, -0.00738525390625
  50,  0.00000000000000
  51,  0.00427246093750
  52,  0.00372314453125
  53,  0.00088500976563
  54, -0.00119018554688
  55, -0.00140380859375
  56, -0.00054931640625
  57,  0.00018310546875
  58,  0.00036621093750
  59,  0.00021362304688
  60,  0.00000000000000


Filter Kernel
-------------------------------------------------------------------
 0.96                               X

 0.90

 0.83                              X X

 0.77

 0.71

 0.65

 0.58

 0.52
                                  X   X
 0.46

 0.40

 0.34

 0.27

 0.21

 0.15
                            X    X     X    X
 0.09                      X                 X
                             X             X
 0.02 XXXXXXXXXXXX-XXXX---X-------------------X---XXXX-XXXXXXXXXXXX
                  X    X                         X    X
-0.04                   XX                     XX
                              X           X
-0.10
                                X       X
-0.16                          X         X

-0.23

-------------------------------------------------------------------


FFT
---------------------------------------------------------------------
   5
     XXXXXXXXXXXXXXXXXX
  -5 XXXXXXXXXXXXXXXXXXXX
     XXXXXXXXXXXXXXXXXXXXX
 -15 XXXXXXXXXXXXXXXXXXXXXX
     XXXXXXXXXXXXXXXXXXXXXX
 -25 XXXXXXXXXXXXXXXXXXXXXXX
     XXXXXXXXXXXXXXXXXXXXXXX
 -35 XXXXXXXXXXXXXXXXXXXXXXXX
     XXXXXXXXXXXXXXXXXXXXXXXX
 -45 XXXXXXXXXXXXXXXXXXXXXXXX
     XXXXXXXXXXXXXXXXXXXXXXXXX
 -55 XXXXXXXXXXXXXXXXXXXXXXXXX
     XXXXXXXXXXXXXXXXXXXXXXXXX
 -65 XXXXXXXXXXXXXXXXXXXXXXXXX
     XXXXXXXXXXXXXXXXXXXXXXXXX
 -75 XXXXXXXXXXXXXXXXXXXXXXXXXX
     XXXXXXXXXXXXXXXXXXXXXXXXXX
 -85 XXXXXXXXXXXXXXXXXXXXXXXXXX  X       XX  X      X X          X
     XXXXXXXXXXXXXXXXXXXXXXXXXXXXX   X   XX  XX   X X XXXX    X XX
 -95 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  XX  XX  XXXXXXXX XXXX  X XXXXXX
     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  XXXXXX XXXXXXXXXXXXXX  XXXXXXXXX
-105 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX XXXXXXXXXXXXXXXXXXXXXXXXX
     XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
-115 XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
---------------------------------------------------------------------
     0.0                                                          0.5


C:\Users\Hamster\Documents\Visual Studio 2010\Projects\sinc_window\Debug>

C source

/************************************************************************
 * windowed_sinc.c - Calculate and display a Windowed Sinc filter kernel
 *
 * Author : Mike Field <hamster@snap.net.nz>
 *
 ************************************************************************/
#include "stdafx.h"        /* For Windows */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <assert.h>
#include <string.h>

/*******************
 * Graphing details 
 */
#define FFT_GRAPH_SCALE      5  /* db per line on the frequency graph */
#define FFT_GRAPH_HEIGHT    25  /* How tall the FFT graph is */
#define KERNEL_GRAPH_HEIGHT 40  /* How tall the kernel graph is */

#define PI 3.141592653589793

/* Begin DFT taken from http://paulbourke.net/miscellaneous/dft/ */
int DFT(int dir,int m,double *x1,double *y1)
{
   long i,k;
   double arg;
   double cosarg,sinarg;
   double *x2=NULL,*y2=NULL;

   x2 = (double *)malloc(m*sizeof(double));
   y2 = (double *)malloc(m*sizeof(double));
   if (x2 == NULL || y2 == NULL)
      return(0);

   for (i=0;i<m;i++) {
      x2[i] = 0;
      y2[i] = 0;
      arg = - dir * 2.0 * 3.141592654 * (double)i / (double)m;
      for (k=0;k<m;k++) {
         cosarg = cos(k * arg);
         sinarg = sin(k * arg);
         x2[i] += (x1[k] * cosarg - y1[k] * sinarg);
         y2[i] += (x1[k] * sinarg + y1[k] * cosarg);
      }
   }

   /* Copy the data back */
   if (dir == 1) {
      for (i=0;i<m;i++) {
         x1[i] = x2[i] / (double)m;
         y1[i] = y2[i] / (double)m;
      }
   } else {
      for (i=0;i<m;i++) {
         x1[i] = x2[i];
         y1[i] = y2[i];
      }
   }

   free(x2);
   free(y2);
   return(1);
}
/* End DFT from http://paulbourke.net/miscellaneous/dft/ */


/*******************************
 * Ouput the FFT specturm 
 ******************************/
static void graphFFT(double *freq, int fft_width)
{
  int i;
  /* Graph and display the FFT spectrum */
  printf("\n\nFFT\n");
  for(i = 0; i < fft_width/2+5; i++)
      putchar('-');
  putchar('\n');

  for(i = 0; i < FFT_GRAPH_HEIGHT; i++)
  {
      int j;
      if((i&1) == 0)
        printf("%4i ",-i*FFT_GRAPH_SCALE + FFT_GRAPH_SCALE);
      else
        printf("     ");

      for(j = 0; j < fft_width/2; j++)
      {
          if(freq[j] > -i*FFT_GRAPH_SCALE + FFT_GRAPH_SCALE/2.0)
              putchar('X');
          else
              putchar(' ');
      }
      putchar('\n');
  }
  for(i = 0; i < fft_width/2+5; i++)
      putchar('-');
  putchar('\n');
  printf("     0.0%*.*s0.5\n", fft_width/2-6, fft_width/2-6, "");
}

/*******************************
 * Graph the filter kernel
 ******************************/
static void graphKernel(double *values, int width)
{
  int i;
  double step, max, min;

  /* Graph and display the FFT spectrum */
  printf("\n\nFilter Kernel\n");
  for(i = 0; i < width+6; i++)
      putchar('-');
  putchar('\n');


  max = values[0];
  min = values[0];

  for(i = 1; i < width; i++) {
     if(max < values[i]) max = values[i];
     if(min > values[i]) min = values[i];
  }
  step = (max - min) / (KERNEL_GRAPH_HEIGHT-4);

  for(i = 0; i < KERNEL_GRAPH_HEIGHT; i++)
  {
     double level = max - i * step + step/2;
      int j;
      if((i&1) == 0)
        printf("%5.2f ",level);
      else
        printf("      ");

      for(j = 0; j < width; j++)
      {
          if(values[j] >= level - step && values[j] < level)
              putchar('X');
          else  if(0.0 >= level - step && 0.0 < level)
              putchar('-');
          else
              putchar(' ');
      }
      putchar('\n');
  }
  for(i = 0; i < width+6; i++)
      putchar('-');
  putchar('\n');
}

/**********************************
 * Calculate the sinc function
 *********************************/
static void sincFunction(double *sinc, int width, double cutoff)
{
  int i;
  for(i = 0; i < width; i++)
  {
    if(i-width/2 == 0)
        sinc[i] = 2 * PI * cutoff;
    else
        sinc[i] = sin(2*PI * cutoff * ((double)i-width/2)) /((double)(i-width/2));
  }
}
/**********************************
 * Calculate the blackman window
 *********************************/
static void blackmanWindow(double *window, int width)
{
  int i;
  /* Calculate the Blackman window */
  for(i = 0; i < width; i++)
    window[i] = 0.42659 - 0.49656 * cos(2 * PI * i / (width-1)) + 0.076849 * cos(4 * PI * i / (width-1));
}

/**********************************
 * Calculate the Hamming window
 *********************************/
static void hammingWindow(double *window, int width)
{
  int i;
  for(i = 0; i < width; i++)
    window[i] = 0.54 - 0.46 * cos(2 * PI * i / (width-1));
}

/**********************************
 * Calculate the Null window
 *********************************/
static void noWindow(double *window, int width)
{
  int i;
  for(i = 0; i < width; i++)
    window[i] = 1.0;
}


/*********************************************
 * truncate the number of bits in the filter
 ********************************************/
static void truncateBits(double *kernel, int width, int bits)
{
  int i;
  for(i = 0; i < width; i++)
    if(kernel[i]<0)
      kernel[i] = (int)(kernel[i]*(1<<bits)-0.5) / (double)(1<<bits);
    else
      kernel[i] = (int)(kernel[i]*(1<<bits)+0.5) / (double)(1<<bits);
}

/****************************************************
 * Scale the filter so it totals to the required gain
 ***************************************************/
void setGain(double *kernel, int width, double gain)
{
  double total = 0.0;
  int i;

  for(i = 0; i < width; i++)
    total += kernel[i];

  for(i = 0; i < width; i++)
    kernel[i] = kernel[i]/total;

  for(i = 0; i < width; i++)
    kernel[i] *= gain;

  total = 0.0;
  for(i = 0; i < width; i++)
    total += kernel[i];

}

/**************************************/
void printKernel(double *kernel, int width) {
  int i;
  for(i = 0; i < width; i++)
    printf("%4i, %17.14f\n", i, kernel[i]);
}

/**************************************/
void calcFreqResponse(double *kernel, int width) {
  double *freq;
  int i, fft_width = 1;
  double *fft_x, *fft_y;

  /* Work out the size */
  while(fft_width < width)
      fft_width *= 2;
  if(fft_width < 128)
	  fft_width	= 128;

  freq   = (double *)malloc(sizeof(double)*fft_width);
  fft_x  = (double *)malloc(sizeof(double)*fft_width);
  fft_y  = (double *)malloc(sizeof(double)*fft_width);
  assert(freq != NULL && fft_x != NULL && fft_y != NULL);

  /* Set up the FFT data */
  for(i = 0; i < fft_width; i++) {
    fft_x[i] = (i < width) ? kernel[i] :0.0;
    fft_y[i] = 0.0;
  }

  /* Perform the DFT */
  DFT(1, fft_width, fft_x, fft_y);

  /* Combine the real and imaginary back into absolute */
  for(i = 0; i < fft_width/2; i++)
    freq[i] = sqrt(fft_x[i] * fft_x[i] + fft_y[i] * fft_y[i]);

  /* Find the maximum FFT output and scale it to be 1.0 */
  double max = 0.0;
  for(i = 0; i < fft_width/2; i++)
  {
    if(freq[i] > max)
      max = freq[i];
  }
  for(i = 0; i < fft_width/2; i++)
       freq[i] /= max;

  /* Convert to db */
  for(i = 0; i < fft_width/2; i++) {
    freq[i] = (freq[i] == 0.0) ? -99.0 : log(freq[i])/log(10.0)*20;
  }

  graphFFT(freq,fft_width);

}

void usage(void)
{
  fprintf(stderr, "Usage: %s [-w width] [-b bits] [-f cutoff] [-t type]\n");
  fprintf(stderr, "\n");
  fprintf(stderr, "   width   Width of the kernel\n");
  fprintf(stderr, "   bits    Truncate down to 'bits' accuracy in the kernel\n");
  fprintf(stderr, "   cutoff  Cutoff of the filter (0.0 < cutoff < 0.5)\n");
  fprintf(stderr, "   type    Window type: 'h' = Hamming, 'b' = Blackman, 'n' = none\n");
}
/*******************************
 * Calculate the a Windowed Sync 
 * filter and display it.
 ******************************/
int main(int argc, char *argv[])
{
  int i;
  double total = 0, gain = 1.0, freq=0.2;;

  int width = 51;
  int bits = 0;
  int current_arg = 1;

  char windowType = 'b';
  double *sinc, *window, *values;

  while(current_arg != argc)
  {
      if(argv[current_arg][0] ==  '-' || argv[current_arg][0] ==  '/')
      {
		  /* Make sure that it only has one letter in the switch */
		  if(argv[current_arg][1] != '\0' && argv[current_arg][2] != '\0' || argv[current_arg+1] == NULL)
		  {
              usage();
              return 0;
		  }

          switch(argv[current_arg][1]) {
          case 'w':
                width = atoi(argv[current_arg+1]);
                current_arg += 2;
                break;
          case 'b':
                bits = atoi(argv[current_arg+1]);
                current_arg += 2;
                break;
          case 'f':
                sscanf(argv[current_arg+1],"%lf",&freq);
                current_arg += 2;
                break;
          case 't':
                if(strcmp(argv[current_arg+1],"b") == 0)
                    windowType = 'b';
                else if(strcmp(argv[current_arg+1],"h") == 0)
                    windowType = 'h';
                else if(strcmp(argv[current_arg+1],"n") == 0)
                    windowType = 'n';
                else {
                    usage();
                    return 0;
                }
                current_arg += 2;
                break;
          default:
              usage();
              return 0;
          }
      }
      else {
          usage();
          return 0;
      }
          
  }

  printf("Windowed Sinc filter\n");
  printf("  Width              : %i\n",    width);
  printf("  Critical frequency : %5.3f\n", freq);
  if(bits > 0)
    printf("  Bits               : %i\n",  bits);
  printf("\n\n");

  /* Allocate storage */
  sinc   = (double *)malloc(sizeof(double)*width);
  window = (double *)malloc(sizeof(double)*width);
  values = (double *)malloc(sizeof(double)*width);
  assert(sinc != NULL && window != NULL && values != NULL);


  /* Generate the Sinc function */
  sincFunction(sinc, width, freq);

  /* Calculate window*/
  if(windowType == 'h')
    hammingWindow(window,width);
  else if(windowType == 'b')
    blackmanWindow(window,width);
  else
    noWindow(window,width);

  /* Apply Window */
  for(i = 0; i < width; i++)
      values[i] = window[i] * sinc[i];

  setGain(window, width, gain);

  /* Truncate the number of bits (less one for the sign!)*/
  if(bits > 0)
    truncateBits(values, width, bits-1);


  printKernel(values, width);
  graphKernel(values, width);

  calcFreqResponse(values, width);

  /* Pause for the user to press a key */
  getchar();
}

Personal tools