Magnitude and Phase response of a Lowpass filter
65 次查看(过去 30 天)
显示 更早的评论
I have a problem with understanding the phase response of lowpass filter in MATLAB(I'm writing my own code not using inbuilt functions to find phase response & Matlab). I am trying to pass sine signals of different frequencies into a lowpass filter with a certain passband frequency. Later, magnitude response is obtained by the change in the output amplitude divided by input amplitude . Whereas phase response is the angle of (output amplitude/input amplitude)???. I am getting correct magnitude response, but my phase response is zero. I don't understand why ??
I mean even though there isnt a phase term in the test signal ,there should be some phase change due to the lowpass filter at least. I almost spent more than two weeks on this and still unable to understand .
Im attaching the graphs for the amplitude response and phase response ,along with the code.
fs=1000;
f=1:15;
t = 0:1/fs:10-1/fs;
x= sin(2*pi*t'*f);
[Y,d]= lowpass(x,10,fs);
figure;
plot(t,x);
hold on
plot(t,Y);
for i=1:15
figure(i);
plot(t,x(:,i));
hold on
plot(t,Y(:,i));
title ('Original/ Filtered Data')
%%%%%%%%% peak to peak amp %%%%%%%%%%%%
input(i)=(rms( x(:,i))/0.707)*2;%%%% inp peak to peak
output(i)=(rms(Y(:,i))/0.707)*2 ;%%%%%% out peak to peak
change(i)=output(i)./input(i);%%%%%% change in out peak to peak to input
phase(i)=angle(change(i));
end
figure;
subplot(2,1,1);
plot(f,mag2db(abs(change)));
title('Frequency response of the filter using abs(Out/Inp)')
xlabel('Frequency');
ylabel ('Change in output/input P-P' );
subplot(2,1,2);
plot(f,phase);
title('Phase response of the filter using angle(Out/Inp)')
xlabel('Frequency');
ylabel ('phase' );
2 个评论
Paul
2020-9-3
Are you trying to understand why the lowpass function in the Signal Processing Toolbox does not introduce a phase shift in the output?
Or are you trying to develop a procedure to determine the frequency response of an unknown system by injecting a series of sine waves of various frequencies?
采纳的回答
Paul
2020-9-4
编辑:Paul
2023-9-7
it sounds like there are two issues here: 1. estimating a frequency response using test inputs, 2. what does lowpass do. Let's tackle them in sequence.
Estimating the frequency response: You were on the right track with injecting sine waves of different frequencies and comparing the input to the output. A couple of things to keep in mind. You need to make sure to either delete the initial transient from the output (and the corresponding input time window in the input) or ensure that the initial transient is a small proportion of the overall response. Keep in mind that the response you're trying to estimate is for a particular frequency. Taking the ratio of the RMS's to estimate the magnitude response might be o.k. for linear systems, but in general you have to be careful that you're not RMS'ing components of the output at frequencies other than the input test frequency.
To get the phase response in the time domain you need to estimate the delay between the input and the output at the test frequency and then convert to phase. For sure, taking the angle of the RMS ratio will not yield that delay. If you really want to estimate the delay in the time domain, there is a function called finddelay that may be of use.
However, the alternative approach is to work in the frequency domain, i.e., use the ratio of the frequency response of the output to the frequency resposne of the input to directly estimate the magnitude and phase response of the system. The code that follows shows how to do that for a simple low pass filter. Keep in mind that any approach you use may begin to suffer as your test frequency gets close to the Nyquist frequency. You can experiment with this code to see how close you can get before this simple approach begins to break down.
% filter parameters
fs = 250; % Hz
Ts = 1/fs;
[b,a] = butter(3,0.25); % butterworth filter
% frequency response for plotting, Hz
fplot = logspace(-2,pi,1000)/pi*fs/2;
hplot = freqz(b,a,fplot,fs);
% test signal input, with lots of cycles and time
f = 30; % Hz
t = 0:Ts:(100/f); t = t(:);
u = sin(2*pi*f*t);
% filtered output
y = filter(b,a,u);
% FFT's, interpolated to the test frquency
f_fft = (0:length(t)-1)/(length(t)-1)*fs;
Y_fft = interp1(f_fft,fft(y),f);
U_fft = interp1(f_fft,fft(u),f);
% true filter response at the test frequency
h = freqz(b,a,[0 f],fs); h = h(2);
% compare true gain to FFT ratio
m_comp = 20*log10([abs(h) abs(Y_fft/U_fft) (rms(y)/rms(u))])
% compare true phase to FFT ratio
p_comp = 180/pi*angle([h Y_fft/U_fft])
% make a plot
figure
subplot(211);
semilogx(fplot,db(abs(hplot)),f,m_comp(2),'o'),grid,set(gca,'xlim',[10 100]);
subplot(212);
semilogx(fplot,180/pi*angle(hplot),f,p_comp(2),'o'),grid,set(gca,'xlim',[10 100]);
As you discovered, the low pass filter returned from the lowpass function you called is indeed a linear phasek,FIR filter. However, the lowpass function applies that filter in such a way to yield something very close to zero-phase filtering, which is what I think the doc means by "compensates for the delay introduced by the filter," which is not a very clear statement and is not explained in the doc. The filtfilt function does something similar, and is explained in doc filtfilt. Let's compare the output of lowpass for a particular input to the output of filtfilt using the FIR filter returned by lowpass. Note that the effective gain of filtfilt is the magnitude squared of the filter, so the output of filtfilt has to be divided by the filter gain to match the result from lowpass. If you zoom in on the plot, you'll see that lowpass and filtfilt must use different approaches near the intial and final times of the response for a FIR filter. I believe that lowpass does a simpe shift for a FIR filter and makes call to filtfilt for an IIR filter.
fs = 1000;
f = 60;
t = 0:1/fs:10-1/fs;
x = sin(2*pi*t'*f);
[Y,d]= lowpass(x,10,fs);
y = filtfilt(d.Coefficients,1,x(:,end));
h = freqz(d,[0 f],fs); h = h(2); % need to specify at least two frequencies, see doc freqz
figure
plot(t,Y(:,end),t,y/abs(h(end))),grid
copyobj(gca,figure)
xlim([0 1])
2 个评论
Paul
2020-9-12
There should be button to click to accept an answer. Don’t know exactly where, but it should be straightforward.
更多回答(1 个)
Robert U
2020-9-3
Hi Suvvi Kuppur Narayana Swamy:
The function angle() is calculating the phase angle of a complex number. Since all values are real there would always be a phase angle of zero. The function lowpass() compensates for delays according to documentation. If you would calculate the phase angle between the input and the output function it should remain close to zero anyway.
Example without using signal toolbox including calculations on transfer function:
Kind regards,
Robert
2 个评论
Robert U
2020-9-3
Hi Suvvi Kuppur Narayana Swamy:
I assume that actually you are interested in analogue system transfer functions.
f=1:15;
% second order filter with 10 Hz pass band
H = @(f,s) 1./(1./(2*pi*f).^2*s.^2 + 2./(2*pi*f)*s + 1);
H10 = @(s) H(10,s);
% magnitude estimation
magn = abs(H10(1j*2*pi*f));
% phase estimation
phase = rad2deg( angle(H10(1j*2*pi*f)) );
figure;
subplot(2,1,1);
semilogx(f,mag2db(magn));
title('Frequency response of the filter using abs(Out/Inp)')
xlabel('Frequency');
ylabel ('Change in output/input P-P' );
subplot(2,1,2);
semilogx(f,phase);
title('Phase response of the filter using angle(Out/Inp)')
xlabel('Frequency');
ylabel ('phase' );
% check against inbuilt functions
H_tf = @(f) tf(1,[1./(2*pi*f).^2 2./(2*pi*f) 1]);
figure;
bodeplot(H_tf(10),2*pi*f)
If you want to calculate the time discrete transfer functions you would have to convert the transfer functions from Laplace-domain into z-Domain. Inbuilt is a function c2d().
Kind regards,
Robert
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!