CCS C Software and Maintenance Offers
FAQFAQ   FAQForum Help   FAQOfficial CCS Support   SearchSearch  RegisterRegister 

ProfileProfile   Log in to check your private messagesLog in to check your private messages   Log inLog in 

CCS does not monitor this forum on a regular basis.

Please do not post bug Reports on this forum. Send them to support@ccsinfo.com

Digital Filter Design Implementation (dsPIC, FIR Filter)

 
Post new topic   Reply to topic    CCS Forum Index -> Code Library
View previous topic :: View next topic  
Author Message
newguy



Joined: 24 Jun 2004
Posts: 1766

View user's profile Send private message

Digital Filter Design Implementation (dsPIC, FIR Filter)
PostPosted: Tue Sep 14, 2021 11:26 am     Reply with quote

If you've never tackled DSP filter design, pretty much everything you need to know can be found here. The filter below is a FIR (finite impulse response) type, as in my application I needed to have a very well defined group delay.

Filter Design
Download and install gnu Octave (https://www.gnu.org/software/octave/download). It's a free/open source alternative to Matlab and is generally fantastic. The code below can be downloaded and run under Octave. Name the file whatever you please with a .m extension and Octave will run it directly. It is based upon the article by Tim Youngblood (https://www.allaboutcircuits.com/technical-articles/design-of-fir-filters-design-octave-matlab/)

Code:
% save this as "something".m
% it has been developed with gnu Octave (freeware)
% this is based on Tim Youngblood's All About Circuits article:
% https://www.allaboutcircuits.com/technical-articles/design-of-fir-filters-design-octave-matlab/
close all;
clear all;
clf;
pkg load signal;

% low pass filter characteristics
f1 = 500; # corner frequency
f2 = 575; # freq at which the attenuation below is in effect
delta_f = f2 - f1;
Fs = 5000; # sampling frequency
dB = 8; # attenuation in dB
N = dB * Fs / (22 * delta_f); # filter order

gw = gausswin(round(N)); # gaussian window to be applied to the filter coefficients

f = [f1] / (Fs/2);
hc = fir1(round(N) - 1, f, 'low', gw); % this generates the DSP filter coefficients and they are a floating point
% number. The sum of all the coefficients adds up to 1.0000

step = ones(1, size(hc, 2)); % this is required to see the step response of the filter later if things like
% over/under shoot are important to you

% the dsPIC requires filter coefficients which are integers (but it will treat them as signed Q15)
% simply take the coefficients in hc above and multiply them by 32768. Be careful as they must be integers
% and their sum must total 32768 otherwise you're going to be introducing some gain/attenuation to your filter
h_q1 = [32, 21, -44, -185, -370, -474, -273, 465, 1830, 3620, 5347, 6415, 6415, 5347, 3620, 1830, 465, -273, -474, -370, -185, -44, 21, 32];

plot((-0.5:1/4096:0.5-1/4096)*Fs, 20 * log10(abs(fftshift(fft(hc, 4096))))) % this plots the original filter
% coefficients (floating point)
hold on;
plot((-0.5:1/4096:0.5-1/4096)*Fs, 20 * log10(abs(fftshift(fft(h_q1 / 32768, 4096))))) % this plots our
% quantized coefficients (integer), and is useful to ensure that we didn't drastically alter the filter's
% response when we did quantize it
axis([0 Fs/2 -60 5])
title('Filter Frequency Response')
grid on

x = 10 * sin(50*pi*[1:1000]*2/Fs) + 12 * sin(60*pi*[1:1000]*2/Fs) + 2 * sin(300*pi*[1:1000]*2/Fs) + sin(500*pi*[1:1000]*2/Fs) + 20 * sin(1000*pi*[1:1000]*2/Fs) + 30 * sin(2000*pi*[1:1000]*2/Fs);
% generate an input waveform to our filter with signals of interest and a bunch of high power high frequency noise
sig = 20 * log10(abs(fftshift(fft(x, 4096))));
%xf = filter(hc, 1, x); % this uses the original float version of the filter
xf = filter(h_q1 / 32768, 1, x); % this uses our quantized filter
figure
subplot(211);
plot(x);
title('Sinusoid with multiple freq components');
subplot(212);
plot(xf);
title('Filtered Signal');
xlabel('time')
ylabel('amplitude')

x = x / 512;
sig = 20 * log10(abs(fftshift(fft(x, 4096))));
%xf = filter(hc, 1, x); % original float filter
xf = filter(h_q1 / 32768, 1, x); % quantized filter

figure
subplot(211)
plot((-0.5:1/4096:0.5-1/4096) * Fs, sig);
hold on
%plot((-0.5:1/4096:0.5-1/4096) * Fs, 20 * log10(abs(fftshift(fft(hc, 4096)))),'color','r');
plot((-0.5:1/4096:0.5-1/4096) * Fs, 20 * log10(abs(fftshift(fft(h_q1 / 32768, 4096)))),'color','r');
hold off
axis([0 Fs/2 -60 30]);
title('Input to filter');
grid on
subplot(212)
plot((-0.5:1/4096:0.5-1/4096) * Fs, 20 * log10(abs(fftshift(fft(xf, 4096)))));
axis([0 Fs/2 -60 25]);
title('Filter output');
xlabel('Hz');
ylabel('dB');
grid on;

figure
plot(xcorr(hc, step)); % step response of the filter in the time domain
grid on;
title('Step Response');
axis([0 size(hc, 2) -.2, 1.2]);


The design is for a FIR LPF with about a 300Hz cutoff. You can experiment with the various parameters to fine tune the response to your liking. I added the gaussian window to the filter because it helped reduce the over/undershoot in the filter's step response.

dsPIC Implementation
The demo program below is a very stripped down implementation of the real-time filter. It simulates an analog waveform being acquired point-by-point and subsequently loaded into a circular buffer. The circular buffer's size needs to be the same size as the filter and in this case is 24 elements.

As each new analog reading (simulated in this case) arrives, it is stuffed into the next available slot in the circular buffer. In this example, the first sample is placed in input_wave[0], then the next sample is placed in input_wave[1], and so on up to [23]. At that point, filling wraps back to [0] and increases from there, and so on.

The filtering algorithm multiplies the most recent analog sample (input_wave[]) with filter coefficient [0]. The next most recent sample is multiplied with filter coefficient [1], and so on for all 24 elements of both the input_wave[] and dsp_coeffs[] arrays. This math is handled automatically by the dsPIC's modulo addressing capability which automatically handles properly wrapping when either the max or min array element is encountered. Once all 24 input_wave[] samples have been multiplied with all 24 dsp_coeffs[], the result of all those multiplications is returned. The example program doesn't currently do anything with the filtered result, but it should be quite obvious how you would take advantage of this in your own program.

The actual filter function takes just 51 instruction cycles. At 70MIPS this equates to about 729ns. The filter order (i.e. number of filter coefficients) will influence the execution time, but in general it's very fast and more than sufficient for any "real time" filtering you may need, including filtering of more than one A/D channel simultaneously.

Code:
#include <33EP512GP806.h>
#use delay(clock=140mhz)

#word CORCON = getenv("SFR:CORCON")
#word MODCON = getenv("SFR:MODCON")
#word XMODSRT = getenv("SFR:XMODSRT")
#word XMODEND = getenv("SFR:XMODEND")
#word YMODSRT = getenv("SFR:YMODSRT")
#word YMODEND = getenv("SFR:YMODEND")

#define FILTER_SIZE 24
#BANKY
signed int16 input_wave[FILTER_SIZE];
#BANKX
signed int16 dsp_coeffs[FILTER_SIZE];

signed int16 apply_fir_filter(signed int16 *waveform, signed int16 *filter_coeffs);

void main(void) {
   unsigned int16 j = 0;
   signed int result;
   signed int16 *ptr;
   
   // initialize the input waveform to 0
   for (j = 0; j < FILTER_SIZE; j++) {
      input_wave[j] = 0;
   }
   // initialize the DSP coeffs
   dsp_coeffs[0] = 32;
   dsp_coeffs[23] = 32;
   dsp_coeffs[1] = 21;
   dsp_coeffs[22] = 21;
   dsp_coeffs[2] = -44;
   dsp_coeffs[21] = -44;
   dsp_coeffs[3] = -185;
   dsp_coeffs[20] = -185;
   dsp_coeffs[4] = -370;
   dsp_coeffs[19] = -370;
   dsp_coeffs[5] = -474;
   dsp_coeffs[18] = -474;
   dsp_coeffs[6] = -273;
   dsp_coeffs[17] = -273;
   dsp_coeffs[7] = 465;
   dsp_coeffs[16] = 465;
   dsp_coeffs[8] = 1830;
   dsp_coeffs[15] = 1830;
   dsp_coeffs[9] = 3620;
   dsp_coeffs[14] = 3620;
   dsp_coeffs[10] = 5347;
   dsp_coeffs[13] = 5347;
   dsp_coeffs[11] = 6415;
   dsp_coeffs[12] = 6415;

   YMODEND = &input_wave[0];
   YMODEND += 2 * (FILTER_SIZE - 1); // end must be odd: the end address is automatically made to be odd here
   YMODSRT = &input_wave[0]; // Set the End address
   
   XMODEND = &dsp_coeffs[0];
   XMODEND += 2 * (FILTER_SIZE - 1); // end must be odd: the end address is automatically made to be odd here
   XMODSRT = &dsp_coeffs[0];
   ptr = &input_wave[0];

   j = 0; // reset our index
 
   while (TRUE) {
      // insert new dummy sample into input_wave[] (this simulates a new A/D reading being acquired)
      *ptr = 0x7ff0;
      result = apply_fir_filter(ptr, &dsp_coeffs[0]); // filter it (you'll have to store the result obviously)
      if (++j == FILTER_SIZE) { // increment our index
         ptr = &input_wave[0]; // wrap if necessary
         j = 0;
      }
      else {
         *ptr++;
      }
   }
}

signed int16 apply_fir_filter(signed int16 *addr_waveform, signed int16 *addr_filter_coeffs) {
   signed int16 result;

   // waveform is filled in the traditional manner for a circular buffer - in ascending order; buffer wraps back to
   // index 0 after highest index has been filled
   // most recent sample is at the current buffer index, and oldest is at the "next" index, be that one higher than the current
   // or 0 if the current index is the highest
   // algorithm below takes advantage of the dsPIC's built in modulus functionality whereby buffer wrapping is handled automatically
   // whether you're stepping through the buffer in ascending or descending order so we step through the filter coeffs in ascending order
   // and multiply each of them with our circular buffer of waveform samples in descending order (most recent - oldest)

   #ASM
   PUSH W0
   PUSH W1
   PUSH W4
   PUSH W5
   PUSH W8
   PUSH W10
   PUSH CORCON   
   PUSH MODCON
   MOV #0x00b0,W0
   MOV W0,CORCON
   // fractional mode enabled
   // unbiased (convergent) rounding enabled
   // stack frame not active
   // CPU priority 7 or less
   // accumulator = 9.31 (super saturation)
   // data space write saturation is enabled
   // accumulator B saturation disabled
   // accumulator A saturation enabled
   // DO loops are active
   // DSP multiplies are signed
   // fixed exception processing
   MOV #0xcfa8,W0
   MOV W0,MODCON
   // bank X uses W8
   // bank Y uses W10
   // no bit-reversed addressing
   // modulo addressing enabled for Y bank
   // modulo addressing enabled for X bank
   MOV addr_waveform,W10  // store address of waveform in W10 (bank Y)
   MOV addr_filter_coeffs,W8 // store address of filter coeffs in W8 (bank X)
   MOV #22, W1 // filter size - 2
   CLR   A,[W8]+=2,W4,[W10]-=2,W5 // clear accumulator and preload first waveform sample and first filter coeff

   REPEAT W1
   MAC   W4*W5,A,[W8]+=2,W4,[W10]-=2,W5
   MAC   W4*W5,A

   SAC.R   A,W0 // Store result in W0
      MOV W0, result  // Save the result
   POP MODCON
   POP CORCON
   POP W10
   POP W8
   POP W5
   POP W4
   POP W1
   POP W0
   #ENDASM

   return (result); // return the result
}
Display posts from previous:   
Post new topic   Reply to topic    CCS Forum Index -> Code Library All times are GMT - 6 Hours
Page 1 of 1

 
Jump to:  
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum


Powered by phpBB © 2001, 2005 phpBB Group