200 views (last 30 days)

Show older comments

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

Paul
on 4 Sep 2020

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

Paul
on 12 Sep 2020

Robert U
on 3 Sep 2020

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

Robert U
on 3 Sep 2020

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

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

Start Hunting!