How to make the frequency vector to analyze a complex signal?

I have this example and I'm puzzled with the frequency vector (hz) if I change any parameter in it, the result won't be the frequency components. How to calculate the frequency vector, generally?

srate = 1000;

time = -2: 1/srate : 2;

pnts = length(time);

hz = linspace(0,srate/2,(pnts-1)/2);

freq = [2 10 25 40];

amp = [3 5 7 9];

phas = [-pi/4 pi/2 pi 2*pi];

signal = zeros(size(time));

for i = 1:length(freq)
    
    signal = signal + amp(i)*sin( 2*pi*freq(i)*time + phas(i));
    
end

subplot(211)
plot(time,signal,'k-','linew',0.5)
title("Time Domain")

y = fft(signal)/length(signal);

pwr = abs( y ).^2;
pwr(1) = [];

subplot(212)
plot(hz,pwr(1:length(hz)),'r-')
title("Frequency Domain")

[Michael Soskind]

I'm a bit puzzled by "if I change any parameter in it, the result won't be the frequency components". If you mean that changing the sample rate gives an incorrect frequency axis, or that the points do not always work. What you do is what should be done, but the only change I advise is replacing "(pnts-1)/2" with "floor(points/2)", in order to account for even numbers of points in the time/signal vectors. Other than that, the single sided FFT that you show in this example is reasonable, although the powers that you are getting should be pwr = 2*abs(y), rather than pwr = abs(y).^2, based on the MATLAB description of fft here.

[Muhamed Sewidan]

Thank you for your response.

I mean by changing srate/2 or (pnts-1)/2, I don't get in the end the values of the frequencies.

Pardon me, but my second question why should I use pwr = 2*abs(y) and not pwr = abs(y).^2 ?

Thank you

[Michael Soskind]

I see what you mean, you are referring to the small errors in frequency value compared to the ideal frequency values. If you increase the sample rate, you will find that the approximaton of the DFT should improve. If I increase srate to srate = 10000, I get the pairs:

freq: 2, amp: 3

freq: 10, amp: 5

freq: 25, amp: 6.99

freq: 40, amp: 9

With respect to the second question, pwr = 2*abs(y) should be used since this is really the folding of the FFT on itself in a single sided Fourier Transform. In order to get the correct signal amplitude as you had built in your composite signal, you should just double the magnitude, which will get you the correct signal amplitudes.

[Muhamed Sewidan]

Thank you you've been a great help to me.

Answer 1

[Cris LaPierre]

I think I would explain it this way. Due to Nyquist theory, the highest frequency you can detect is half your sample rate (must have at least 2 points in one period). The lowest is obviously 0. The assumption here is that the detectable frequencies increase linearly from 0 to srate/2.

If you inspect the peaks in your plot, you'll notice this approach is an approximation, as the peaks are about 1 point off.

[Muhamed Sewidan]

Thank you for your help.

This makes sense about the second argument in linspace (srate/2) but must the number of points be this value (pnts-1)/2 ?

[Cris LaPierre]

I believe that was because there is an odd number of points in time (4001). So (pnts-1)/2 is a nice round number of 2000.

[Muhamed Sewidan]

Yes I know but rather than using (pnts-1)/2 which equals 2000. Can I use 1000 for example rather than 2000?

[Cris LaPierre]

Not without having to change some other things. The number of points is set by your sample rate. Changing it directly impacts the results, as you've observed. To get the frequencies to line up, you'd have to half your sample rate, which halves the frequencies you can resolve in the signal.

Better to zoom in on your axis if you want to focus on the peaks. Use xlim.

srate = 1000;
time = -2: 1/srate : 2;
pnts = length(time);

hz = linspace(0,srate/2,(pnts-1)/2);

freq = [2 10 25 40];
amp = [3 5 7 9];
phas = [-pi/4 pi/2 pi 2*pi];

signal = zeros(size(time));
for i = 1:length(freq)    
    signal = signal + amp(i)*sin( 2*pi*freq(i)*time + phas(i));    
end

y = fft(signal)/length(signal);
pwr = abs( y ).^2;
pwr(1) = [];

plot(hz,pwr(1:length(hz)),'r-')
title("Frequency Domain")
xlim([0 250])

[Muhamed Sewidan]

Thank you. Now I get it.