matlab中morl小波,小波变换后的各频率分量的功率谱，

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

下面程序是用来提取脑电波信号的&＃xff0c;利用小波提取四种频率分量的波α、β、δ、θ&＃xff0c;现在不明白如何求频率分量的功率谱&＃xff1f;

其中一段代码是&＃xff1a;

%dalt

fs1&＃61;fft(dalt,1024);%快速付氏变换

pp1&＃61;fs1.*conj(fs1)/1024;%计算功率谱

ff1&＃61;((0:511)/1024)*173.61;%计算各点对应的频率值

其中173.61是什么意思

完整代码如下&＃xff1a;

load singal.mat

s&＃61;sig;

fs&＃61;1024;%采样频率

dalt&＃61;1/fs;%采样间隔

t&＃61;0:length(s)-1*dalt;

figure(1);

subplot(311);plot(s);title(&＃39;原始信号&＃39;);

axis([-10 4000 -200 200]);

fs&＃61;fft(s,1024);%快速付氏变换

pp&＃61;fs.*conj(fs)/1024;%计算功率谱

ff&＃61;(0:511)/1024/dalt;%计算各点对应的频率值

subplot(312);plot(ff,pp(1:512));ylabel(&＃39;功率谱密度&＃39;);xlabel(&＃39;频率&＃39;);title(&＃39;信号功率谱&＃39;);axis([0 55 -10 110000]);

l&＃61;wmaxlev(s,&＃39;db4&＃39;);%取db4分解变量s的最大尺度l 实际分解尺度小于l

sd&＃61;wden(s,&＃39;minimaxi&＃39;,&＃39;s&＃39;,&＃39;mln&＃39;,real(l),&＃39;db4&＃39;);%除噪 dbN 对称 minimaxi用极大极小原理进行阈值选择 &＃39;s&＃39;为软阈值 &＃39;mln&＃39;在不同层估计噪声层并调整阈值

subplot(313);plot(sd);xlabel(&＃39;滤除噪声信号&＃39;);axis([-10 4000 -200 200]);

[LD,HD,LR,HR]&＃61;wfilters(&＃39;db4&＃39;);

figure(2);

subplot(221);stem(LD);title(&＃39;LD 低通分解滤波器&＃39;);grid;

subplot(222);stem(HD);title(&＃39;HD 低通重构滤波器&＃39;);grid;

subplot(223);stem(LR);title(&＃39;LR 高通分解滤波器&＃39;);grid;

subplot(224);stem(HR);title(&＃39;HR 高通重构滤波器&＃39;);grid;

[C,L]&＃61;wavedec(s,7,&＃39;db4&＃39;); %用小波db4对信号进行多尺度分解(7层)

C7&＃61;appcoef(C,L,&＃39;db4&＃39;,7);%尺度128

D7&＃61;detcoef(C,L,7);%细节系数

D6&＃61;detcoef(C,L,6);%尺度64

D5&＃61;detcoef(C,L,5);%尺度32

D4&＃61;detcoef(C,L,4);%尺度16

D3&＃61;detcoef(C,L,3);%尺度8

D2&＃61;detcoef(C,L,2);%尺度4

D1&＃61;detcoef(C,L,1);%尺度2

figure(3);

subplot(711);plot(D7);ylabel(&＃39;D7&＃39;);title(&＃39;细节系数&＃39;);

subplot(712);plot(D6);ylabel(&＃39;D6&＃39;);

subplot(713);plot(D5);ylabel(&＃39;D5&＃39;);

subplot(714);plot(D4);ylabel(&＃39;D4&＃39;);

subplot(715);plot(D3);ylabel(&＃39;D3&＃39;);

subplot(716);plot(D2);ylabel(&＃39;D2&＃39;);

subplot(717);plot(D1);ylabel(&＃39;D1&＃39;);

[C,L]&＃61;wavedec(s,6,&＃39;db4&＃39;);%近似系数 尺度64

C6&＃61;appcoef(C,L,&＃39;db4&＃39;,6);

[C,L]&＃61;wavedec(s,5,&＃39;db4&＃39;);%尺度32

C5&＃61;appcoef(C,L,&＃39;db4&＃39;,5);

[C,L]&＃61;wavedec(s,4,&＃39;db4&＃39;);%尺度16

C4&＃61;appcoef(C,L,&＃39;db4&＃39;,4);

[C,L]&＃61;wavedec(s,3,&＃39;db4&＃39;);%尺度8

C3&＃61;appcoef(C,L,&＃39;db4&＃39;,3);

[C,L]&＃61;wavedec(s,2,&＃39;db4&＃39;);%尺度4

C2&＃61;appcoef(C,L,&＃39;db4&＃39;,2);

[C,L]&＃61;wavedec(s,1,&＃39;db4&＃39;);%尺度2

C1&＃61;appcoef(C,L,&＃39;db4&＃39;,1);

figure(4);

subplot(711);plot(C7);ylabel(&＃39;C7&＃39;);title(&＃39;近似系数&＃39;);

subplot(712);plot(C6);ylabel(&＃39;C6&＃39;);

subplot(713);plot(C5);ylabel(&＃39;C5&＃39;);

subplot(714);plot(C4);ylabel(&＃39;C4&＃39;);

subplot(715);plot(C3);ylabel(&＃39;C3&＃39;);

