chirp信号-信号处理——短时傅里叶变换

内容概括

在故障诊断中，信号频域分析的基本手段公众号基本涵盖，有遗漏的相关方法和技术，知道的小伙伴可以后台私信我，我继续添加！

这次给大家带来的是短时傅里叶变换（Short-Time ，简称STFT），它是傅里叶变换在时频域的扩展，它是为了分析频域随时间变化的非平稳信号。

本文的流程如下图所示，模拟了一个啁啾信号（一个线性调频的信号），并借助的短时傅里叶变换函数stft，分析了其时频特性，绘制了时频图。同时，复现了的stft函数的计算过程，供大家参考。

该内容主要参考了一些资料：

代码采用了 2024a进行运行，欢迎大家测试和提出问题！

具体案例

啁啾(zhōu jiū，英文中为chirp)信号一般是指调频信号，即随着时间的变化，信号的主频在不断发生变化（或增加、或减少）。当这种信号当做音频输出时，听起来会像鸟的唧唧声，所以叫啁啾。

利用的chirp函数生成了一个1s内从100Hz到的线性调频信号（啁啾信号），采样频率为，具体如下：

从啁啾信号的局部细节图能发现，随着时间的增长，信号波形越来越密集，即信号的频率逐渐增大。在啁啾信号频谱中，低频占据信号的主频，无法发现信号频率的时变特征。

采用stft函数，设置矩形窗，窗长256，滑动步长设置为128，FFT分析的信号长度与窗长一致chirp信号，均为256，设置stft函数的为,即仅分析信号的傅里叶变换的单边频谱，获得的短时傅里叶变换结果如下图。

上图分别为短时傅里叶变换的不同画图方式，分别为mesh函数和函数。从上图中能发现，啁啾信号的频率从100Hz到线性增长，这表明短时傅里叶变换能较好地分析该啁啾信号，对比FFT的结果，表明了短时傅里叶变换在信号时频分析方面的优势。

stft函数的计算流程如下，它利用窗函数，将原始信号划分为了若干个子信号，然后分别采用傅里叶变换，获得其对应的频谱，最后将这些频谱按照不同的时间先后顺序进行拼接，绘制了时频图。

因此，在上述原理基础上，手动复现了stft的计算过程，获得的时频图如下：

该图与stft获得的时频图是一致的，计算二者的平均偏差结果：（平均偏差为两数差的绝对值的平均数）

从定量结果来看，复现过程结果与函数计算结果一致，验证了复现过程的正确性。

注意：值得注意是，stft函数计算的是子信号的FFT的单边频谱的绝对值，它并未进行信号长度的标准化，就是描述信号的频频幅值，不是能量也不是功率。

具体代码

代码主要分为两个：

1、main.m（主函数）

2、.m（幅值谱和相位谱计算函数）

说明：

.m函数

调用形式：

[freq,P1,Theta]=(y,fs,tol)

输入：

信号y，矩阵，行X列=单个信号的采样索引X信号数，比如信号的大小为，表示一个有12个信号的数据矩阵，每个信号长度为8192。注意，如果仅有一个信号，则y应该是一个列向量。同时，y的行数尽量为偶数，奇数的话会引起程序索引的警告。

fs:采样频率，标量

tol:相位阈值参数，标量（可缺省chirp信号，下面讲解）

输出

freq:表示幅值谱的横轴，向量，HzX1

P1:表示幅值谱的纵轴，矩阵，单个信号的采样索引X信号数，类比信号y

Theta:表示相位谱的纵轴，矩阵，单个信号的采样索引X信号数，类比信号y

main.m（主函数）

