本帖最后由 happy 于 2010-10-19 16:48 编辑
本程序来自网络,版权归原作者所有,由于原作者无法考究,故此处没有标明
如原作者看到本贴,请和管理员联系添加相关说明或做其他处理
- %%======================%%
- %%小波包分析程序
- %%=======================%%
- clc;clear all;close all;
- % %-----------------输入实际信号---------------
- N=8192;
- fs=20480;
- [DATAfile DATApath]=uigetfile('.txt','输入信号');%显示一个取文件的窗口
- FILENAME=[DATApath,DATAfile]
- DATA=load(FILENAME);
- Data=DATA(1:N);
- DataLen=length(Data);%数据长度
- t=[0:(DataLen-1)]'/fs;%时间数组
- % -------------------小波包分解---------------
- n=3;%分解层数
- rr=wpdec(Data,n,'db15');%小波包分解
- % ---------对特征频带波形进行分析-----
- for i=1:2^n
- Datar(i,:)=wprcoef(rr,[n,i-1]);%求取每个频带的小波包系数
- end
- Datar=[Datar(1,:);Datar(2,:);Datar(4,:);Datar(3,:);Datar(7,:);Datar(8,:);Datar(6,:);Datar(5,:)];%数据调整
- f=(0:N/2-1)*fs/N;
- for i=1:2^n
- DataFFT(i,:)=abs(fft(Datar(i,:)))/(N/2); %%fft
- end
- % %-----------小波包分解计算能量-------
- for i=0:(2^n-1)
- rcfs=wprcoef(rr,[n i]);
- r(i+1,:)=rcfs;
- end
- r_check=[r(1,:);r(2,:);r(4,:);r(3,:);r(7,:);r(8,:);r(6,:);r(5,:)];
- for i=1:2^n
- xx(i)=sum(r_check(i,:).^2)/N;%求能量
- end
- temp=sum(xx);
- for i=1:2^n
- xx(i)=xx(i)/temp;%%能量归一化
- end
- %---------画小波分频带图---------
- % figure()
- % for i=1:2^n/2
- % subplot(4,1,i);
- % plot(t,Datar(i,:));ylabel(i);axis([min(t) max(t) min(Datar(i,:)) max(Datar(i,:))]);
- % if i==1
- % title('小波包分解')
- % end
- % end
- % xlabel('time');
- % grid on
- % figure()
- % for i=2^n/2+1:2^n
- % subplot(4,1,i-2^n/2);
- % plot(t,Datar(i,:));ylabel(i);axis([min(t) max(t) min(Datar(i,:)) max(Datar(i,:))]);
- % if i== 2^n/2+1
- % title('小波包分解')
- % end
- % end
- % xlabel('time');
- % grid on
- %%-----------对每个频带求取频谱--------
- figure()
- set(gcf,'Name','ff1')
- for i=1:2^n/2
- subplot(4,1,i);
- plot(f,DataFFT(i,1:N/2));ylabel(i);axis([0 10240 0 max(DataFFT(i,1:N/2))]);grid on
- if i==1
- title('小波包各频带频谱图')
- end
- end
- xlabel('frequency')
- figure()
- set(gcf,'Name','ff2')
- for i=2^n/2+1:2^n
- subplot(4,1,i-2^n/2);
- plot(f,DataFFT(i,1:N/2));ylabel(i);axis([0 10240 0 max(DataFFT(i,1:N/2))]);grid on
- if i==2^n/2+1
- title('小波包各频带频谱图')
- end
- end
- xlabel('frequency')
- %-------------画出能量谱-----------
- figure()
- set(gcf,'Name','能量谱')
- width=0.5;
- bar(xx,width);
- colormap(cool);
- title('小波包能量谱图');
- xlabel('频 带');
- ylabel('能量归一化');
- axis([0 9 0 1]);
- grid on
复制代码 |