Same filters but using different methods give different results

조회 수: 3 (최근 30일)
Giggs B.
Giggs B. 2022년 3월 24일
댓글: Star Strider 2022년 3월 24일
Hi,
I am getting different results when using the [b,a] filter method and [z,p,k] method. My code goes like this:
for loop for cut-off frequecy (starts from 50 Hz and goes upto 5000 Hz with step size of 50 Hz) ---------------------
Read couple of csv files in a loop (1 to 7) -> Filtering out each of them using a Elliptical band/high pass 4th order filter for "x" cut-off frequency -> Calculating amplitude (which means adding the absolute values of this filtered data) of each file -> Then taking average of the amplitude of these 1-6 files -> dividing the amplitude of 7th file by the average amplitude of 1-6 files -> store this value -> loop starts again for cut-off frequency of "x+50" upto 5000 Hz
------------------------ end loop
Now, when I plot the graph of the values using [b,a] method and by using [z,p,k] method, I can see whole different value. Why would that be the case?
fl=0;
fh=9000;
fs=40000;
array_aW=zeros(100,1);
y1=readmatrix('1.csv');
y2=readmatrix('2.csv');
y3=readmatrix('3.csv');
y4=readmatrix('4.csv');
y5=readmatrix('5.csv');
y6=readmatrix('6.csv');
y7=readmatrix('7.csv');
for e=1:100
fl=fl+50;
% [b,a]=ellip(2, 1, 40, [fl,fh]/(fs/2),'bandpass');
[z,p,k]=ellip(2, 1, 40, [fl,fh]/(fs/2),'bandpass');
[sos,g]=zp2sos(z,p,k);
% y_n1=filter(b,a,y1);
y_n1=filtfilt(sos,g,y1);
y_h1 = normalize(y_n1, 'range', [-1 1]);
amp1=sum(abs(y_h1));
% y_n3=filter(b,a,y2);
y_n2=filtfilt(sos,g,y2);
y_h2 = normalize(y_n2, 'range', [-1 1]);
amp2=sum(abs(y_h2));
y_n3=filtfilt(sos,g,y3);
% y_n3=filter(b,a,y3);
y_h3 = normalize(y_n3, 'range', [-1 1]);
amp3=sum(abs(y_h3));
y_n4=filtfilt(sos,g,y4);
% y_n4=filter(b,a,y4);
y_h4 = normalize(y_n4, 'range', [-1 1]);
amp4=sum(abs(y_h4));
y_n5=filtfilt(sos,g,y5);
% y_n5=filter(b,a,y5);
y_h5 = normalize(y_n5, 'range', [-1 1]);
amp5=sum(abs(y_h5));
y_n6=filtfilt(sos,g,y6);
% y_n6=filter(b,a,y6);
y_h6 = normalize(y_n6, 'range', [-1 1]);
amp6=sum(abs(y_h6));
%calculations
ampNoiseWeightedAvg=((0.033*amp1)+(0.033*amp2)+(0.033*amp3)+(0.05*amp4)+(0.05*amp5)+(0.8*amp6))/6;
%%%%%%%%%%%%7 file%%%%%%%%%%%%%
y_n7=filtfilt(sos,g,y7);
% y_n7=filter(b,a,y7);
y_h7 = normalize(y_n7, 'range', [-1 1]);
ampW=sum(abs(y_h7));
%%%%%%Figure of Merit%%%%%%
FigureOfMerit_ampW= ampW/ampNoiseWeightedAvg;
array_aW(e,1)=FigureOfMerit_ampW;
end
Thanks,
Gagan
  댓글 수: 2
Jan
Jan 2022년 3월 24일
"My code goes like this" - do not paraphrase, what you code does, but post the code. If it contains a bug, the paraphrasation will conceal it.

댓글을 달려면 로그인하십시오.

채택된 답변

Star Strider
Star Strider 2022년 3월 24일
I am getting different results when using the [b,a] filter method and [z,p,k] method.’
I am not surprised. Filters using transfer function ‘[b,a]’ notation can be unstable. Ues the freqz function to explore the filter characteristics:
fl=0;
fh=9000;
fs=40000;
[b,a]=ellip(2, 1, 40, eps+[fl,fh]/(fs/2),'bandpass');
figure
freqz(b, a, 2^16, fs)
[z,p,k]=ellip(2, 1, 40, eps+[fl,fh]/(fs/2),'bandpass');
[sos,g]=zp2sos(z,p,k);
figure
freqz(sos, 2^16, fs)
The phase characteristics are significantly different.
This is the optimal approach for the best filter implementation:
[z,p,k]=ellip(2, 1, 40, [fl,fh]/(fs/2),'bandpass');
[sos,g]=zp2sos(z,p,k);
Note that ‘fl=0’ designs a lowpass filter. It would likely be more efficient to design the filter as a lowpass filter rather than a bandpass filter.
Also if ‘y1’ through ‘y7’ are equal-length column vectors (or can be made to be equal length column vectors), concatenate them in a matrix and call filtfilt once with your filter and the matrix. The ooutput will be a matrix of the filtered signals.
.
  댓글 수: 2
Giggs B.
Giggs B. 2022년 3월 24일
@Star Strider, thank you, I understand now. Yes, all vectors are of same length. I will try your suggested way.
Star Strider
Star Strider 2022년 3월 24일
As always, my pleasure!

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Filter Design에 대해 자세히 알아보기

태그

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by