当前位置：首页 > article >正文

MATLAB与Python中的快速傅里叶变换频谱分析

article 2025/1/14 9:28:17

MATLAB与Python中的快速傅里叶变换频谱分析

Python和MATLAB都可以方便地实现FFT频谱分析

Python Code

# -*- coding: utf-8 -*- 
"""
@File    :   python_fft.py
@Time    :   2025/01/12 19:46:39
@Version :   
@Desc    :   频谱分析——FFT
"""



import numpy as np                              # 使用numpy库中的fft
import matplotlib.pyplot as plt                 # 绘制图像支持

from scipy.signal import butter, filtfilt       # 信号处理: 滤波


# 生成信号
fs = 1000  # 采样频率 (Hz) 大于奈奎斯特采样频率

# 生成等间隔时间向量, 包含fs个等间隔时间点, 并且不包含最后一个点
# 设置endpoint=False可以确保时间向量的长度与采样点数fs一致
t = np.linspace(0, 1, fs, endpoint=False)       # 时间向量 (1秒)
f1 = 50  # 信号频率 (Hz)
f2 = 120  # 信号频率 (Hz)
signal = 0.7 * np.sin(2 * np.pi * f1 * t) + 1.0 * np.sin(2 * np.pi * f2 * t)  # 合成信号

# 添加噪声
noise = 0.5 * np.random.normal(size=len(t))     # 标准正态分布的噪声
signal += noise

# 计算 FFT
n = len(signal)  # 信号长度
fft_result = np.fft.fft(signal)  # FFT 计算
frequencies = np.fft.fftfreq(n, d=1/fs)  # 频率向量

# 取正频率部分 截取
positive_freq = frequencies[:n//2]      
magnitude = np.abs(fft_result[:n//2])  # 取幅值

# 绘制原始信号图像
plt.figure(figsize=(10, 6))
plt.plot(t, signal)
plt.title('Original Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid()
plt.show()



# 绘制频谱图
plt.figure(figsize=(10, 6))
plt.plot(positive_freq, magnitude)
plt.title('Spectrum Analysis (FFT)')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.grid()
plt.show()


# 滤波处理

# 设计低通滤波器
def butter_lowpass(cutoff, fs, order=5):
    nyquist = 0.5 * fs
    normal_cutoff = cutoff / nyquist
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = filtfilt(b, a, data)  # 零相位滤波
    return y

# 应用低通滤波器
cutoff = 80  # 截止频率 (Hz)
filtered_signal = lowpass_filter(signal, cutoff, fs)

# 绘制滤波后波形图
plt.figure(figsize=(10, 6))
plt.plot(t, filtered_signal, label='Filtered Signal')
plt.title('Filtered Signal (Lowpass)')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.grid()
plt.legend()
plt.show()

MATLAB Code


% 编码格式 : ---> -*- utf-8 -*-
% @Time       : 2025-01-12 19:46:44
% @Author     : 不说(bushuo)
% @Contact    : 3513566868@qq.com
% Description : 


% 生成信号
fs = 1000;  % 采样频率 (Hz)
t = 0:1/fs:1-1/fs;  % 时间向量 (1秒)
f1 = 50;  % 信号频率 (Hz)
f2 = 120;  % 信号频率 (Hz)
signal = 0.7 * sin(2 * pi * f1 * t) + 1.0 * sin(2 * pi * f2 * t);  % 合成信号

% 添加噪声
noise = 0.5 * randn(size(t));
signal = signal + noise;                         % matlab 不支持 += 这种语法

% 计算 FFT
n = length(signal);  % 信号长度
fft_result = fft(signal);  % FFT 计算
frequencies = (0:n-1) * (fs/n);  % 频率向量

% 取正频率部分
positive_freq = frequencies(1:n/2);
magnitude = abs(fft_result(1:n/2));  % 取幅值




% 绘制原始信号图像
figure;
plot(t, signal);
title('Original Signal');
xlabel('Time (s)');
ylabel('Amplitude');
grid on;


% 绘制频谱图
figure;
plot(positive_freq, magnitude);
title('Spectrum Analysis (FFT)');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
grid on;


% 滤波处理

% 设计低通滤波器
cutoff = 80;  % 截止(截断)频率 (Hz)
order = 5;  % 滤波器阶数
[b, a] = butter(order, cutoff/(fs/2));  % 巴特沃斯低通滤波器

% 应用低通滤波器
filtered_signal = filtfilt(b, a, signal);  % 零相位滤波


% 绘制滤波后信号
figure();
subplot(3, 1, 3);
plot(t, filtered_signal);
title('Filtered Signal (Lowpass)');
xlabel('Time [s]');
ylabel('Amplitude');
grid on;