How can I get Frequency Response from Impulse Freponse In Matlab?

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?

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

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

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