Amplitude of signal not matching up after filtering when calculating it using different methods.

Hi,

I am running a code where I am generating a signal. I am doing processing on the signal using 2 ways mentioned below:

1. I am filtering the signal and then calculating the amplitude of the signal (Sum1).
2. I am breaking the signal into chunks of 2000 data samples, filtering each chunk seperately and calculating amplitude of each chunk of data and then summing up the amplitude of all the chunks (Sum2).

I have tried this code with both [b,a] method and [z,p,k] method. The problem is I am not getting both the sums matched up. There's always a difference between them, however I want to mention that when using [z,p,k] method, the difference is very low (<2) as opposed to (~30) when using [b,a] method. I am not sure why that is happening. I am not using filtfilt since I have to implement my code in hardware and hence my hardware won't be doing processing the way as filtfilt does. One reason I can think of is maybe related to phase shift? Can someone put their thoughts into this? Thank you!

fs=44100;
[b,a]=ellip(4, 1, 40, 100/(fs/2),'high'); %Elliptical filter
% [z,p,k] = ellip(4,1,40,2000/(fs/2),'high'); %z,p,k method
% [sos,g] = zp2sos(z,p,k);

y = rand (10000,1);
y_h2=filter(b,a,y);
Sum1 = sum(abs(y_h2)); %amplitude of whole signal together 

s2_sections = buffer(y,2000);
len=length(y);
[rows, columns] = size(s2_sections);
s2_filtered = zeros(rows,columns);
for col = 1:columns
    thiscolumn = s2_sections(:,col);
    s2_filtered(:,col) = filter(b,a,thiscolumn);
end

[rows2,columns2] = size(s2_filtered);
y2_amplitude = zeros(1,columns2);
for col2=1:columns2
    thiscolumn2 = s2_filtered(:,col2);
    y2_amplitude(:,col2)=sum(abs(thiscolumn2));
end

Sum2 = sum(abs(y2_amplitude)); %amplitude of signal when processing done on chunks of data 
diff=Sum1-Sum2;

Answer 1

[Walter Roberson]

Every separate filter() call initializes the filter state to zeroes before processing the figure. Every window that you process this way effectively introduces transcients.

You can avoid this by outputing the filter state and using that in the next filter() call:

for col = 1:columns
    thiscolumn = s2_sections(:,col);
    if col == 1
      [s2_filtered(:,col), zi] = filter(b, a, thiscolumn);
    else
        [s2_filtered(:,col), zi] = filter(b, a, thiscolumn, zi);
    end
end

[Giggs B.]

Thank you @Walter Roberson. I am getting output now with difference of <1. I checked what is zi and i saw is is "initial conditions for filter delay". Did you mean that everytime filter() is called, it introduces a new delay? Also, could you explain what does this "delay" actually mean when using filter()?

[Walter Roberson]

The output of a filter for any particular step is always defined in terms of the inputs before that step. For the particular kind of filter you are using, that would look something like

y[n] = c1*x[n] + c2*x[n-1] + c3*x[n-2] (etc)

It would not be uncommon for simple filters to look something like

y[n] = 1/3*x[n] - 1/6*x[n-3] - 1/6*x[n - 7]

but filters such as butterworth might well have longer representations.

Each reference to a the input at an earlier time requires implicitly "delaying" the corresponding signal. In the above example, the signal needs to be delayed 3 samples and 7 samples to be available to be mixed in to calculate the current value.

We then run into the question of what to do when we are processing near the beginning of the signal, before we have seen enough samples to be able to delay as much as is requested. In the above example, what should be done for the first 3 samples when x[n-3] is not yet available, and what should be done for the first 7 samples when x[n-7] is not yet available?

What filter() without any extra parameters does is initialize those slots to 0. In the above example it would keep a list of the previous 7 signals, initialized to 0, and when a new sample comes it, it uses the new sample and the memory to calculate an output, then it drops the first (oldest) entry off the list and puts the newest sample at the end of the list.

If you interrupt the filtering and resume it (because you are windowing) then each time you end the filter() call, it discards the signal history, and the next call would have to start over with 0's in the buffer.

But if you ask for the second output, then you get a copy of the signal history, and if you pass that in to filter() then filter() can pick up where it left off, just as if it had not been interrupted.

[Giggs B.]

Oh, I get it now, Thank you so much for the great explanation :) :) :)