Matlab和Python中使用fft函数绘制频谱(个人学习记录)

夜曲2021_

已于 2024-04-22 15:59:19 修改

阅读量1.2w

点赞数 27

文章标签： matlab 开发语言

于 2022-08-23 10:31:43 首次发布

本文链接：https://blog.csdn.net/fightclub_/article/details/126477236

版权

省流：最简洁的版本：

Matlab代码：

%输入：需要求单边频谱的序列y及其采样频率fs
%输出：单边频谱图（幅值已修正）
y_fft=fft(y);
P2 = abs(y_fft/n); %n为序列的长度
P1 = P2(1:n/2+1);  %取单边谱
P1(2:end-1) = 2*P1(2:end-1); %幅值修正，第一个0频率点和最后一个点n/2+1均不需要乘2
f = fs*(0:(n/2))/n; %横轴频率向量
plot(f,P1);

Python代码：

# 输入：信号data，采样频率fs
# 输出：单边频谱（幅值已修正）
n = len(data)
y_fft = np.fft.fft(data)
P2 = np.abs(y_fft / n)
P1 = P2[0:int(n / 2) + 1]  # 取单边谱
P1[1:-1] = 2 * P1[1:-1]  # 幅值修正，第一个0频率点和最后一个点n/2+1均不需要乘2
f = fs * np.arange(0, (n / 2) + 1) / n  # 横轴频率向量
plt.plot(f, P1)
plt.show()

详细推导版本(以matlab代码为例)：

%%
%生成序列信号
Fs = 1000;            % 信号的采样频率(1s采样的点数)，根据采样定理应大于信号中最高频率的2倍                    
T = 1/Fs;             % 相邻采样点之间的采样间隔       
L = 1500;             % 序列长度，也即采样点数
t = (0:L-1)*T;        % 时间向量。此序列对应的采样时长=采样间隔*采样点数
t1=0:1/Fs:1-1/Fs;     % 这种常见的时间向量表示指的是单位时间(1s)对应的时间向量。等价于T*(0:Fs-1) 
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); %序列信号，频率成分分别为50hz和120hz
S1=0.7*sin(2*pi*50*t1) + sin(2*pi*120*t1); 
figure
plot(t,S)
figure
plot(t1,S1)

%求真实单边频谱
L=length(S);
S_fft=fft(S);   %默认做S的序列长度1500点的FFT，结果为1*1500点的复数，a+bi
                %该复数的模为幅值(但不是真实幅值)，相位为arctan(b/a)
S_shuangbianpu=abs(S_fft/L);%abs求复数的摸，得到双边幅度谱。且进行了幅值修正1(从公式角度)：由IDFT的公式可知，
                %公式前面有个系数1/N，N为DFT的点数，而fft直接得到的结果没有乘这个系数，因此需要乘上
S_danbianpu=S_shuangbianpu(1:L/2+1);%得到单边谱，多取一个点是因为第一个点为0hz频率点，即直流分量，真实频率点后延一个

%幅值修正2(单边谱和双边谱幅值的关系角度)：已知双边谱的幅度等于单边谱幅度的一半，故单边谱的真实幅度还需要乘2，
%但是0频率点处的双边谱幅值=单边谱幅值，不需要乘。且第L/2+1个点也不需要乘：此点被正负频率点共用，幅值已经变成了2倍，正好等于单边谱的幅值
S_danbianpu(2:end-1)=2*S_danbianpu(2:end-1);
S_danbianpu(1)=0; %去直流分量，可以直接置零也可以预处理时令序列减去其均值
f=Fs*(0:L/2)/L; %横轴频率向量，单位hz，Fs/L为频率分辨率(每两个频点间的频率间隔)，此例为1000/1500=2/3 hz，(0:750)*2/3 = (0:500)hz
figure
plot(f,S_danbianpu) 
title('单边谱')
xlabel('频率/f (Hz)')
ylabel('幅值')

序列长度为1.5s：