Amplitude of signal not matching up after filtering when calculating it using different methods.
6 次查看(过去 30 天)
显示 更早的评论
Hi,
I am running a code where I am generating a signal. I am doing processing on the signal using 2 ways mentioned below:
- I am filtering the signal and then calculating the amplitude of the signal (Sum1).
- 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;
0 个评论
采纳的回答
Walter Roberson
2021-8-5
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
3 个评论
Walter Roberson
2021-8-5
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.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Digital Filter Design 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!