MATLAB - Digital Filtering in Second Order Sections - Scale Values

Hi Everyone,
I am struggling to manually apply the scale values to a digital filter I
designed in MATLAB using fdatool.

I generated the filter (IIR band stop filter consisting of four second
order sections) and exported it as an object.  I can use MATLAB's "filter"
function to filter my data correctly.

As a precursor to implementing the filter in C and embedding it onto a
hardware device, I wanted to manually write the code in MATLAB.

I can successfully apply the coefficients and the waveform shape matches
that obtained using MATLAB's "filter" function however the values have
changed because I am not applying any scaling.

I have a scale vector consisting of 5 values (this is as I expected as I
have four second order sections in my filter).

My question is how to correctly apply the scale values to my filter?

My (obviously incorrect) belief was that you apply the first scale value
to the input and the subsequent scale values to the outputs of each second
order section.

However when I do this, the shape of the data is changed and it is centred
around 0, as though I applied a high pass filter.  To get the correct
waveform I need to ignore the scale values provided by MATLAB and apply
only the coefficients and then "fudge" the output by applying a scale
value at the end which I came up with by trial and error.

Please see my code below which attempts to apply MATLAB's scale values as
per my previous "understanding":

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

  % 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];


  %% Apply the filter to the data called "serial"
  for i=10:size(serial,1)
    
    % Section 1
    output_section_1(i)=Scale(1)*(Section1(1)*serial(i) +
Section1(2)*serial(i-1) + Section1(3)*serial(i-2) -
Section1(5)*output_section_1(i-1) - Section1(6)*output_section_1(i-2));
    
    % Section 2
    output_section_2(i)=Scale(2)*(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(3)*(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(4)*(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 advice or input would be very much appreciated.

Kind Regards,
David Graham.
---------------------------------------
Posted through http://www.DSPRelated.com

The scale factor should only be applied to the numerator operations, so change this 

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


To this

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

Just moved one parentheses!

Alternatively you can just pre-multiply all the numerator coefficients times the scale factor for each filter, especially in this case where the scale factors are close to 1.0 (in many cases the scale factors could be very small in which case you might want to apply them separately to avoid coefficient quantization issues). 

Bob

>Just moved one parentheses!
>
>Alternatively you can just pre-multiply all the numerator coefficients
>times the scale factor for each filter, especially in this case where
the
>scale factors are close to 1.0 (in many cases the scale factors could be
very
>small in which case you might want to apply them separately to avoid
>coefficient quantization issues). 
>

Hi Bob,
Thank you very much for your helpful response!

I have taken your advice and now applying the scale values (as per your
description) no longer changes the waveform shape which is a great step
forward.

However, I think I am still missing something because the scale still
isn't right.  As you can see from my scale values, they are all close to 1
and so multiplying by these doesn't have a massive effect on the output.

With no scale applied, my average output is around 5,000 counts.  With
scale applied it is around 4,500 counts but my target is 1875 counts.

I can fix this by multiplying my output by 0.415 (implemented by changing
my last scale value) but this seems a "fudge" because MATLAB's filter
function just uses the filter object (where the output scale value is 1)
and gives me an average output of 1875 counts.

Thanks again for your response, I feel a step closer to getting it right
now.

Kind Regards,
David Graham.
---------------------------------------
Posted through http://www.DSPRelated.com

Where did you get your "target" value from? If you are just looking at the peaks of a complicated waveform, the action of the filter will change the peak values so this is not a reliable way of determining whether or not you have the proper gain factors. 
In some cases if you specify large amounts of passband ripple, the gain factors may be set so that the peak value of the frequency response is 0 dB, which may change the "nominal" passband gain. 
You can also check the gain at dc directly from the coefficients. Just calculate the sum of numerator coefficients divided by the sum of denominator coefficients. 


Bob

>Where did you get your "target" value from? If you are just looking at
the
>peaks of a complicated waveform, the action of the filter will change
the
>peak values so this is not a reliable way of determining whether or not
you
>have the proper gain factors. 
>In some cases if you specify large amounts of passband ripple, the gain
>factors may be set so that the peak value of the frequency response is 0
dB,
>which may change the "nominal" passband gain. 
>You can also check the gain at dc directly from the coefficients. Just
>calculate the sum of numerator coefficients divided by the sum of
>denominator coefficients. 
>
>
>Bob

Hi Bob,
Thanks again for getting in contact.

My "target" values come from using MATLAB's "filter" function using the
same filter coefficients and gain values.

So I designed the filter using fdatool and filtered it using MATLAB as
below

output=filter(Hd,input);

I then implemented the filter manually, not using MATLAB's "filter"
function and I am comparing the outputs.

So I hardcoded the coefficients and the scale values and I am attempting
to get my output to match that of MATLAB's "filter" function.

I hope that makes sense.

I know the filter will change the output, where I am struggling is in
getting my values to match those of MATLAB's filter function.  Currently
the waveform shape matches so I know I am correctly applying the
coefficients but the scale does not match.

Kind Regards,
David.
---------------------------------------
Posted through http://www.DSPRelated.com

I suspect you didn't copy enough digits into your coefficient array. The filter function is getting full precision. Try setting "format long" and then copy your coefficients with full precision. 

Bob

I suspect you didn't copy enough digits into your coefficient array. The filter function is getting full precision. Try setting "format long" and then copy your coefficients with full precision. 

Bob

>I suspect you didn't copy enough digits into your coefficient array. The
>filter function is getting full precision. Try setting "format long" and
>then copy your coefficients with full precision. 
>
>Bob

Hi Bob,
Thanks again, you are quite correct.  Once I changed the format and
included the full coefficients and scale values, the outputs matched.

Sincere thanks for your help!

David.
---------------------------------------
Posted through http://www.DSPRelated.com