%% 信号处理——短时傅里叶变换%% 作者：冷漠%% 时间：2024年9月28日%% 关注公众号 ："故障诊断与寿命预测工具箱",每天进步一点点clcclear allclose all
% 定义时间向量和啁啾信号参数fs = 24000; % 采样率T = 1/fs:1/fs:1; % 1秒长的时间向量f0 = 100; % 起始频率100Hzf1 = 5000; % 结束频率5000Hz
% 生成啁啾信号y = chirp(T, f0, 1, f1);
% 播放啁啾信号sound(y, fs);
filename = 'chirpSignal.mp3';% 使用audiowrite函数保存为MP3格式audiowrite(filename, y, fs, 'Quality', 95);
%[freq1,P1,~]=frequ_am_phase(y',fs);figure;subplot(311);plot(T,y,'b-');xlabel('Time/s');ylabel('Amplitude/g');axis tight;title("啁啾信号的时域波形")subplot(312);plot(T(1:5000),y(1:5000),'b-');xlabel('Time/s');ylabel('Amplitude/g');axis tight;title("时域波形的局部细节")subplot(313);plot(freq1,P1,'b-');xlabel('Frequency/Hz');ylabel('Amplitude/g');title("啁啾信号的频域波形")

[s,f,t] = stft(y,fs,Window=rectwin(256),OverlapLength=128,FFTLength=256,FrequencyRange="onesided");

figuremesh(t,f,abs(s));colorbar;xlabel('Time/s');ylabel('Frequency/Hz');zlabel('Amplitude/g');axis tight;title("matlab的stft函数获得时频图")
figureimagesc(t,f,abs(s));colorbar;set(gca, 'YDir', 'normal');xlabel('Time/s');ylabel('Frequency/Hz');axis tight;title("matlab的stft函数获得时频图")
%%  复现短时傅里叶变换结果rect_win=256;  %矩形窗口的长度FFT_Length=256;  %FFT分析的长度  overlap_len=128; %分割信号重叠的长度seg_counts=floor((length(y)-overlap_len)/(rect_win-overlap_len));stft_abs=(abs(s));for i=1:seg_counts    index=(i-1)*(rect_win-overlap_len)+1:(i-1)*(rect_win-overlap_len)+rect_win;    sub_y=y(index);    %复现过程    Y=fft(sub_y,FFT_Length);          % FFT 快速傅里叶变换    freq=(0:FFT_Length/2)*fs/FFT_Length;   % 设置频率刻度  横轴Hz    P2 = abs(Y);    P1 = P2(1:FFT_Length/2+1)';    P1(:,3)=stft_abs(:,i);    y_matrix__man(:,i)=P1(:,1);end
figure;imagesc(t,freq,y_matrix__man);colorbar;set(gca, 'YDir', 'normal');xlabel('Time/s');ylabel('Frequency/Hz');axis tight;title("手动计算绘制的短时傅里叶变换时频图")
fprintf('matlab的stft函数获得时频图与手写复现之间的差为%f n',sum(stft_abs-y_matrix__man,'all'));

.m（幅值谱和相位谱计算函数）

function [freq,P1,Theta]=frequ_am_phase(y,fs,tol)%% 作者：冷漠%% 时间：2024年6月01日%% 关注公众号 ："故障诊断与寿命预测工具箱",每天进步一点点%% 绘制信号频域的幅值谱和相位谱%% 参数解释： %     y: 表示输入信号，它可以为一个矩阵，行X列,具体为单个信号的采样索引X信号数%        比如y的大小为8192X12，表示一个有12个信号的数据矩阵，每个信号长度为8192%        注意，如果仅有一个信号，则y应该是一个列向量%        同时，y的行数尽量为偶数，奇数的话会引起程序索引的警告%     fs:表示采样频率%     tol:相位阈值参数%     freq:表示幅值谱的横轴%     P1:表示幅值谱的纵轴%     Theta:表示相位谱的纵轴
if nargin==2    tol=1e-6;  %计算误差的默认阈值end
L=size(y,1);         % 信号长度% Y=fft(y,2^nextpow2(L));          % FFT 快速傅里叶变换Y=fft(y,L);          % FFT 快速傅里叶变换freq=(0:L/2)*fs/L;   % 设置频率刻度  横轴Hz%幅值谱P2 = abs(Y/L);P1 = P2(1:L/2+1,:);P1(2:end-1,:) = 2*P1(2:end-1,:);  %纵轴 幅值
%相位谱P2(2:end-1,:)=2*P2(2:end-1,:);for i=1:size(Y,2)    Y(P2(:,i)<tol,i) = 0;    theta(:,i) = angle(Y(:,i))/pi;endTheta=theta(1:L/2+1,:);end

细节说明

1、关于啁啾信号的具体公式，不太了解，后续有机会单独出一期讲讲。

2、从复现stft函数来看，它绘制的子信号的FFT的单边频谱的绝对值。

3、本文采用了矩形窗函数，确保了stft和手动计算的一致性。

总结

本文模拟了一个啁啾信号，他是一个线性调频信号，并利用短时傅里叶变换分析了它频率的时变特性。同时，复现了的stft函数的计算过程，结果表明复现代码的正确性。

相关文件

1、上述源码

.m(主函数)；

.m（幅值谱和相位谱计算函数）

啁啾信号的mp3文件（.mp3）

上述相关资料，请关注本公众号，回复：信号处理-短时傅里叶变换，获取网盘链接。

关注公众号“故障诊断与寿命预测工具箱”，每天进步一点点。

欢迎大家提意见哈，公众号的定位是一个工具箱，冷漠会尽量满足大家需求！

你们点点“分享”，给我充点儿电吧~

admin会员

打赏 收藏 海报 链接