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
Author Message
newguy

Joined: 24 Jun 2004
Posts: 1900

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 }
 Display posts from previous: All Posts1 Day7 Days2 Weeks1 Month3 Months6 Months1 Year Oldest FirstNewest First
 All times are GMT - 6 Hours Page 1 of 1

 Jump to: Select a forum Software----------------General CCS C DiscussionCode LibraryEZ App LynxBest OfKnown Issues Hardware----------------CCS ICD / Mach X / Load-n-Go
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