matlab将示波器得到的信号进行FFT处理

jayson.h

已于 2024-12-16 13:26:44 修改

阅读量1.5k

点赞数 19

文章标签： matlab 数据库开发语言

于 2024-12-13 18:05:04 首次发布

本文链接：https://blog.youkuaiyun.com/fenxiang_Hao/article/details/144457027

版权

1FFT

1.1什么是fft

FFT（Fast Fourier Transformation）就是快速傅里叶变换的意思。输入的是离散数据，输出的也是离散频率。

在matlab中具体常用的使用方法为X=fft(x)或X=fft(x,Ns)。

1.2fft信号长度Ns

计算信号的fft，长度取信号长度。

Ns=length(t)

信号分析长度：一般取信号长度，或者2的整数次方长度。如果大于信号长度，

matlab默认在信号后面补0。

1.3计算频率索引号

fft得到的数据具有左右对称的特点，通常只采用左半边的数据作为信号数据。

同时提一句，fft得到的数据长度同Ns。

1.4计算真实幅值

abs(X)*2/Ns这个是固定的，记住就行了，如果想知道为什么可以了解fft的具体算法。这里Ns取得是信号的有效长度，如果在信号末尾补0的话，不算做有效长度。

1.5设置频域刻度

freq = (0:Ns/2)*fs/Ns;这个也是固定的用法。唯一注意的就是要和频率索引对应上。

通过计算可以看到，频率的最大值为fs/2，频率的分辨率为fs/Ns=1/T。

1.6信号采样频率，计算前两个时间点的差的倒数

fs=1/(t(2)-t(1));

信号采样频率：就是每秒钟信号的采样点，与信号采样间隔成倒数关系。

所以为了提高频率的最大值，需要提高信号的采样频率，这个主要是针对信号采样仪器设备的要求。

所以为了提高频率的分辨率，需要增加测量时间T，把信号尽可能测量的更长

7 信号加窗避免信号频谱泄露

如果信号在没有能在时间窗口内出现整数个周期，就会出现信号的泄露。
下图是一个1hz的信号，没有在末尾处没有出现完整，得到的频谱图。

信号泄露会导致该频率峰值的减小，同时也会导致频率峰宽度的增加。如上图，频率峰值不是1，而是比1小，而且1hz在两侧很长一段距离衰减的都很慢。

由于实际测量中，很难测到正好整数个波长，所以我们可以考虑给信号做一定的处理来减小这个泄漏问题。

代码

写代码的逻辑

 假设 t 和 x 已经定义，其中 t 是时间向量，x 是对应的信号值

将示波器的数据导出，只保存成t和x两列数据

matalb里包含了读取该.m文件下的表格的能力

接下来就是：

1.计算数据长度

2 对纵坐标加汉宁窗

3对横坐标进行采样频率的计算，（横坐标相邻两个点时间差的倒数）

4 生成频率轴

5 对横坐标进行FFT转换

计算幅度值，将幅度转换为功率谱密度，

然后再将功率谱密度转换为dBm

6 只取单边带频谱

7作图，sbu(2,1,1)和sbu(2,1,2)，这样就可以将原始信号和FFT处理后的图放到一个表格里

代码1：

% 数据选择界面
[specname,pathname,filterindex] = uigetfile('*.*','选择参数文件','MultiSelect', 'on');
% vm=[];rsd=[];len =[];specnum=[];num=0;numt=0;yy=[];cc=[];list=[];FWHMh=[];vvm=[];
 
specname = fullfile(pathname,specname);
data = readtable(specname);
data = table2array(data);
database = data;
t=data(:,1); %示波器横坐标
x=data(:,2); %示波器纵坐标
 
%先把示波器的数据分两列导入工作区，横坐标命名为t，纵坐标命名为x
subplot(2,1,1);
plot(t,x); %画出原始信号
 
%求FFT
N=length(t); %计算信号长度
Fs=1/(t(2)-t(1));%信号采样频率，计算前两个时间点的差的倒数
 
