Matlab信号处理:短时傅里叶变换

ops/2024/11/19 21:25:46/

短时傅里叶变换(简称STFT)是傅里叶变换在时频域的扩展,它是为分析频域随时间变化的非平稳信号。本文模拟一个啁啾信号(一个线性调频的信号),借助matlab的短时傅里叶变换函数stft,分析其时频特性并绘制时频图,同时复现了matlab的stft函数的计算过程,供大家参考。

1.具体案例

啁啾信号一般是指调频信号,即随着时间的变化,信号的主频在不断发生变化(或增加、或减少)。当这种信号当做音频输出时,听起来会像鸟的唧唧声,所以叫啁啾。

利用matlab的chirp函数生成了一个1s内从100Hz到5000Hz的线性调频信号(啁啾信号),采样频率为24000Hz,具体如下:

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

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

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

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

在上述原理基础上,手动复现stft的计算过程,获得的时频图如下:

该图与stft获得的时频图是一致的,计算二者的平均偏差结果。从定量结果来看,复现过程结果与函数计算结果一致,验证了复现过程的正确性。 

2.具体代码

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

输出:freq表示幅值谱的横轴,向量形式;P1表示幅值谱的纵轴,矩阵形式,单个信号的采样索引X信号数;Theta表示相位谱的纵轴,矩阵形式,单个信号的采样索引X信号数。

主函数main.m代码:

%% 信号处理——短时傅里叶变换clc
clear all
close all% 定义时间向量和啁啾信号参数
fs = 24000; % 采样率
T = 1/fs:1/fs:1; % 1秒长的时间向量
f0 = 100; % 起始频率100Hz
f1 = 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");figure
mesh(t,f,abs(s));colorbar;
xlabel('Time/s');ylabel('Frequency/Hz');zlabel('Amplitude/g');axis tight;
title("matlab的stft函数获得时频图")figure
imagesc(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_countsindex=(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;   % 设置频率刻度  横轴HzP2 = abs(Y);P1 = P2(1:FFT_Length/2+1)';P1(:,3)=stft_abs(:,i);y_matrix__man(:,i)=P1(:,1);
endfigure;
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'));

幅值谱和相位谱计算函数frequ_am_phase.m代码:

function [freq,P1,Theta]=frequ_am_phase(y,fs,tol)%% 绘制信号频域的幅值谱和相位谱
%% 参数解释: 
%     y: 表示输入信号,它可以为一个矩阵,行X列,具体为单个信号的采样索引X信号数
%        比如y的大小为8192X12,表示一个有12个信号的数据矩阵,每个信号长度为8192
%        注意,如果仅有一个信号,则y应该是一个列向量
%        同时,y的行数尽量为偶数,奇数的话会引起程序索引的警告
%     fs:表示采样频率
%     tol:相位阈值参数
%     freq:表示幅值谱的横轴
%     P1:表示幅值谱的纵轴
%     Theta:表示相位谱的纵轴if nargin==2tol=1e-6;  %计算误差的默认阈值
endL=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;
end
Theta=theta(1:L/2+1,:);
end

http://www.ppmy.cn/ops/135067.html

相关文章

HTML5实现剪刀石头布小游戏(附源码)

文章目录 1.设计来源1.1 主界面1.2 皮肤风格1.2 游戏中界面 2.效果和源码源码下载万套模板&#xff0c;程序开发&#xff0c;在线开发&#xff0c;在线沟通 作者&#xff1a;xcLeigh 文章地址&#xff1a;https://blog.csdn.net/weixin_43151418/article/details/143798520 HTM…

python学习_2.去除字符strip方法

.strip() 是 Python 字符串的一个方法&#xff0c;用于去除字符串首尾的空白字符&#xff08;包括空格、制表符、换行符等&#xff09;。这个方法非常有用&#xff0c;特别是在处理从文件或用户输入中读取的字符串时&#xff0c;可以确保字符串没有多余的空白字符。 示例 假设…

Vue之el-date-picker日期选择器标签—选择日期范围,数据格式:yyyy-MM-dd HH:mm:ss,设置默认时间:HH:mm:ss

Vue之el-date-picker日期选择器标签—选择日期范围&#xff0c;数据格式:yyyy-MM-dd HH:mm:ss&#xff0c;设置默认时间:HH:mm:ss 需求是选择日期范围&#xff0c;即只能选择日期&#xff0c;但是想要的数据格式带有时间&#xff1a;yyyy-MM-dd HH:mm:ss&#xff0c;而且开始时…

shell命令统计文件行数之和

你可以使用以下 shell 命令来统计每个 .txt 文件的行数,并将其加和在一起: find . -name "*.txt" -not -name "*.json" -exec wc -l {} + | awk {sum += $1} END {print sum} 解释: find . -name "*.txt" -not -name "*.json": f…

PostgreSQL 函数与存储过程及调用

PostgreSQL 随着云服务的盛行&#xff0c;越发被广泛的应用&#xff0c;免费开源且有丰富的特性支持&#xff0c;加上性能也很不错&#xff0c;因而备受青睐。PostgreSQL 的函数与存储过程区别并不太大&#xff0c;不像某些数据库的函数与存储过程必须是无副作用或有副作用&…

自建k8s集群,利用开源的GitLab、Jenkins和Harbor实现CI/CD和DevOps的过程回顾

使用自己部署的Kubernetes集群&#xff0c;结合GitLab、Jenkins和Harbor实现CI/CD和DevOps的过程大致如下&#xff1a; 1.代码管理&#xff08;GitLab&#xff09;&#xff1a; - 开发者在GitLab上创建代码仓库&#xff0c;编写代码并提交变更。 - 每次代码提交都会触发GitL…

创建游戏云存档功能的完整指南

✅作者简介&#xff1a;2022年博客新星 第八。热爱国学的Java后端开发者&#xff0c;修心和技术同步精进。 &#x1f34e;个人主页&#xff1a;Java Fans的博客 &#x1f34a;个人信条&#xff1a;不迁怒&#xff0c;不贰过。小知识&#xff0c;大智慧。 &#x1f49e;当前专栏…

大模型时代,呼叫中心部门如何建设一套呼出机器人系统?

大模型时代&#xff0c;呼叫中心部门如何建设一套呼出机器人系统&#xff1f; 作者&#xff1a;开源呼叫中心系统 FreeIPCC&#xff0c;Github地址&#xff1a;https://github.com/lihaiya/freeipcc 在大模型时代&#xff0c;呼叫中心部门建设一套呼出机器人系统需要综合考虑技…