How to get from discrete time impulse response vector to plot of frequency spectrum

Hello,

I'm very new to DSP concepts and signal processing in MATLAB so please go easy on me. I have the following code that I've begun working on which is the BZT of an RLC filter. I want to take the DFT of my impulse response y[n] and end up with a magnitude plot Y(f). I know that I can use the fft() and plot() functions for this.

I've seen examples using linspace/logspace vectors which I'll need to define before being able to plot this thing, and some such things as NFFTs (to enable a faster fft?) but I'm not really sure how to put it all together here.

Any help/insight would be really appreciated.

Thank you, Joe.

CODE:

R=0.1; % resistor

L=0.0016; % inductor

C=0.0016; % capacitor

resF=1/sqrt(L*C); % resonant frequency

sampRate=44100; % sampling rate (Hz)

T=1/sampRate; % sample period

prewarpF=(2/T)*tan(resF*T/2); % prewarp factor

a=(2/T)*R*C;

b=(2/(T*prewarpF))^2;

e0=a/(a+b+1);

e1=(2-2*b)/(a+b+1);

e2=(1+b-a)/(a+b+1);

N=2048; % number of taps

x=zeros(1,N); % input vector

x(5) = 1; % impulse at delta(n-5)

y=zeros(1,N); % output vector

for n = 5 : N

    y(n) = e0*x(n) - e0*x(n-2) - e1*y(n-1) - e2*y(n-2);

end

[Joseph]

If anyone is able to respond about how to define my frequency vector according to the parameters already defined in my code, that would be really useful to me. Insight would be as nice as just the code.

Answer 1

[Wayne King]

Since you have a discrete-time sequence with data sampled at 44.1 kHz, your frequency vector will run from -22050 Hz to 22050 Hz for a double-sided spectrum, or just 0 to 22050 for a single-sided.

Do you have the Signal Processing Toolbox? If so, then use freqz() with the coefficients of your difference equation. From this line:

 y(n) = e0*x(n) - e0*x(n-2) - e1*y(n-1) - e2*y(n-2);

You difference equation has the following form as a rational transfer function.

A = [1 e1 e2];
B = [e0 -e0];

So the frequency response is

[H,F] = freqz(B,A,[],44100);
plot(F,20*log10(abs(H)));
xlabel('Hz');

You can see from the above that the response is lowpass.

[Joseph]

Beautiful, thanks.

[Abdullah Ashraf]

@Wayne King I logged in just to thank you! it helped me so much with my assignment, so again, THANK YOU!!

I have another question, can you tell me what does the square bracket notation in this expressions mean please?

[H,F] = freqz(B,A,[],44100);