小波變換程序

1、小波濾波器構(gòu)造和消噪程序（2個）1重構(gòu)% mallet_wavelet.m% 此函數(shù)用于研究Mallet算法及濾波器設(shè)計(jì)% 此函數(shù)僅用于消噪a=pi/8; %角度賦初值b=pi/8;%低通重構(gòu)FIR濾波器h0(n)沖激響應(yīng)賦值h0=cos(a)*cos(b);h1=sin(a)*cos(b);h2=-sin(a)*sin(b);h3=cos(a)*sin(b);low_construct=h0,h1,h2,h3;L_fre=4; %濾波器長度low_decompose=low_construct(end:-1:1); %確定h0(-n)，低通分解濾波器for i_high=1:L_fre; %

2、確定h1(n)=(-1)n,高通重建濾波器if(mod(i_high,2)=0);coefficient=-1;elsecoefficient=1;endhigh_construct(1,i_high)=low_decompose(1,i_high)*coefficient;endhigh_decompose=high_construct(end:-1:1); %高通分解濾波器h1(-n)L_signal=100; %信號長度n=1:L_signal; %信號賦值f=10;t=0.001;y=10*cos(2*pi*50*n*t).*exp(-20*n*t);figure(1);plot(y)

3、;title('原信號');check1=sum(high_decompose); %h0(n)性質(zhì)校驗(yàn)check2=sum(low_decompose);check3=norm(high_decompose);check4=norm(low_decompose);l_fre=conv(y,low_decompose); %卷積l_fre_down=dyaddown(l_fre); %抽取,得低頻細(xì)節(jié)h_fre=conv(y,high_decompose);h_fre_down=dyaddown(h_fre); %信號高頻細(xì)節(jié)figure(2);subplot(2,1,1)pl

4、ot(l_fre_down);title('小波分解的低頻系數(shù)');subplot(2,1,2);plot(h_fre_down);title('小波分解的高頻系數(shù)');l_fre_pull=dyadup(l_fre_down); %0差值h_fre_pull=dyadup(h_fre_down);l_fre_denoise=conv(low_construct,l_fre_pull);h_fre_denoise=conv(high_construct,h_fre_pull);l_fre_keep=wkeep(l_fre_denoise,L_signal); %

