How to plot magnitude spectrum of a signal?

I want to plot magnitude spectrum. Let's say I want to generate two input signals with 100 Hz and 200 Hz.

x1 = cos(2*pi*100*[0:1/fsampling:1.23]);
x2 = cos(2*pi*200*[0:1/fsampling:1.23]);
x = x1 + x2;
x(end) = [];
[b,a] = butter(2,[0.6 0.7],'bandpass');
filtered_noise = filter(b,a,randn(1, length(x)*2));
y = (x + 0.5*filtered_noise(500:500+length(x)-1))/length(x)*2;

%plot first half of DFT (normalised frequency)
Y_mags = abs(fft(y));
num_bins = length(Y_mags);
plot([0:1/(num_bins/2 -1):1], Y_mags(1:num_bins/2)),grid on;
title('Magnitude spectrum of noisy signal');
xlabel('Normalised frequency (\pi rads/sample)');
ylabel('Magnitude');

My problem now is I don't really understand how the y-axis works. Why both 100 Hz and 200 Hz signals have magnitude of 1? Please help me!

Answer 1

[Star Strider]

They have magnitudes of 1 because looking at your code, you defined them that way:

x1 = cos(2*pi*100*[0:1/fsampling:1.23]);
x2 = cos(2*pi*200*[0:1/fsampling:1.23]);

[Nur Fauzira Saidin]

Ahh okay I got it. Thank you much for replying!

[Nur Fauzira Saidin]

But I have another question. Let's say I want to generate another 50Hz input signal. I use the same code.

x1 = cos(2*pi*50*[0:1/fsampling:1.23]);
x2 = cos(2*pi*100*[0:1/fsampling:1.23]);
x3 = cos(2*pi*200*[0:1/fsampling:1.23]);

Why is that the magnitude for 50Hz signal only 0.6?

[Star Strider]

Since we are dealing with discretely-sampled signals with non-linearities that the sampling process introduces, (rather than continuous signals), some of the energy of your signals (that are all harmonically related) ‘leak’ into your other signals. The presence of the noise also affects the total amount of ‘energy’ in the signal, so the amplitudes of the signals is reduced accordingly. (It is likely possible to demonstrate this analytically with the appropriate trigonometric identities and sufficient maths. It may also be in the better signal processing textbooks.) There are probably also computational effects in calculating the fft that would require more analysis to investigate and verify.

If you eliminate the noise (as an experiment), and use signals that not harmonically-related, all the signal amplitudes are equal to 1, as they should be. You will get different results for different frequency signals.

The code that demonstrates this:

fsampling = 1000;
x1 = cos(2*pi*50*[0:1/fsampling:1.23]);
x2 = cos(2*pi*100*[0:1/fsampling:1.23]);
x3 = cos(2*pi*200*[0:1/fsampling:1.23]);
x(1,:) = x1 + x2 + x3;
x4 = cos(2*pi*30*[0:1/fsampling:1.23]);
x5 = cos(2*pi*100*[0:1/fsampling:1.23]);
x6 = cos(2*pi*160*[0:1/fsampling:1.23]);
x(2,:) = x4 + x5 + x6;
[b,a] = butter(2,[0.6 0.7],'bandpass');
filtered_noise = filter(b,a,randn(1, length(x)*2));
% y = (x + 0.5*filtered_noise(500:500+length(x)-1))/length(x)*2;
y = x';

Fn = fsampling/2;
Fy = fft(y)/length(x);
Fv = linspace(0, 1, fix(length(y)/2) + 1)*Fn;
Iv = 1:length(Fv);

figure(1)
subplot(2,1,1)
plot(Fv, abs(Fy(Iv,1))*2)
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Harmonically-Related Signals (No Noise)')
axis([xlim    0  1.5])
grid
subplot(2,1,2)
plot(Fv, abs(Fy(Iv,2))*2)
xlabel('Frequency (Hz)')
ylabel('Magnitude')
title('Harmonically-Unrelated Signals (No Noise)')
axis([xlim    0  1.5])
grid

I used a slightly different fft call and plot than you did, because I am used to it and can interpret it more easily.

Answer 2

[xiaodong lu]

And with zero-padding, one can limit the spectrum leakage effect. For example, the total signal sample points is N, one can do FFT with 2^fftpoints > or ~= 4N. the modified code are:

fsampling = 1000;

x1 = cos(2*pi*50*[0:1/fsampling:1.23]);

x2 = cos(2*pi*100*[0:1/fsampling:1.23]);

x3 = cos(2*pi*150*[0:1/fsampling:1.23]);

x(1,:) = x1 + x2 + x3;

x4 = cos(2*pi*30*[0:1/fsampling:1.23]);

x5 = cos(2*pi*100*[0:1/fsampling:1.23]);

x6 = cos(2*pi*160*[0:1/fsampling:1.23]);

x(2,:) = x4 + x5 + x6;

y = x';

Fn = fsampling/2;

fftpoints=4096;

Fy = fft(y,fftpoints)/length(x);

Fv = linspace(0, 1, fix(fftpoints/2) + 1)*Fn;

Iv = 1:length(Fv);

figure(1)

subplot(2,1,1)

plot(Fv, abs(Fy(Iv,1))*2)

xlabel('Frequency (Hz)')

ylabel('Magnitude')

title('Harmonically-Related Signals (No Noise)')

axis([xlim 0 1.5])

grid

subplot(2,1,2)

plot(Fv, abs(Fy(Iv,2))*2)

xlabel('Frequency (Hz)')

ylabel('Magnitude')

title('Harmonically-Unrelated Signals (No Noise)')

axis([xlim 0 1.5])

grid

The sample number is still 1231, the FFT length is 4096 (one can use 8192 either). The output spectrum is much better.