信号的时频分析

时频分析是将时域信号转变为在时域和频域上的分布

1.小波

采用小波变换(Wavelet Transform, WT)对各通道信号进行时频分析。其原理如式(1)所示。

                                             (1)

式中α是尺度函数，它标定了小波变换的分析频率；τ是时移参数，它定义了分析时刻,ψ(t)是母小波。尺度函数与实际频率的关系如式(2)所示。

                       (2)

式中FC为母小波的中心频率，n为频率分度，k为归一化频率。根据实验需求，选择母小波为cmor3-3，它的中心频率FC=3，选定频率分度n=2000,即对采样率Fs=1000Hz的信号具有0.25Hz的频率分辨率。

参考文献：：：：：.胡广书,现代信号处理教程[M].第二版.2015,北京:清华大学出版社.421-424（强烈推荐，老教授人超级nice）

大概是这样的图

 

%%image the  time-frequency distribution with wavelet 
clear
load('sig.mat');
[bbw,abw]=cheby1(4,0.5,1.5/180,'high');
sig=filtfilt(bbw,abw,lfp4);  %高通滤波，去基线漂移
N=length(sig);

fs=1000;
f1=50;
f2=100;
t=1/fs:1/fs:N/fs;
F1=100;%设定最高频率值
F2=4;F3=12;%f2-f3为theta band 
F4=30; F5=80;%f4-f5 is gamma band

% figure
% plot(t, sig)
% 连续小波变换
wavename='cmor3-3';%母小波
totalscal=256*50;
Fc=centfrq(wavename); % 小波的中心频率
c=2*Fc*totalscal;
scals=c./(1:totalscal);
f=scal2frq(scals,wavename,1/fs); % 将尺度转换为频率
coefs=cwt(sig,scals,wavename); % 求连续小波系数
Nf=length(f);
len1= floor(F1*Nf/(fs/2));
%image the 0-100Hz time-frequency distribution 
figure
imagesc(t,f(1:len1),abs(coefs(1:len1,:)));
set(gca,'YDir','normal')
 colormap(jet); view(0,90);colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('小波时频图');
%image the 4-12Hz time-frequency distribution 
len2= floor(F2*Nf/(fs/2));
len3= floor(F3*Nf/(fs/2));
figure
imagesc(t,f(len2:len3),abs(coefs(len2:len3,:)));
set(gca,'YDir','normal')
 colormap(jet); view(0,90);colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('theta band 小波时频图');
%image the 30-80Hz time-frequency distribution 
len4= floor(F4*Nf/(fs/2));
len5= floor(F5*Nf/(fs/2));
figure
imagesc(t,f(len4:len5),abs(coefs(len4:len5,:)));
set(gca,'YDir','normal')
 colormap(jet); view(0,90);colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('gamma band 小波时频图');

 

2.STFT

本论文对经过预处理后多通道LFPs进行时频变换。对多通道LFPs中每一个通道的LFP信号进行应用短时傅里叶变换（Short-time Fourier Transform，STFT），获取多通道LFPs频谱的时空分布，研究LFPs在工作记忆过程中的能量在时间和频率上动态变化特性。设单通道LFP为信号x(t)，其STFT定义如式4.3所示：

                (4.3)

其中 是窗函数。选取时间窗口为1s，窗口移动步长为200ms，分别计算sig的动态时频分布，

clear 
 load('sig')

F=100;%视频图的最高频率
fs=1000;%采样频率，可从软件上设定
Ts=1/fs;%时间间隔为采样频率的倒数
NLen=length(sig);%信号长度
f = fs*(1/NLen:1/NLen:1);%频率范围
T1=-2500;%起始时间
T2=1500;%终止时间
t=(1/fs:NLen)/fs;
[bbw,abw]=cheby1(4,0.5,1.5/180,'high');
sig1=filtfilt(bbw,abw,sig);
h=window('hamming',255);
tfr = tfrstft(sig1,1:NLen,NLen,h);
len= floor(F*NLen/fs);
 figure;
 imagesc(t,f(1:len),abs(tfr(1:len,:)));
 axis xy;axis tight;
 colormap(jet); view(0,90);colorbar;
 ylabel('Freq (Hz)');
 xlabel('Time (s)');

还有其他的时频分析方法

各有优劣势

tfrpwv、tfrgabor等等