Wavelet 是蛮复杂的一个信号分析模式,我们主要利用这个做EEG信号的时频分析

时频分析上述讲过可以通过短时间傅里叶变换来做,但是会造谷频率泄露

所以小波分析很好地解决这个问题,即可获得信号频率而且可以知道频率所在的位置

这篇主要分两个部分,

第一个是利用小波分析来画时频图

第二个是利用希尔伯特变换求瞬时频率 相位 还有振幅

@小波变换画时频图

小波有很多种,也可以参考窗函数,其实就是各种小波基

目前效果最高的应属morlet小波,其中心频率越大,体现出的效果越好

 

三个matlab函数

COEFS = cwt(S,SCALES,’wname’)
    说明:该函数能实现连续小波变换,其中S为输入信号,SCALES为尺度,wname为小波名称。

    FREQ = centfrq(‘wname’)
    说明:该函数能求出以wname命名的母小波的中心频率。

    F = scal2frq(A,’wname’,DELTA)
    说明:该函数能将尺度转换为实际频率,其中A为尺度,wname为小波名称,DELTA为采样周期。

 

   设a为尺度,fs为采样频率,Fc为小波中心频率,则a对应的实际频率Fa为
                     
                      Fa=Fc×fs/a                                   

    显然,根据采样定理,为使小波尺度图的频率范围为(0,fs/2),尺度范围应为(2*Fc,inf),其中inf表示为无穷大。在实际应用中,只需取尺度足够大即可。

 

代码实现:

clear;
clc;
fs=1024;                          %采样频率
f1=50;
f2=100;
t=0:1/fs:1;

s=sin(2*pi*78*t)+sin(2*pi*f2*t)+cos(2*pi*200*t);%两个不同频率正弦信号合成的仿真信号
%%%%%%%%%%%%%%%%%小波时频图绘制%%%%%%%%%%%%%%%%%%
wavename=’cmor3-3′;
totalscal=256;                    %尺度序列的长度,即scal的长度
wcf=centfrq(wavename);            %小波的中心频率
cparam=2*wcf*totalscal;           %为得到合适的尺度所求出的参数
a=totalscal:-1:1;
scal=cparam./a;                   %得到各个尺度,以使转换得到频率序列为等差序列
coefs=cwt(s,scal,wavename);       %得到小波系数
f=scal2frq(scal,wavename,1/fs);   %将尺度转换为频率
imagesc(t,f,abs(coefs));          %绘制色谱图
colorbar;
xlabel(‘时间 t/s’);
ylabel(‘频率 f/Hz’);
title(‘小波时频图’);

untitled

幅频分析:

主要利用FFT来运算

fs=100%抽样频率等于100, 单位时间内抽样100个点

t=0:0.01:10 %10秒的总信号长度,注意取值长度要和抽样保持一致
N=10*fs+1;

y=sin(2*pi*34*t)+sin(2*pi*12*t);%仿真信号

f_x=0:fs/(N-1):fs;%确定X轴是,一般0:fs(采样频率),取值间隔恰好是抽样频率除以总点数减去1.

z=fft(y);

figure
plot(t,y);
figure
plot(f_x,abs(z));

 

untitle11d

 

根据娜斯奎特定理,fft能计算的频率最高不超过50hz,因为采样频率我们采用的是100hz, 意味着 超过50hz的频率无法被争取的检测出,所以50hz以后的频率就不要参考了

 

 

幅频分析 指定抽样点数

fs=100%抽样频率等于1000, 单位时间内抽样1000个点

N=128;% total 128 datapoint

t=0:0.01:128/100;
y=sin(2*pi*33*t)+sin(2*pi*12*t);%仿真信号

f_x=0:fs/N:fs;

z=fft(y);

figure
subplot(2,1,1)
plot(t,y);

hold on
subplot(2,1,2)
plot(f_x,abs(z));

 

untitl222ed