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

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

其中一段代码是：

%dalt

fs1=fft(dalt,1024);%快速付氏变换

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

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

其中173.61是什么意思

完整代码如下：

load singal.mat

s=sig;

fs=1024;%采样频率

dalt=1/fs;%采样间隔

t=0:length(s)-1*dalt;

figure(1);

subplot(311);plot(s);title('原始信号');

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

fs=fft(s,1024);%快速付氏变换

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

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

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

l=wmaxlev(s,'db4');%取db4分解变量s的最大尺度l 实际分解尺度小于l

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

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

[LD,HD,LR,HR]=wfilters('db4');

figure(2);

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

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

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

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

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

C7=appcoef(C,L,'db4',7);%尺度128

D7=detcoef(C,L,7);%细节系数

D6=detcoef(C,L,6);%尺度64

D5=detcoef(C,L,5);%尺度32

D4=detcoef(C,L,4);%尺度16

D3=detcoef(C,L,3);%尺度8

D2=detcoef(C,L,2);%尺度4

D1=detcoef(C,L,1);%尺度2

figure(3);

subplot(711);plot(D7);ylabel('D7');title('细节系数');

subplot(712);plot(D6);ylabel('D6');

subplot(713);plot(D5);ylabel('D5');

subplot(714);plot(D4);ylabel('D4');

subplot(715);plot(D3);ylabel('D3');

subplot(716);plot(D2);ylabel('D2');

subplot(717);plot(D1);ylabel('D1');

[C,L]=wavedec(s,6,'db4');%近似系数 尺度64

C6=appcoef(C,L,'db4',6);

[C,L]=wavedec(s,5,'db4');%尺度32

C5=appcoef(C,L,'db4',5);

[C,L]=wavedec(s,4,'db4');%尺度16

C4=appcoef(C,L,'db4',4);

[C,L]=wavedec(s,3,'db4');%尺度8

C3=appcoef(C,L,'db4',3);

[C,L]=wavedec(s,2,'db4');%尺度4

C2=appcoef(C,L,'db4',2);

[C,L]=wavedec(s,1,'db4');%尺度2

C1=appcoef(C,L,'db4',1);

figure(4);

subplot(711);plot(C7);ylabel('C7');title('近似系数');

subplot(712);plot(C6);ylabel('C6');

subplot(713);plot(C5);ylabel('C5');

subplot(714);plot(C4);ylabel('C4');

subplot(715);plot(C3);ylabel('C3');

subplot(716);plot(C2);ylabel('C2');

subplot(717);plot(C1);ylabel('C1');

%带宽 0.5 Hz to 85 Hz;采样频率fs=173.61Hz

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

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

[C,L]=wavedec(s,7,'db4');%尺度128~64

C7=appcoef(C,L,'db4',7);

D7=detcoef(C,L,7);

D6=detcoef(C,L,6);%尺度64

D5=detcoef(C,L,5);%尺度32

D4=detcoef(C,L,4);%尺度16

D3=detcoef(C,L,3);%尺度8

D2=detcoef(C,L,2);%尺度4

D1=detcoef(C,L,1);%尺度2

[C,L]=wavedec(s,7,'db4');

SRD7=wrcoef('d',C,L,'db4',7);

SRD6=wrcoef('d',C,L,'db4',6);

SRD5=wrcoef('d',C,L,'db4',5);

SRD4=wrcoef('d',C,L,'db4',4);

SRD3=wrcoef('d',C,L,'db4',3);

SRD2=wrcoef('d',C,L,'db4',2);

SRD1=wrcoef('d',C,L,'db4',1);

dalt=[SRD5];%δ-wave(1~4Hz)

thalt=[SRD4];%θ-wave(4~8Hz)

alpha=[SRD3];%α-wave(8~13Hz)

belta=[SRD2];%β-wave(14~30Hz)

%dalt

fs1=fft(dalt,1024);%快速付氏变换

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

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

%thalt

fs2=fft(thalt,1024);%快速付氏变换

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

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

%alpha

fs3=fft(alpha,1024);%快速付氏变换

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

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

%belta

fs4=fft(belta,1024);%快速付氏变换

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

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

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

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

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

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

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

figure(6);

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

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

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

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

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

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