subplot(716);plot(C2);ylabel(&＃39;C2&＃39;);

subplot(717);plot(C1);ylabel(&＃39;C1&＃39;);

%带宽 0.5 Hz to 85 Hz;采样频率fs&＃61;173.61Hz

%δ-wave(1~4Hz);θ-wave(4~8Hz);α-wave(8~13Hz);β-wave(14~30Hz);

%******************************

[C,L]&＃61;wavedec(s,7,&＃39;db4&＃39;);%尺度128~64

C7&＃61;appcoef(C,L,&＃39;db4&＃39;,7);

D7&＃61;detcoef(C,L,7);

D6&＃61;detcoef(C,L,6);%尺度64

D5&＃61;detcoef(C,L,5);%尺度32

D4&＃61;detcoef(C,L,4);%尺度16

D3&＃61;detcoef(C,L,3);%尺度8

D2&＃61;detcoef(C,L,2);%尺度4

D1&＃61;detcoef(C,L,1);%尺度2

[C,L]&＃61;wavedec(s,7,&＃39;db4&＃39;);

SRD7&＃61;wrcoef(&＃39;d&＃39;,C,L,&＃39;db4&＃39;,7);

SRD6&＃61;wrcoef(&＃39;d&＃39;,C,L,&＃39;db4&＃39;,6);

SRD5&＃61;wrcoef(&＃39;d&＃39;,C,L,&＃39;db4&＃39;,5);

SRD4&＃61;wrcoef(&＃39;d&＃39;,C,L,&＃39;db4&＃39;,4);

SRD3&＃61;wrcoef(&＃39;d&＃39;,C,L,&＃39;db4&＃39;,3);

SRD2&＃61;wrcoef(&＃39;d&＃39;,C,L,&＃39;db4&＃39;,2);

SRD1&＃61;wrcoef(&＃39;d&＃39;,C,L,&＃39;db4&＃39;,1);

dalt&＃61;[SRD5];%δ-wave(1~4Hz)

thalt&＃61;[SRD4];%θ-wave(4~8Hz)

alpha&＃61;[SRD3];%α-wave(8~13Hz)

belta&＃61;[SRD2];%β-wave(14~30Hz)

%dalt

fs1&＃61;fft(dalt,1024);%快速付氏变换

pp1&＃61;fs1.*conj(fs1)/1024;%计算功率谱

ff1&＃61;((0:511)/1024)*173.61;%计算各点对应的频率值

%thalt

fs2&＃61;fft(thalt,1024);%快速付氏变换

pp2&＃61;fs2.*conj(fs2)/1024;%计算功率谱

ff2&＃61;((0:511)/1024)*173.61;%计算各点对应的频率值

%alpha

fs3&＃61;fft(alpha,1024);%快速付氏变换

pp3&＃61;fs3.*conj(fs3)/1024;%计算功率谱

ff3&＃61;((0:511)/1024)*173.61;%计算各点对应的频率值

%belta

fs4&＃61;fft(belta,1024);%快速付氏变换

pp4&＃61;fs4.*conj(fs4)/1024;%计算功率谱

ff4&＃61;((0:511)/1024)*173.61;%计算各点对应的频率值

figure(5);%特征波形提取

subplot(411);plot(ff1,pp1(1:512));ylabel(&＃39;功率谱密度&＃39;);title(&＃39;dalt信号功率谱&＃39;);%axis([0 30 -10 110000]);

subplot(412);plot(ff2,pp2(1:512));ylabel(&＃39;功率谱密度&＃39;);title(&＃39;thalt信号功率谱&＃39;);%axis([0 55 -10 110000]);

subplot(413);plot(ff3,pp3(1:512));ylabel(&＃39;功率谱密度&＃39;);title(&＃39;alpha信号功率谱&＃39;);%axis([0 55 -10 110000]);

subplot(414);plot(ff4,pp4(1:512));ylabel(&＃39;功率谱密度&＃39;);xlabel(&＃39;频率&＃39;);title(&＃39;belta信号功率谱&＃39;);%axis([0 55 -10 110000]);

figure(6);

subplot(311);xlabel(&＃39;时间 t/s&＃39;);ylabel(&＃39;频率 f/Hz&＃39;);zlabel(&＃39;幅值 A/uV&＃39;);

c &＃61; cwt(s,[64 32 16:-2:2],&＃39;morl&＃39;,&＃39;3Dlvl&＃39;); %morl只能进行连续小波变换 morl具有对称性&＃xff0c;不生成正交系&＃xff0c;不具有有限支集等性质

subplot(312);xlabel(&＃39;时间 t/s&＃39;);ylabel(&＃39;频率 f/Hz&＃39;);zlabel(&＃39;幅值 A/uV&＃39;);

c &＃61; cwt(s,[16 8 4 2 1],&＃39;morl&＃39;,&＃39;3Dlvl&＃39;); %morl进行连续小波变换 &＃39;lvl&＃39; By scale不同尺度对比

subplot(313);xlabel(&＃39;时间 t/s&＃39;);ylabel(&＃39;频率 f/Hz&＃39;);zlabel(&＃39;幅值 A/uV&＃39;);

c &＃61; cwt(s,[16 8 4 2 1],&＃39;mexh&＃39;,&＃39;3Dlvl&＃39;);%用mexihat进行连续小波变换 不同小波对比