5、取結(jié)果的中心部分,消除卷積影響h_fre_keep=wkeep(h_fre_denoise,L_signal);sig_denoise=l_fre_keep+h_fre_keep; %信號重構(gòu)compare=sig_denoise-y; %與原信號比較figure(3);subplot(3,1,1)plot(y);ylabel('y'); %原信號subplot(3,1,2);plot(sig_denoise);ylabel('sig_denoise'); %重構(gòu)信號subplot(3,1,3);plot(compare);ylabel('compare

6、'); %原信號與消噪后信號的比較2消噪% mallet_wavelet.m% 此函數(shù)用于研究Mallet算法及濾波器設(shè)計(jì)% 此函數(shù)用于消噪處理%角度賦值%此處賦值使濾波器系數(shù)恰為db9%分解的高頻系數(shù)采用db9較好,即它的消失矩較大%分解的有用信號小波高頻系數(shù)基本趨于零%對于噪聲信號高頻分解系數(shù)很大，便于閾值消噪處理l,h=wfilters('db10','d');low_construct=l;L_fre=20; %濾波器長度low_decompose=low_construct(end:-1:1); %確定h0(-n)，低通分解濾波器for i_h

7、igh=1:L_fre; %確定h1(n)=(-1)n,高通重建濾波器if(mod(i_high,2)=0);coefficient=-1;elsecoefficient=1;endhigh_construct(1,i_high)=low_decompose(1,i_high)*coefficient;endhigh_decompose=high_construct(end:-1:1); %高通分解濾波器h1(-n)L_signal=100; %信號長度n=1:L_signal; %原始信號賦值f=10;t=0.001;y=10*cos(2*pi*50*n*t).*exp(-30*n*t);z

8、ero1=zeros(1,60); %信號加噪聲信號產(chǎn)生zero2=zeros(1,30);noise=zero1,3*(randn(1,10)-0.5),zero2;y_noise=y+noise;figure(1);subplot(2,1,1);plot(y);title('原信號');subplot(2,1,2);plot(y_noise);title('受噪聲污染的信號');check1=sum(high_decompose); %h0(n),性質(zhì)校驗(yàn)check2=sum(low_decompose);check3=norm(high_decompose

9、);check4=norm(low_decompose);l_fre=conv(y_noise,low_decompose); %卷積l_fre_down=dyaddown(l_fre); %抽取,得低頻細(xì)節(jié)h_fre=conv(y_noise,high_decompose);h_fre_down=dyaddown(h_fre); %信號高頻細(xì)節(jié)figure(2);subplot(2,1,1)plot(l_fre_down);title('小波分解的低頻系數(shù)');subplot(2,1,2);plot(h_fre_down);title('小波分解的高頻系數(shù)')

10、;% 消噪處理for i_decrease=31:44;if abs(h_fre_down(1,i_decrease)>=0.000001h_fre_down(1,i_decrease)=(10-7);endendl_fre_pull=dyadup(l_fre_down); %0差值h_fre_pull=dyadup(h_fre_down);l_fre_denoise=conv(low_construct,l_fre_pull);h_fre_denoise=conv(high_construct,h_fre_pull);l_fre_keep=wkeep(l_fre_denoise,L_s

11、ignal); %取結(jié)果的中心部分,消除卷積影響h_fre_keep=wkeep(h_fre_denoise,L_signal);sig_denoise=l_fre_keep+h_fre_keep; %消噪后信號重構(gòu)%平滑處理for j=1:2for i=60:70;sig_denoise(i)=sig_denoise(i-2)+sig_denoise(i+2)/2;end;end;compare=sig_denoise-y; %與原信號比較figure(3);subplot(3,1,1)plot(y);ylabel('y'); %原信號subplot(3,1,2);plot(

12、sig_denoise);ylabel('sig_denoise'); %消噪后信號subplot(3,1,3);plot(compare);ylabel('compare'); %原信號與消噪后信號的比較小波譜分析mallat算法經(jīng)典程序clc;clear;% 1.正弦波定義f1=50; % 頻率1f2=100; % 頻率2fs=2*(f1+f2); % 采樣頻率Ts=1/fs; % 采樣間隔N=120; % 采樣點(diǎn)數(shù)n=1:N;y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts); % 正弦波混合figure(1)plot(y);tit

13、le('兩個正弦信號')figure(2)stem(abs(fft(y);title('兩信號頻譜')% 2.小波濾波器譜分析h=wfilters('db30','l'); % 低通g=wfilters('db30','h'); % 高通h=h,zeros(1,N-length(h); % 補(bǔ)零（圓周卷積，且增大分辨率變于觀察）g=g,zeros(1,N-length(g); % 補(bǔ)零（圓周卷積，且增大分辨率變于觀察）figure(3);stem(abs(fft(h);title('低通濾波

14、器圖')figure(4);stem(abs(fft(g);title('高通濾波器圖')% 3.MALLET分解算法(圓周卷積的快速傅里葉變換實(shí)現(xiàn))sig1=ifft(fft(y).*fft(h); % 低通(低頻分量)sig2=ifft(fft(y).*fft(g); % 高通(高頻分量)figure(5); % 信號圖subplot(2,1,1)plot(real(sig1);title('分解信號1')subplot(2,1,2)plot(real(sig2);title('分解信號2')figure(6); % 頻譜圖subpl

15、ot(2,1,1)stem(abs(fft(sig1);title('分解信號1頻譜')subplot(2,1,2)stem(abs(fft(sig2);title('分解信號2頻譜')% 4.MALLET重構(gòu)算法sig1=dyaddown(sig1); % 2抽取sig2=dyaddown(sig2); % 2抽取sig1=dyadup(sig1); % 2插值sig2=dyadup(sig2); % 2插值sig1=sig1(1,1:N); % 去掉最后一個零sig2=sig2(1,1:N); % 去掉最后一個零hr=h(end:-1:1); % 重構(gòu)低通g

16、r=g(end:-1:1); % 重構(gòu)高通hr=circshift(hr',1)' % 位置調(diào)整圓周右移一位gr=circshift(gr',1)' % 位置調(diào)整圓周右移一位sig1=ifft(fft(hr).*fft(sig1); % 低頻sig2=ifft(fft(gr).*fft(sig2); % 高頻sig=sig1+sig2; % 源信號% 5.比較figure(7);subplot(2,1,1)plot(real(sig1);title('重構(gòu)低頻信號');subplot(2,1,2)plot(real(sig2);title(

17、9;重構(gòu)高頻信號');figure(8);subplot(2,1,1)stem(abs(fft(sig1);title('重構(gòu)低頻信號頻譜');subplot(2,1,2)stem(abs(fft(sig2);title('重構(gòu)高頻信號頻譜');figure(9)plot(real(sig),'r','linewidth',2);hold on;plot(y);legend('重構(gòu)信號','原始信號')title('重構(gòu)信號與原始信號比較')使用小波包變換分析信號的MATLA

18、B程序%t=0.001:0.001:1;t=1:1000;s1=sin(2*pi*50*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t);for t=1:500;s2(t)=sin(2*pi*50*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t);endfor t=501:1000;s2(t)=sin(2*pi*200*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t);endsubplot(9,2,1)plot(s1)title('原始信號'

19、)ylabel('S1')subplot(9,2,2)plot(s2)title('故障信號')ylabel('S2')wpt=wpdec(s1,3,'db1','shannon');%plot(wpt);s130=wprcoef(wpt,3,0);s131=wprcoef(wpt,3,1);s132=wprcoef(wpt,3,2);s133=wprcoef(wpt,3,3);s134=wprcoef(wpt,3,4);s135=wprcoef(wpt,3,5);s136=wprcoef(wpt,3,6);s13

20、7=wprcoef(wpt,3,7);s10=norm(s130);s11=norm(s131);s12=norm(s132);s13=norm(s133);s14=norm(s134);s15=norm(s135);s16=norm(s136);s17=norm(s137);st10=std(s130);st11=std(s131);st12=std(s132);st13=std(s133);st14=std(s134);st15=std(s135);st16=std(s136);st17=std(s137);disp('正常信號的特征向量');snorm1=s10,s11,

21、s12,s13,s14,s15,s16,s17std1=st10,st11,st12,st13,st14,st15,st16,st17subplot(9,2,3);plot(s130);ylabel('S130');subplot(9,2,5);plot(s131);ylabel('S131');subplot(9,2,7);plot(s132);ylabel('S132');subplot(9,2,9);plot(s133);ylabel('S133');subplot(9,2,11);plot(s134);ylabel(

22、9;S134');subplot(9,2,13);plot(s135);ylabel('S135');subplot(9,2,15);plot(s136);ylabel('S136');subplot(9,2,17);plot(s137);ylabel('S137');wpt=wpdec(s2,3,'db1','shannon');%plot(wpt);s230=wprcoef(wpt,3,0);s231=wprcoef(wpt,3,1);s232=wprcoef(wpt,3,2);s233=wprcoef

23、(wpt,3,3);s234=wprcoef(wpt,3,4);s235=wprcoef(wpt,3,5);s236=wprcoef(wpt,3,6);s237=wprcoef(wpt,3,7);s20=norm(s230);s21=norm(s231);s22=norm(s232);s23=norm(s233);s24=norm(s234);s25=norm(s235);s26=norm(s236);s27=norm(s237);st20=std(s230);st21=std(s231);st22=std(s232);st23=std(s233);st24=std(s234);st25=st

24、d(s235);st26=std(s236);st27=std(s237);disp('故障信號的特征向量');snorm2=s20,s21,s22,s23,s24,s25,s26,s27std2=st20,st21,st22,st23,st24,st25,st26,st27subplot(9,2,4);plot(s230);ylabel('S230');subplot(9,2,6);plot(s231);ylabel('S231');subplot(9,2,8);plot(s232);ylabel('S232');subplot

25、(9,2,10);plot(s233);ylabel('S233');subplot(9,2,12);plot(s234);ylabel('S234');subplot(9,2,14);plot(s235);ylabel('S235');subplot(9,2,16);plot(s236);ylabel('S236');subplot(9,2,18);plot(s237);ylabel('S237');%fftfigurey1=fft(s1,1024);py1=y1.*conj(y1)/1024;y2=fft(s2

26、,1024);py2=y2.*conj(y2)/1024;y130=fft(s130,1024);py130=y130.*conj(y130)/1024;y131=fft(s131,1024);py131=y131.*conj(y131)/1024;y132=fft(s132,1024);py132=y132.*conj(y132)/1024;y133=fft(s133,1024);py133=y133.*conj(y133)/1024;y134=fft(s134,1024);py134=y134.*conj(y134)/1024;y135=fft(s135,1024);py135=y135.

27、*conj(y135)/1024;y136=fft(s136,1024);py136=y136.*conj(y136)/1024;y137=fft(s137,1024);py137=y137.*conj(y137)/1024;y230=fft(s230,1024);py230=y230.*conj(y230)/1024;y231=fft(s231,1024);py231=y231.*conj(y231)/1024;y232=fft(s232,1024);py232=y232.*conj(y232)/1024;y233=fft(s233,1024);py233=y233.*conj(y233)/

28、1024;y234=fft(s234,1024);py234=y234.*conj(y234)/1024;y235=fft(s235,1024);py235=y235.*conj(y235)/1024;y236=fft(s236,1024);py236=y236.*conj(y236)/1024;y237=fft(s237,1024);py237=y237.*conj(y237)/1024;f=1000*(0:511)/1024;subplot(1,2,1);plot(f,py1(1:512);ylabel('P1');title('原始信號的功率譜')subp

29、lot(1,2,2);plot(f,py2(1:512);ylabel('P2');title('故障信號的功率譜')figuresubplot(4,2,1);plot(f,py130(1:512);ylabel('P130');title('S130的功率譜')subplot(4,2,2);plot(f,py131(1:512);ylabel('P131');title('S131的功率譜')subplot(4,2,3);plot(f,py132(1:512);ylabel('P132&#

30、39;);subplot(4,2,4);plot(f,py133(1:512);ylabel('P133');subplot(4,2,5);plot(f,py134(1:512);ylabel('P134');subplot(4,2,6);plot(f,py135(1:512);ylabel('P135');subplot(4,2,7);plot(f,py136(1:512);ylabel('P136');subplot(4,2,8);plot(f,py137(1:512);ylabel('P137');figur

31、esubplot(4,2,1);plot(f,py230(1:512);ylabel('P230');title('S230的功率譜')subplot(4,2,2);plot(f,py231(1:512);ylabel('P231');title('S231的功率譜')subplot(4,2,3);plot(f,py232(1:512);ylabel('P232');subplot(4,2,4);plot(f,py233(1:512);ylabel('P233');subplot(4,2,5);pl

32、ot(f,py234(1:512);ylabel('P234');subplot(4,2,6);plot(f,py235(1:512);ylabel('P235');subplot(4,2,7);plot(f,py236(1:512);ylabel('P236');subplot(4,2,8);plot(f,py237(1:512);ylabel('P237');figure%plottree(wpt)基于小波變換的圖象去噪 Normalshrink算法function T_img,Sub_T=threshold_2_N(img,

33、levels)% reference :image denoising using wavelet thresholdingxx,yy=size(img);HH=img(xx/2+1):xx,(yy/2+1):yy);delt_2=(std(HH(:)2;%(median(abs(HH(:)/0.6745)2;%T_img=img;for i=1:levels temp_x=xx/2i; temp_y=yy/2i;% belt=1.0*(log(temp_x/(2*levels)0.5; belt=1.0*(log(temp_x/(2*levels)0.5; %2.5 0.8 %HL HL=i

34、mg(1:temp_x,(temp_y+1):2*temp_y); delt_y=std(HL(:); T_1=belt*delt_2/delt_y; % T_HL=sign(HL).*max(0,abs(HL)-T_1); % T_img(1:temp_x,(temp_y+1):2*temp_y)=T_HL; Sub_T(3*(i-1)+1)=T_1; %LH LH=img(temp_x+1):2*temp_x,1:temp_y); delt_y=std(LH(:); T_2=belt*delt_2/delt_y; % T_LH=sign(LH).*max(0,abs(LH)-T_2); %

35、 T_img(temp_x+1):2*temp_x,1:temp_y)=T_LH; Sub_T(3*(i-1)+2)=T_2; %HH HH=img(temp_x+1):2*temp_x,(temp_y+1):2*temp_y); delt_y=std(HH(:); T_3=belt*delt_2/delt_y; % T_HH=sign(HH).*max(0,abs(HH)-T_3); % T_img(temp_x+1):2*temp_x,(temp_y+1):2*temp_y)=T_HH; Sub_T(3*(i-1)+3)=T_3;end用小波神經(jīng)網(wǎng)絡(luò)來對時間序列進(jìn)行預(yù)測%File name

36、 : nprogram.m %Description : This file reads the data from its source into their respective matrices prior to % performing wavelet decomposition. % %Clear command screen and variables clc; clear; % % user desired resolution level (Tested: resolution = 2 is best) level = menu('Enter desired resol

37、ution level: ', '1',. '2 (Select this for testing)', '3', '4'); switch level case 1, resolution = 1; case 2, resolution = 2; case 3, resolution = 3; case 4, resolution = 4; end msg = 'Resolution level to be used is ', num2str(resolution); disp(msg); % % us

38、er desired amount of data to use data = menu('Choose amount of data to use: ', '1 day', '2 days', '3 days', '4 days',. '5 days', '6 days', '1 week (Select this for testing)'); switch data case 1, dataPoints = 48; %1 day = 48 points case

39、 2, dataPoints = 96; %2 days = 96 points case 3, dataPoints = 144; %3 days = 144 points case 4, dataPoints = 192; %4 days = 192 points case 5, dataPoints = 240; %5 days = 240 points case 6, dataPoints = 288; %6 days = 288 points case 7, dataPoints = 336; %1 weeks = 336 points end msg = 'No. of d

40、ata points to be used is ', num2str(dataPoints); disp(msg); % %Menu for data set selection select = menu('Use QLD data of: ', 'Jan02',. 'Feb02', 'Mar02 (Select this for testing)', 'Apr02', 'May02'); switch select case 1, demandFile = 'DATA20060

41、1_QLD1' case 2, demandFile = 'DATA200602_QLD1' case 3, demandFile = 'DATA200603_QLD1' case 4, demandFile = 'DATA200604_QLD1' case 5, demandFile = 'DATA200605_QLD1' end % %Reading the historical DEMAND data into tDemandArray selectedDemandFile=demandFile,'.csv&

42、#39;regionArray, sDateArray, tDemandArray, rrpArray, pTypeArray . = textread(selectedDemandFile, '%s %q %f %f %s', 'headerlines', 1, 'delimiter', ','); %Display no. of points in the selected time series demand data demandDataPoints, y = size(tDemandArray); msg = '

43、The no. of points in the selected Demand data is ', num2str(demandDataPoints); disp(msg); %Decompose historical demand data signal dD, l = swtmat(tDemandArray, resolution, 'db2'); approx = dD(resolution, :); %Plot the original demand data signal figure (1); subplot(resolution + 2, 1, 1);

44、 plot(tDemandArray(1: dataPoints) legend('Demand Original'); title('QLD Demand Data Signal'); %Plot the approximation demand data signal for i = 1 : resolution subplot(resolution + 2, 1, i + 1); plot(approx(1: dataPoints) legend('Demand Approximation'); end %After displaying

45、approximation signal, display detail x for i = 1: resolution if( i > 1 ) detail(i, :) = dD(i-1, :)- dD(i, :); else detail(i, :) = tDemandArray' - dD(1, :); end if i = 1 subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints) legendName = 'Demand Detail ', num2str

46、(i); legend(legendName); else subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints) legendName = 'Demand Detail ', num2str(i); legend(legendName); end i = i + 1; end %Normalising approximation demand data maxDemand = max(approx'); %Find largest component normDeman

47、d = approx ./ maxDemand; %Right divison maxDemandDetail = ; normDemandDetail = , ; detail = detail + 4000; for i = 1: resolution maxDemandDetail(i) = max(detail(i, :); normDemandDetail(i, :) = detail(i, :) ./maxDemandDetail(i); end % %Reading the historical historical PRICE data into rrpArray select

48、edPriceFile = demandFile, '.csv' regionArray, sDateArray, tDemandArray, rrpArray, pTypeArray . = textread(selectedDemandFile, '%s %q %f %f %s', 'headerlines', 1, 'delimiter', ','); %Display no. of points in Price data noOfDataPoints, y = size(rrpArray); msg =

49、'The no. of points in Price data data is ', num2str(noOfDataPoints); disp(msg); %Decompose historical Price data signal dP, l = swtmat(rrpArray, resolution, 'db2'); approxP = dP(resolution, :); %Plot the original Price data signal figure (2); subplot(resolution + 3, 1, 1); plot(rrpAr

50、ray(2: dataPoints) legend('Price Original'); title('Price Data Signal'); %Plot the approximation Price data signal for i = 1 : resolution subplot(resolution + 3, 1, i + 1); plot(approxP(2: dataPoints) legend('Price Approximation'); end %After displaying approximation signal,

51、display detail x for i = 1: resolution if( i > 1 ) detailP(i, :) = dP(i-1, :)- dP(i, :); else detailP(i, :) = rrpArray' - dP(1, :); end if i = 1 B,A=butter(1,0.65,'low'); result =filter(B,A, detailP(i, 1: dataPoints); subplot(resolution + 3, 1, resolution - i + 4);plot(result(i, 2: da

52、taPoints) legendName = 'low pass filter', num2str(i); legend(legendName); subplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i, 2: dataPoints) legendName = 'Price Detail ', num2str(i); legend(legendName); else subplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i,

53、2: dataPoints) legendName = 'Price Detail ', num2str(i); legend(legendName); end i = i + 1; end %Normalising approximation Price data maxPrice = max(approxP'); %Find largest component normPrice = approxP ./ maxPrice; %Right divison maxPriceDetail = ; normPriceDetail = , ; detailP = detailP + 40; for i = 1: resolution maxPriceDetail(i) = max(detailP(i, :); normPriceDetail(i, :) = detailP(i, :) ./maxPriceDetail(i); en