matlab 作出信号频谱图

如题所述

先对信号等时间采样得到一组时域信号然后做傅里叶变换。特殊情况下可以看出数据点所满足的解析式,使用拟合,然后对拟合得到的函数进行傅里叶变换,用matlab的fourier函数即可。一般情况下得到的离散的数据点没有明显的拟合函数,这时候可以考虑用离散傅里叶变换。matlab中的fft函数可以完成这个功能。


由于一般情况中的fft更具有应用性,下面着重举例说明fft。

引用一段matlab帮助文件提供的代码作说明:(%后面是中文或英文注释)

clc;clear;

Fs = 1000;                    % Sampling frequency,取样频率

T = 1/Fs;                      % Sample time,采样时间间隔

L = 1000;                     % Length of signal,总时间

t = (0:L-1)*T;                % Time vector,时间向量

% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid 信号函数,提供50Hz和120Hz的主频率

x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); 

y = x + 2*randn(size(t));     % Sinusoids plus noise 信号函数加上模拟的噪音

plot(Fs*t(1:100),y(1:100))    % 信号图

title('Signal Corrupted with Zero-Mean Random Noise')

xlabel('time (milliseconds)')

NFFT = 2^nextpow2(L); % Next power of 2 from length of y 扩充采样点,由1000变为1024

Y = fft(y,NFFT)/L;%除一个L,使归一化,可以不除,不影响对主频率的判断

f = Fs/2*linspace(0,1,NFFT/2+1);%这里除以2是因为fft的对称性,因此只画一半

% Plot single-sided amplitude spectrum.

figure

plot(f,2*abs(Y(1:NFFT/2+1))) %乘2是为了归一化,因为右边一半的fft图像没画;不乘,不影响对主频率的判断

title('Single-Sided Amplitude Spectrum of y(t)')

xlabel('Frequency (Hz)')

ylabel('|Y(f)|')


这行代码“f = Fs/2*linspace(0,1,NFFT/2+1);”如果不理解可以写成

“f = Fs*linspace(0,1,NFFT);”然后把后面的plot行的乘2去掉,NFFT/2+1也改成NFFT,这就等于没有折叠的状态。

至于为什么对称、为什么表达式是这样,就需要去做DFT数学推导了,这里不做推导。


得图如下:

折叠了的fft图。在50与120Hz处有明显的主峰。


未折叠的fft图。右边两个峰值并没有实际意义,只是由于对称性而得到的。

温馨提示:答案为网友推荐,仅供参考
第1个回答  推荐于2018-07-24
看看这个————matlab如何做频谱分析% ----------------------------------------------------------------------------------------------------------------
% 方案1:“x = a*cos(2*pi*w*t)”的形式:
% ----------------------------------------------------------------------------------------------------------------
% 注意:1.时域的持续时间范围应较大;
% 2.频率w与序列k的对应关系(N为序列总长度):w = 1/dt * k/N;
% 3.采样频率1/dt应大于 w 的2倍
% 4.结果曲线的峰值的横坐标对应的就是w和-w值
% ----------------------------------------------------------------------------------------------------------------
a = 0.75;
w = pi/3;
dt = 0.2;
t = [-30*pi:dt:30*pi];
N = size(t, 2);
x = a*cos(2*pi*w*t);
y = fft(x);
y = fftshift(y);
figure;
subplot(2,1,1);
plot(t,x);
subplot(2,1,2);
plot(1/dt*(-N/2+1:N/2)/N, abs(y));

% ----------------------------------------------------------------------------------------------------------------
% 方案2:“x = a*cos(w*t)”的形式:
% ----------------------------------------------------------------------------------------------------------------
% 注意:1.时域的持续时间范围应较大;
% 2.频率w与序列k的对应关系(N为序列总长度):w = 1/dt *2*pi* k/N;
% 3.采样频率1/dt应大于 w/(2*pi) 的2倍
% 4.结果曲线的峰值的横坐标对应的就是w和-w值
% ----------------------------------------------------------------------------------------------------------------
a = 0.75;
w = pi/3;
dt = 1;
t = [-20*pi:dt:20*pi];
N = size(t, 2);
x = a*cos(w*t);
y = fft(x);
y = fftshift(y);
figure;
subplot(2,1,1);
plot(t,x);
subplot(2,1,2);
plot(1/dt*2*pi*(-N/2:N/2-1)/N, abs(y));

备注:由于使用了fftshift,所以得到的频谱序列关于原点对称,如果不需要负半轴的话自行修改一下就ok了function f=frequency(x,fs)
dtlen=length(x);
t=(0:dtlen-1)/fs;
subplot(211);
plot(t,x);
axis tight;
y=abs(fft(x))*2/dtlen;
ff=(0:dtlen/2-1)*fs/dtlen;
subplot(212);
f=y(1:floor(dtlen/2));
plot(ff,f);
axis tight;
return;本回答被网友采纳
相似回答