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

Started by September 6, 2004
```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

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,

[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