Hi! I want to get FRF of IIR Filter from Impulse Response. So, I got FRF using two method. One is FREQZ, another is FFT of impulse response. But two results are different. Why? My source code are below : clear all close all clc fs=2048; % Sampling Frequency N=65536; % Buffer Size t=[0:N-1]'./fs; % Time f = [0:N/2]'*fs/N; % Frequency x = zeros(N,1); % Impulse Input Signal x(1)=1; % 1 pole butterworth highpass filter with 2.7 Hz Cutoff [B, A] = butter(1, 2*2.7/fs, 'high'); [h, f] = freqz(B, A, N, fs); h = abs(h(1:N/2+1)); h_mag = 20*log10(h); f = f(1:N/2+1); % Impulse Response y = filter(B, A, x); % Calculate FFT with impulse Response h_fft = fft(y, N); h_fft = abs(h_fft); h_fft = h_fft(1:N/2+1); figure, semilogx(f, h, 'b', 'LineWidth', 2), hold on, semilogx(f, h_fft, 'r-.', 'LineWidth', 2) xlabel('Frequency(Hz)'), ylabel('Magnitude') grid, legend('From FREQZ Function', 'From Impulse Function') h_fft = 20*log10(h_fft); figure, semilogx(f, h_mag, 'b', 'LineWidth', 2), hold on, semilogx(f, h_fft, 'r-.', 'LineWidth', 2), grid, xlabel('Frequency(Hz)'), ylabel('dB'), title('Compare h1'), legend('FREQZ Function', 'Impulse Function') figure, freqz(B, A, N, fs); In my opinion, the results of FREQZ are correct. And the results of FFT are not correct. Could you say what I wrong in the cource code?
How can I get Frequency Response from Impulse Freponse In Matlab?
Started by ●September 6, 2004
Reply by ●September 7, 20042004-09-07
kimjinki@shinbiro.com (KIM) wrote in message news:<29a47753.0409061839.6388eb2a@posting.google.com>...> Hi!Hi.> I want to get FRF of IIR Filter from Impulse Response. So, I got FRF > using two method. One is FREQZ, another is FFT of impulse response. > But two results are different. Why? > > My source code are below :Let's just see that I understand what you do. 1) Clear the workspace and initialize variables.> clear all > close all > clc > > fs=2048; % Sampling Frequency > N=65536; % Buffer Size > t=[0:N-1]'./fs; % Time > f = [0:N/2]'*fs/N; % Frequency > > x = zeros(N,1); % Impulse Input Signal > x(1)=1;2) Find the parameters of a high-pass IIR filter. The filter is supposed to be a high-pass fiilter of order 2, but what is the cut-off frequency supposed to be? I do not understand the way you specify the Wn parameter to BUTTER. Would this be the cause of the problems?> % 1 pole butterworth highpass filter with 2.7 Hz Cutoff > [B, A] = butter(1, 2*2.7/fs, 'high');3) Use FREQZ toplot its frequency response...> [h, f] = freqz(B, A, N, fs);... and prepare to plot it in logarithmic scale.> h = abs(h(1:N/2+1)); > h_mag = 20*log10(h); > f = f(1:N/2+1);4) Now for the impulse response version. Find the impulse response by filtering the impulse sequence.> % Impulse Response > y = filter(B, A, x); > % Calculate FFT with impulse Response > h_fft = fft(y, N); > h_fft = abs(h_fft); > h_fft = h_fft(1:N/2+1);The rest is plotting.> figure, > semilogx(f, h, 'b', 'LineWidth', 2), hold on, > semilogx(f, h_fft, 'r-.', 'LineWidth', 2) > xlabel('Frequency(Hz)'), ylabel('Magnitude') > grid, legend('From FREQZ Function', 'From Impulse Function') > > h_fft = 20*log10(h_fft);The reference level for the dB scale is left unspecified here. This could account for some of the difference, but it should be a constant offset over the whole frequency band. Your code produces a 6 dB offset in the transition band and no offset in the pass-band, indicating that the scaling is not part of the problem.> figure, > semilogx(f, h_mag, 'b', 'LineWidth', 2), hold on, > semilogx(f, h_fft, 'r-.', 'LineWidth', 2), > grid, xlabel('Frequency(Hz)'), ylabel('dB'), > title('Compare h1'), legend('FREQZ Function', 'Impulse Function') > > figure, > freqz(B, A, N, fs); > > In my opinion, the results of FREQZ are correct. And the results of > FFT are not correct. > > Could you say what I wrong in the cource code?I am a bit confused by the way you specify the parameters to the BUTTER function. This is probably due to the fact that I have never used the BUTTER function myself, and do not know how it works. So I would have checked out that function call first, to see that is is correct. First, you specify the order to be 1. Is this the order of the analog prototype or the digital equivalent? If the order applies to the analog prototype, everything is fine. If it applies to the discrete equivalent, you may have a problem, since the bilinear transform of a 1st order analog protoype produces a 2nd order discrete equaivalent. Second, I am a bit confused about the way you specify the Wn parameter to the BUTTER function. I did run your script to design a low-pass filter, and saw the same type of effect, so whatever the case is, it is not due to aliasing. I did look into FREQZ.m, and on line 82 the function POLYVAL to compute the frequency response. A few lines further down (line 94) the frequency response is computed by spectral division of the DFTs of the numerator and denumerator in the transfer function. I don't understand how POLYVAL works, so I'll refrain from comment on that. The spectral division, however, seems a bit awkward and I wouldn't use that technique myself. All in all, check that you use the BUTTER filter design function correctly. If you find no mistakes there, I'd trust the fft(impulse) answer and discard the FREQZ answer. Rune
Reply by ●September 7, 20042004-09-07
kimjinki@shinbiro.com (KIM) wrote in message news:<29a47753.0409061839.6388eb2a@posting.google.com>...> f = [0:N/2]'*fs/N; % Frequency > f = f(1:N/2+1);> Could you say what I wrong in the cource code?On second thought, I thing you should check out the frequency axis. It appears you half the frequency range twice, first from fs to fs/2, and then a second time. Once is enough. Rune
Reply by ●September 8, 20042004-09-08
On 6 Sep 2004 19:39:51 -0700, kimjinki@shinbiro.com (KIM) wrote:>Hi! > >I want to get FRF of IIR Filter from Impulse Response. So, I got FRF >using two method. One is FREQZ, another is FFT of impulse response. >But two results are different. Why? > >My source code are below : >clear all >close all >clc > >fs=2048; % Sampling Frequency >N=65536; % Buffer Size >t=[0:N-1]'./fs; % Time >f = [0:N/2]'*fs/N; % Frequency > >x = zeros(N,1); % Impulse Input Signal >x(1)=1; > >% 1 pole butterworth highpass filter with 2.7 Hz Cutoff >[B, A] = butter(1, 2*2.7/fs, 'high'); > >[h, f] = freqz(B, A, N, fs); > >h = abs(h(1:N/2+1)); >h_mag = 20*log10(h); >f = f(1:N/2+1); > >% Impulse Response >y = filter(B, A, x); >% Calculate FFT with impulse Response >h_fft = fft(y, N); >h_fft = abs(h_fft); >h_fft = h_fft(1:N/2+1); > >figure, >semilogx(f, h, 'b', 'LineWidth', 2), hold on, >semilogx(f, h_fft, 'r-.', 'LineWidth', 2) >xlabel('Frequency(Hz)'), ylabel('Magnitude') >grid, legend('From FREQZ Function', 'From Impulse Function') > >h_fft = 20*log10(h_fft); > >figure, >semilogx(f, h_mag, 'b', 'LineWidth', 2), hold on, >semilogx(f, h_fft, 'r-.', 'LineWidth', 2), >grid, xlabel('Frequency(Hz)'), ylabel('dB'), >title('Compare h1'), legend('FREQZ Function', 'Impulse Function') > >figure, >freqz(B, A, N, fs); > >In my opinion, the results of FREQZ are correct. And the results of >FFT are not correct. > >Could you say what I wrong in the cource code?Hello, replace your first [h, f] = freqz(B, A, N, fs); command with [h, f] = freqz(B, A, N, fs,'whole'); and your code will then work. (Ah yes, please don't forget the thank Rune for tryin' to help. Anytime you have a guy with his skills help you with a problem, you're lucky.) [-Rick-]