I want to show the Magnitude response of three filters, as one response

조회 수: 38 (최근 30일)
Carl Eranio
Carl Eranio 2019년 12월 6일
댓글: Star Strider 2019년 12월 10일
So, I am to design three butter filters, a low pass, band pass, and a high pass, with the biggest wrinkle being that they must each have a different order (other wise, I'd just create a band stop filter with two stops). (Technically 7, 8, and 10, but I'll get to that later.) I have no problem with that, but I want to get the freq response as a single line that shows what happens when all three filters are inline with one another. (Right now I can just show how each one campares, which is mostly useless, from my perspective. Below is where I am so far. Thank you for the help.
clear;
freqSample = 10000;
freqCutoffLow = 100;
bandPassLower = 500;
bandPassUpper = 1000;
freqCutoffHigh = 2000;
%as = [1/(freqCutoff*2*pi),1];
nLow = 8; % orderOfFilter;
nBandPass = 8;
nHighPass = 10;
wNLow= freqCutoffLow / (freqSample/2); %normalizedCutOffFreq
[lowA,lowB] = butter(nLow, wNLow);
[passA,passB] = butter(nBandPass,[bandPassLower bandPassUpper]/(freqSample/2));
[highA, highB] = butter(nHighPass, freqCutoffHigh/(freqSample/2),'high');
fvt = fvtool(lowA, lowB, passA, passB, highA, highB);
legend(fvt,'butterLow', 'butterHigh', 'butterPass')

답변 (3개)

Max Murphy
Max Murphy 2019년 12월 6일
Could you use freqz and take the product of the frequency responses, assuming the filters will be cascaded in a linear system?
H_low = freqz(lowA,lowB,512);
H_pass = freqz(passA,passB,512);
H_high = freqz(highA,highB,512);
freqz will return the frequency response at a specified number of points, but since you have already determined your filter coefficients, the filter order will not be affected. Then, you can take the element-wise product of those responses.

Star Strider
Star Strider 2019년 12월 6일
It’s always worth experimenting.
It seems that convolving the numerators and denominators of the individual filters is appropriate:
freqSample = 10000;
freqCutoffLow = 100;
bandPassLower = 500;
bandPassUpper = 1000;
freqCutoffHigh = 2000;
%as = [1/(freqCutoff*2*pi),1];
nLow = 8; % orderOfFilter;
nBandPass = 8;
nHighPass = 10;
wNLow= freqCutoffLow / (freqSample/2); %normalizedCutOffFreq
[lowA,lowB] = butter(nLow, wNLow);
[passA,passB] = butter(nBandPass,[bandPassLower bandPassUpper]/(freqSample/2));
[highA, highB] = butter(nHighPass, freqCutoffHigh/(freqSample/2),'high');
Av = conv(conv(lowA,passA),highA); % Convolve
Bv = conv(conv(lowB,passB),highB); % Convolve
figure
freqz(Av, Bv, 2^14, freqSample)
title('Convolved')
figure
freqz(lowA, lowB, 2^14, freqSample)
title('Low')
figure
freqz(passA, passB, 2^14, freqSample)
title('Bandpass')
figure
freqz(highA, highB, 2^14, freqSample)
title('High')
s = zeros(1,10001);
s(5000) = 1;
LowOut = filtfilt(lowA, lowB, s); % First Filter
BandOut = filtfilt(passA, passB, LowOut); % Cascade
AllOut = filtfilt(highA, highB, BandOut); % Cascade
FTs = fft(s);
FTA = fft(AllOut);
H = FTA ./ FTs; % Transfer Function
Fv = linspace(0, 1, fix(numel(s)/2)+1)*(freqSample/2);
Iv = 1:numel(Fv);
figure
plot(Fv, 20*log10(abs(H(Iv))/max(abs(H(Iv)))))
grid
The serial output of the three filters (the last figure) appears to give essentially the same result as the convolution (the first figure).
  댓글 수: 4
Carl Eranio
Carl Eranio 2019년 12월 10일
편집: Carl Eranio 2019년 12월 10일
Thank you, that's my bad. And your answer makes sense. (I did learn alot about that from your code, and it really helps to have convolution means series in my head now.)
Star Strider
Star Strider 2019년 12월 10일
My pleasure.
No worries.
If my Answer helped you solve your problem, please Accept it!

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


Carl Eranio
Carl Eranio 2019년 12월 8일
So, my teammate worked out the answer we were looking for. It apears the answer, graphically speaking, was to only use the 'h' terms. I would like to thank the two of you for your help.
%% LAB2 EEET2447 Task 1
clear all
close all
clc; %clear
% 8th order Butterworth low-pass filter
%reference*****************************
% http://mirlab.org/jang/books/audioSignalProcessing/filterDesign.asp?title=11-2%20Filter%20Design%20(%C2o%AAi%BE%B9%B3%5D%ADp)
%**************************************
fc1 = 100; % cutoff frequency
fs = 1e4; % sampling frequency
order1 = 8;
[b,a] = butter(order1, 2*fc1/fs, 'low');
[sos,g] = tf2sos(b,a);
% === Plot frequency response
[hLow, wLow]=freqz(b, a);
%Get the frequencies in Hz
T = 1/fs;
Hertz = wLow/(2*pi*T);
plot(Hertz, abs(hLow));
title('Magnitude frequency response');
xlabel('Normalized Frequency');
xlim([0 2500]);
grid on
hold on
% 7th order band-pass butterworth filter
fmin = 500; %lower cut-off
fmax = 1000; %higher cut-off
order2 = 7;
[d,c] = butter(order2, [fmin/(fs/2), fmax/(fs/2)], 'bandpass');
[hPass, wPass]=freqz(d, c);
plot(wPass/pi*fs/2, abs(hPass));
grid on
hold on
% 10th order high-pass butterworth filter
fc3 = 2000; % cutoff frequency
order3 = 10;
[f,e] = butter(order3, 2*fc3/fs, 'high');
[sos,g] = tf2sos(f,e);
% === Plot frequency response
[hHigh, wHigh]=freqz(f, e);
plot(wHigh/pi*fs/2, abs(hHigh));
grid on
hold on
  댓글 수: 1
Max Murphy
Max Murphy 2019년 12월 8일
No problem. Looks like the filters were in parallel, apologies for the misleading answer! Star strider's comment below is the correct answer.

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

카테고리

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

제품


릴리스

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by