How to plot magnitude spectrum of a signal?

213 views (last 30 days)
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!

Accepted Answer

Star Strider
Star Strider on 30 Nov 2015
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]);
  3 Comments
Nur Fauzira Saidin
Nur Fauzira Saidin on 30 Nov 2015
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
Star Strider on 30 Nov 2015
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.

Sign in to comment.

More Answers (1)

xiaodong lu
xiaodong lu on 8 Aug 2017
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.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!