% 生成频率轴（以Hz为单位），用于后续绘图展示
f = (0:N-1)*(Fs/N);
 
% 应用汉明窗
win = hamming(N);
x_win = x.*win';
 
% 进行FFT
X = fft(x_win);
 
% 计算幅度谱（取绝对值，并归一化，这里的归一化是常规的一种处理方式）
X_mag = abs(X)/N;
 
% 将幅度转换为功率谱密度（根据 Parseval 定理等相关原理）
Pxx = (X_mag.^2);
 
% 将功率谱密度转换为dBm形式（参考功率为1mW，这里假设信号是在50欧姆负载上，可根据实际调整）
R = 50; % 负载电阻（欧姆）
Pref = 1e-3; % 参考功率（W），1mW转换为W
dBm = 10*log10(Pxx*R/Pref);
 
% 只取单边频谱（因为实际信号的频谱是对称的，通常我们只看单边就够了）
f_single_side = f(1:N/2+1);
dBm_single_side = dBm(1:N/2+1);
 
% 绘制频谱图
subplot(2,1,2);
plot(f_single_side, dBm_single_side);
xlabel('Frequency (Hz)');
ylabel('Power (dBm)');
title('Frequency Spectrum');
grid on;

代码2：

% 数据选择界面
[specname,pathname,filterindex] = uigetfile('*.*','选择参数文件','MultiSelect', 'on');
% vm=[];rsd=[];len =[];specnum=[];num=0;numt=0;yy=[];cc=[];list=[];FWHMh=[];vvm=[];
 
specname = fullfile(pathname,specname);
data = readtable(specname);
data = table2array(data);
database = data;
t=data(:,1); %示波器横坐标
x=data(:,2); %示波器纵坐标
 
%先把示波器的数据分两列导入工作区，横坐标命名为t，纵坐标命名为x
subplot(2,1,1);
plot(t,x); %画出原始信号
 
% %求FFT
% N=length(t); %计算信号长度
% Fs=1/(t(2)-t(1));%信号采样频率，计算前两个时间点的差的倒数
% 
% % 生成频率轴（以Hz为单位），用于后续绘图展示
% f = (0:N-1)*(Fs/N);
 
% 假设 t 和 x 已经定义，其中 t 是时间向量，x 是对应的信号值
 
% 采样频率，假设我们知道这个值
Fs = 1 / (t(2) - t(1)); % 如果 t 是等间隔的，可以通过两个相邻的时间差来计算
 
% 汉明窗大小等于信号长度
window = hamming(length(x));
 
% 对信号应用汉明窗
x_windowed = x .* window;
 
% 执行FFT并取绝对值得到单边频谱
X_fft = fft(x_windowed);
 
% 计算功率谱密度（PSD）
Pxx = (abs(X_fft).^2) / length(x); % 功率谱
 
% 转换到单边频谱，并考虑能量归一化
n = length(x);
if mod(n, 2) == 0
    Pxx = Pxx(1:n/2+1);
    Pxx(2:end-1) = 2*Pxx(2:end-1);
else
    Pxx = Pxx(1:(n+1)/2);
    Pxx(2:end-1) = 2*Pxx(2:end-1);
end
 
% 将功率谱转换为 dBm 单位
% 注意：这里假设输入信号是以伏特为单位的电压信号，
% 并且负载电阻是50欧姆。如果情况不同，请调整公式。
R = 50; % 负载电阻，通常为50欧姆
Pxx_dBm = 10*log10(Pxx * (1e3 / R)); % 将瓦特转换为毫瓦再转为dBm
 
% 频率轴
f = linspace(0, Fs/2, length(Pxx));
 
% 绘制频谱图
subplot(2,1,2);
plot(f, Pxx_dBm);
xlabel('Frequency (Hz)');
ylabel('Power/Frequency (dBm/Hz)');
title('Single-Sided Amplitude Spectrum of x(t)');
grid on;