How to plot magnitude spectrum of a signal?

165 次查看（过去 30 天）
Nur Fauzira Saidin
Nur Fauzira Saidin 2015-11-30
回答： xiaodong lu 2017-8-8
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!

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2015-11-30
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 个评论
显示 1更早的评论
Nur Fauzira Saidin
Nur Fauzira Saidin 2015-11-30
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 2015-11-30
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.

请先登录，再进行评论。

更多回答（1 个）

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

类别

在 Help Center 和 File Exchange 中查找有关 Digital Filter Analysis 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by