SOS Filter - How To Manually Apply Scale Values?

Hi, I used fdatool to design a bandstop filter. I exported this as an object and can successfully filter my data using it.

As a pre-cursor to implementing this filter in C, I wanted to manually apply the filter using my own code, rather than using MATLAB's filter function.

I have found that I can successfully apply the coefficients (without using the scale values) and I get the correct waveform shape, but obviously the scale has changed. I can live with this, but for the sake of correctness I want to apply the scale values so that the both waveforms (MATLAB's filter function and my own implementation) match. I can then run an automated script to detect any differences between MATLAB's filter function, my implementation and ultimately my embedded C implementation.

However I am struggling to apply the scale values manually. My understanding was that you apply scale value 1 to the input signal and then apply the subsequent scale values to each output section. I have 4 Second Order Sections in total.

The result is that the data is centered around 0 and the waveform has changed shape, as though I applied a high pass filter which removes the DC. But MATLAB's filter function is centered around 1500 (as per my input data).

The relevant sections of my code is below...

"serial" is my input data. You can see that I multiply my input data ("serial") by Scale(1) and each output by Scale(2-5)

% Bandstop Filter Coefficients
Section1 = [1 -1.9996 1 1 -1.9579 0.96158];
Section2 = [1 -1.9996 1 1 -1.9959 0.99592];
Section3 = [1 -1.9996 1 1 -1.9868 0.98685];
Section4 = [1 -1.9996 1 1 -1.9101 0.91275];

% Bandstop Scale Values
Scale = [0.98896 0.98896 0.97448 0.97448 1];

% Apply Bandstop filter
% Filter according to own code
for i=10:size(serial,1)

    % Section 1
    output_section_1(i)=Scale(2)*(Section1(1)*serial(i)*Scale(1) + Section1(2)*serial(i-1)*Scale(1) + Section1(3)*serial(i-2)*Scale(1) - Section1(5)*output_section_1(i-1) - Section1(6)*output_section_1(i-2));

    % Section 2
    output_section_2(i)=Scale(3)*(Section2(1)*output_section_1(i) + Section2(2)*output_section_1(i-1) + Section2(3)*output_section_1(i-2) - Section2(5)*output_section_2(i-1) - Section2(6)*output_section_2(i-2));

    % Section 3
    output_section_3(i)=Scale(4)*(Section3(1)*output_section_2(i) + Section3(2)*output_section_2(i-1) + Section3(3)*output_section_2(i-2) - Section3(5)*output_section_3(i-1) - Section3(6)*output_section_3(i-2));

    % Section 4
    output_section_4(i)=Scale(5)*(Section4(1)*output_section_3(i) + Section4(2)*output_section_3(i-1) + Section4(3)*output_section_3(i-2) - Section4(5)*output_section_4(i-1) - Section4(6)*output_section_4(i-2));
end

Any comments or advice would be very much apprciated.

Kind Regards, David Graham.

[David Graham]

Hi everyone, In case anyone is interested, I have an answer to my question.

A very helpful person on dsprelated.com helped answer it. There were two problems:

1. The scale values should be applied only to the feed forward values and not to the feed back sections so I needed to move my closing brackets.

2. MATLAB's filter function uses the full resolution of the coefficients whereas I was using rounded off values. By changing the format to long and using the full resolution of the coefficients, the output of my manual implementation matched MATLAT's filter function.

Regards, David.

[Brennan Kilty]

Thank you for sharing David. This has been helpful to understanding how to apply these filters outside of Matlab. Do you or anyone else know why these scalar values are extracted from the SOS filter b coefficients in the first place?

[Jason Debono]

Thank you for having posted the solution to your own query David Graham. I found it useful.

[Bob Farrell]

I will also add to the chorus of thank-you's for the answer to this question. I found myself here by way of a google search to figure out how to apply the gain factors. I could not find anything in the Matlab documentation that explained (or showed by example) what to do with the gains for second-order-sections.

I am also curious about @Brennan Kitty's additonal question: Can anybody explain why these scalar gain values are extracted from the SOS b coefficients?

I also have the following questions:

1. Why is the final value in the gain vector always = 1, when the overall filter clearly does not have unity gain?
2. I'd think that you when you export the filter values from the Filter Design Tool to your workspace when exporting as coefficients (which includes the SOS matrix and the G array (default names)), that you could feed these straight in to the Filter Visualization Tool (fvtool). But while fvtool will accept an SOS matrix, it will not accept an accompanying G array. Why not?

Regarding question #2, this seems like an omission to me. It took me a long while to understand why I would see the magnitude response plots of the Filter Design Tool and fvtool be offset from each other by 20*log10(1/prod(G(1:end-1))) dB.

Answer 1

[Diarmaid Cualain]

@Brennan Kilty, I had the same issue and I believe I figured it out. I used Matlabs filterbuilder() GUI to create a filter struct (Hd). I tested it out with a vector of my data with:

y = filter(Hd, xData);

Happy with the results, I used the Code generation tab to save the code:

% Design filter
N     = 3;     % Order
Fpass = 0.05;  % Passband Frequency
Apass = 1;     % Passband Ripple (dB)
h = fdesign.lowpass('n,fp,ap', N, Fpass, Apass);
Hd = design(h, 'cheby1');

To create a C implementation, or a real-time implementation in Matlab, I needed the a and b (feedforward and feedback) coefficients.

% Get feedforward and feedbackward coefficients
[b1,a1] = sos2tf(Hd.sosMatrix);

Unfortunately, they are not scaled correctly and will cause a large gain in your data. To fix this, I got the product of all the scale values in the filter:

% Multiply them all scaling vales against each other
productScaleValue = prod(Hd.scaleValues);

And then scaled the feedforward coefficients

% Calculate filter feedforward coefficients with no gain 
b = b1 .* productScaleValue;

This resulted in values of:

1.0e-03 * [0.2206    0.6618    0.6618    0.2206]

Which you can confirm, with the older Matlab function that creates properly scaled values:

[b,a] = cheby1(N,Apass,Fpass)

b = 1.0e-03 * [0.2206    0.6618    0.6618    0.2206]

Note: the a values remain the same.