newguy
Joined: 24 Jun 2004
Posts: 1908
|
| Digital Filter Design Implementation (dsPIC, FIR Filter) |
 Posted: Tue Sep 14, 2021 11:26 am |
 |
|
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
} |
|
|
FAQ 
Search
Register
Profile 
Log in to check your private messages 
Log in
