Magnitude and Phase response of a Lowpass filter

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]

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?

[Suvvi Kuppur Narayana Swamy]

Hi Paul,

Honestly ,when i started my aim was to determine the freqeuncy response of an unknown system by pasing series of test signals of varying freqeuncies using my own code and compare that with MATLAB command freqz,phasez.My observation was magnitude response using ibuilt command and own code were same ,whereas the phase response output were different.

If we assume the matlab "Lowpass" compensates for the delay in the filter ,thats the reason for zero shift using code.Then i dont understand how using phasez matlab command produces a linear phase shift when the passing the "lowpass". Its the same lowpass used in both cases .

This lead to understanding lowpass function in matlab ,there msut be more to it & why its producing two different results.??

I have attached the results frequency response of lowpass filter using code in my question above . Frequency response results using Matlab freqz and phasez ,i will attach in here.

fs=100;
f=1:15;
t = 0:1/fs:2-1/fs; 
x= sin(2*pi*t'*f);
[Y,d]= lowpass(x,10,fs);
figure;
[h1,f] = freqz(d,2^16,fs);
subplot(2,1,1);
plot(f,mag2db(abs(h1)));
title('Frequency response of the filter using MATLAB function')
xlabel('Frequency');
ylabel ('Magnitude in db' );
subplot(2,1,2);
[q,f]=phasez(d,2^16,fs);
plot(f,q);
title('Phase response of the filter using MATLAB function')
xlabel('Frequency');
ylabel ('phase in radians' );

Thank you paul,

Hope you can help me in undertanding this !

Answer 1

[Paul]

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))])

m_comp = 1×3

-2.4618 -2.5030 -2.4806

% compare true phase to FFT ratio
p_comp = 180/pi*angle([h Y_fft/U_fft])

p_comp = 1×2

-128.5455 -128.2826

% 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])

[Suvvi Kuppur Narayana Swamy]

I am very sorry for the late reply.But this explanation really has solved my problem and provided a better understanding of frequency response of filters. I am very new to the matlab forum,Could you please guide me on how to accept your anwer here.

Thanks a ton ,

Best Regards,

Suvvi

[Paul]

There should be button to click to accept an answer. Don’t know exactly where, but it should be straightforward.

Answer 2

[Robert U]

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:

https://de.mathworks.com/matlabcentral/answers/508696-make-a-bode-plot-without-bode?s_tid=srchtitle

Kind regards,

Robert

[Suvvi Kuppur Narayana Swamy]

Thank you for your answer. i now understand why the phase is always zero .I read the example attached above. You have used the frequency vector and i tried same in my code ,inorder to find the phase . I get an error like this

"Array indices must be positive integers or logical values".

But could you please tell me how do i now get the phase response without using matlab functions.

Thank you,

Suvvi

[Robert U]

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