MATLAB: Magnitude and Phase response of a Lowpass filter

lowpassphaseresponse

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' );

Best Answer

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 frequncy domain, i.e., use the ratio of the frequency response of the output to the frequency resposne of the input to directly stimate 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 suffere as you 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
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 phase 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 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.
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
plot(t,Y(:,end),t,y/abs(h(end))),grid
Related Question