0

Suppose I have successfully generated a Single-sided Power spectrum as follows:

X_mags = abs(fft(signal));
bin_vals = [0 : N-1];
fax_Hz = bin_vals*fs/N;
N_2 = ceil(N/2);
plot(fax_Hz(1:N_2), 20*log10(X_mags(1:N_2)));`

now I want to plot a second graph on top of the first one:

hold on;

Lastly I perform an LPC analysis on the signal and calculate the frequency response. The frequency response should be plotted on top of the Power spectrum, so:

[a, g] = lpc(signal,N);
[h,f] = freqz(b,a,N,fs);
plot(?);

For the sake of simplicity, let's assume the parameters are all correctly given, how should I write the plot function for having a correct display of the frequency response? A simple plot(f) doesn't work.

Can someone explain why? Thanks

user3488736
  • 109
  • 1
  • 3
  • 13

2 Answers2

0

The response h is complex, so you'll want to get the magnitude of the response via multiplication by its complex conjugate.:

plot(f, 10*log10(h.*conj(h)));

Note the use of 10*log10 because the operation above squares the amplitude response held in h.

Or, if you're less interested in being pedantic about working with complex numbers, you could just take the absolute value of the complex values, being sure to 20*log10 because the abs does not square the values

plot(f, 20*log10(abs(h)));
chipaudette
  • 1,655
  • 10
  • 13
0

A simple plot(f) tries to plot frequency vector, isn' t it?

Check the code below:

X_mags   = abs(fft(signal));
bin_vals = [0 : N-1];
fax_Hz   = bin_vals*fs/N;
N_2      = ceil(N/2);

[a, g]   = lpc(signal,N);
[h, f]   = freqz(b, a, N, fs);

figure,

hold on,
plot(f, 20*log10(abs(h)), 'r');
plot(fax_Hz(1:N_2), 20*log10(X_mags(1:N_2)));
title('Frequency Spectrum');
xlabel('Frequency (Hz)'); 
ylabel('Amplitude (dB)');
legend('Frequency Response', 'Single-sided Power spectrum')

Btw, there is a MATLAB function db() for decibel calculation. That might be useful.

mehmet
  • 1,631
  • 16
  • 21
  • I tried it but the plot of the lpc is untill fs/2? What is the reason of it? – user3488736 Dec 13 '14 at 15:29
  • And There is a command like that: `plot(fax_Hz(1:N_2), 20*log10(X_mags(1:N_2)));` – mehmet Dec 13 '14 at 15:37
  • I know, so that I have a single side of the power spectrum. But I do not use N_2 for the lpc, so why is it affected by it? EDIT: nevermind, because of the to low order of the FIR of the LPC. By increasing the prediction order, I get better results (which is logical). Though is it pure coincidence that if N=2 the frequency response is untill fs/2? :) – user3488736 Dec 13 '14 at 15:38
  • 1
    Oh sorry. That is caused by `freqz()` function. It says **In all cases, the frequency response is evaluated at N points equally spaced around the upper half of the unit circle. If N isn't specified, it defaults to 512.** You can add `'whole'` to `freqz()` function as an argument.. – mehmet Dec 13 '14 at 15:49
  • I just figured the same indeed. Thanks Mehmet :) – user3488736 Dec 13 '14 at 15:56