matlab 作出信号频谱图

Posted 2023-03-17

先对信号等时间采样得到一组时域信号然后做傅里叶变换。特殊情况下可以看出数据点所满足的解析式，使用拟合，然后对拟合得到的函数进行傅里叶变换，用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图。右边两个峰值并没有实际意义，只是由于对称性而得到的。

参考技术A 看看这个————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;