# DTFT of sine wave using freqz

As mentioned in the title, is it possible to use freqz to find the DTFT of a sine wave? I am confused about what the 'a' and 'b' vectors would look like, since there are only impulses in the numerator.

• hm, why do you think freqz is the tool of choice here? – Marcus Müller Feb 17 at 14:27
• Yes, it's possible. No, it's not a good idea. Needlessly complicated and not intended for this application – Hilmar Feb 17 at 18:00

## 2 Answers

The MATLAB / OCTAVE function freqz(b,a) is typically used for displaying the magnitude and the phase of a Frequency Response function $$H(e^{j\omega})$$ associated with a discrete-time, linear time-invariant (LTI), system with rational transfer function $$H(z)$$ (having a corresponding LCCDE representation) given by

$$H(z) = \frac{ \sum\limits_{k=0}^{M} b_k z^{-k} }{\sum\limits_{k=0}^{N} a_k z^{-k} } \tag{1}$$

where the numerator and denominator coefficients $$b_k$$ and $$a_k$$ represent such systems.

It's not meant to display spectrum of signals. Of course, an impulse response $$h[n]$$ associated with an LTI system can also be considered as a signal, and therefore its Frequency Response function $$H(e^{j\omega})$$ can also be considered as a Fourier transform of a signal ($$h[n]$$ in this case) too.

However, for such an interpretation to make sense, the signal $$h[n]$$ should be the impulse response of an LTI system. Stated in other words, to display the spectrum of $$h[n] = \cos(\omega_0 n )$$, you should find the coefficients $$a_k$$ and $$b_k$$ that constructs a LCCDE whose impulse response solution is $$h[n]$$.

Consider a first order LCCDE (with initial rest):

$$y[n] - A y[n-1] = x[n] \tag{2}$$

taking its Z-transform

$$Y(z)(1 - A z^{-1}) = X(z) \tag{3}$$

and the System Function is

$$H(z) = \frac{Y(z)}{X(z)} = \frac{1 }{1 - A z^{-1}} \tag{4}$$

and the impulse response is (assuming causal)

$$h[n] = A^n u[n] \tag{5}$$

If we let $$A = e^{j\omega_0}$$ , then the impulse response becomes

$$h[n] = e^{j \omega_0 n} u[n] \tag{6}$$

From Eq-6, we have a single complex-exponential, by the choice of the coefficient $$A =e^{j\omega_0}$$. When combined with Euler identity of a causal cosine wave, we can see that

$$h[n] = \cos(\omega_0 n) u[n]= 0.5( e^{j \omega_0 n} + e^{-j \omega_0 n} ) u[n] \tag{7}$$

Which can also be written with $$A$$ as

$$h[n] = \cos(\omega_0 n) = 0.5( A^n + (A^*)^n ) u[n] \tag{8}$$

where $$A^*$$ is the complex-conjugate of $$A$$. Using the relation between Eqs.4&5, we can find the Z-transform of $$h[n]$$ as

$$H(z) = \frac{0.5}{1- A z^{-1}} + \frac{0.5}{1 - A^* z^{-1} } \\ \tag{9}$$

equivalently

$$H(z) = \frac{1 - 0.5(A+A^*) z^{-1}}{ 1- (A + A^*) z^{-1} + |A|^2 z^{-2} } \\ \tag{10}$$

we can construct back the differential equation that corresponds to this $$H(z)$$ as

$$y[n] - 2\mathcal{Re}\{A\} y[n-1] + |A|^2 y[n-2] = x[n] -\mathcal{Re}\{A\} x[n-1] \\ \tag{11}$$

Plugging $$A = e^{j\omega_0 n}$$ into Eq-11 yields

$$y[n] - 2 \cos(\omega_0) y[n-1] + y[n-2] = x[n]-\cos(\omega_0) x[n-1] \tag{12}$$

from Eq-12, you find the corresponding coefficients $$a_k$$ and $$b_k$$ as:

$$b_k = [1, -\cos(\omega_0)] \tag{13}$$

and

$$a_k = [1, -2 \cos(\omega_0), 1] \tag{14}$$

Now, you can display the Fourier transform of $$h[n] = \cos(\omega_0 n)$$ by using the freqz(b,a) function, with the coefficeints in Eqs-13,14.

Note that, $$h[n]$$ is only a causal sine wave valid for $$n>0$$. It's not the ideal two sided sine wave. Furthermore, I hope you are aware that the ideal sine wave's Fourier transform is a dirac impulse function. And what you can view freqz(b,a) or any other function is actually the spectrum of teh windowed cosine wave of a finite length.

The following MATLAB/ OCTAVE code shows you the implementation of above:

% spectrum of h[n] = cos(w0*n) to be displayed by freqz(b,a) function

N = 128;               % number of samples to evaluate the impulse response

w0 = 3*(2*pi/N);       % frequency of the cosine h[n] = cos(w0*n)
b = [1,-cos(w0)];      % coefficients b[k] of the LCCDE
a = [1,-2*cos(w0),1];  % coefficients a[k] of the LCCDE

h = impz(b,a,N);       % N samples of h[n] in n=[0:N-1] evaluated

figure,stem([0:N-1],h);
title('imp-resp of y[n] - 2 cos(\omega_0) y[n-1] + y[n-2] = (1-cos(\omega_0)) x[n]');

K = 1024;              % number of frequency points to be evaluated.
figure,freqz(b,a,K);
title('Freq-Response of the system: y[n] - 2 cos(\omega_0) y[n-1] + y[n-2] = (1-cos(\omega_0)) x[n]');


The outputs look like :

and

Note that I have chosen a perfect frequency for the cosine wave, so that it displayed almost like an impulse. It will not be the case in general.

In short, and importantly, YES you can absolutely and simply use freqz if desired to compute the samples on the DTFT of a signal. The following produce identical results:

DTFT_1 = freqz(signal, 1, N, 'whole')
DTFT_2 = fft(signal, N)


Note the second form is simpler and is the zero padding of the waveform signal out to N samples, which produces N samples on the DTFT. freqz as shown produces the identical result.

But the form

DTFT_3 = freqz(signal)


just returns the first half of DTFT_2 for the default case of 512 samples so can be just as simple when that is sufficient.

As for answering does it make sense to use freqz for that purpose, in some simple profiling, the freqz() approach takes approximately 3 times longer than fft(), so I would say no if computation time is of any concern. And similarly for all the same reasons described, if you need to compute the frequency response for an FIR system with coefficients a, and don't need the convenient formatted plots freqz can provide, then the fir(a, 512) approach would be the preferred method over using freqz(a) for that purpose as well (Or whatever number of samples desired in the more general freqz(a,1,N) form.)

Bottom line, they are indeed equivalent when used as shown in case that point was missed, and given the convenience of the formatted plots returned by the freqz command, this may be a good enough and practical reason to use freqz for a quick inspection of the DTFT.