小波變換及其應(yīng)用2014

1、小波變換及其應(yīng)用 一 什么是小波傅里葉變換：能表示信號(hào) f 的各頻率成分的大小，且但不能體現(xiàn)這個(gè)頻率成分是哪個(gè)時(shí)刻或哪個(gè)時(shí)間段的。即不能提供時(shí)-頻定位信息。因此，傅里葉變換是分析平穩(wěn)信號(hào)的有力工具。讓我們從信號(hào) 的分析談起！ 假設(shè)信號(hào)的持續(xù)時(shí)間為 ，含最高頻率 。采樣間隔： 傅里葉變換離散型DFT（或FFT）： 一 什么是小波 對(duì)于給定的非平穩(wěn)信號(hào) ,人們感興趣的是信號(hào)在特定的時(shí)間的頻率成分。就像在音樂(lè)演奏中，演奏者需要知道在什么時(shí)候演奏什么音調(diào)一樣。給出了信號(hào) f 的時(shí)-頻局部定位信息，且 為了進(jìn)行時(shí)-頻定位分析，引入窗函數(shù) ，形成加窗傅里葉變換： 一 什么是小波窗的中心：窗半徑：例：窗函數(shù)

2、 的作用 一 什么是小波平面上 的窗： 可見(jiàn)，當(dāng)窗函數(shù)確定后，窗不隨時(shí)間和頻率的變化而變化。 一 什么是小波表達(dá)了信號(hào)在窗 內(nèi)的信息. 的減小，時(shí)間分辨率變低，頻率分辨率變高； 的增加，時(shí)間分辨率變高，頻率分辨率變低。 一 什么是小波加窗傅里葉變換的離散型： 一 什么是小波演示：function yanshi(A,B,p)fs=3000;Tp=50;T=1/fs;N=Tp/T;F=fs/N;t=-(N-1)*T:T:(N-1)*T;f=0:F:(N-1)*F;y=exp(-t).*(sin(2*pi*100*t)+sin(2*pi*1000*t);%信號(hào)z=A*exp(-B*(t-p).2);

3、%窗函數(shù)x=z.*y;%時(shí)刻p時(shí)窗下信號(hào)s=x(N:2*N-1);c=fft(s);subplot(3,1,1);plot(t,x);subplot(3,1,2);c1=c(4975:5025)/300;plot(f(4975:5025),abs(c1);subplot(3,1,3);c2=c(49950:50050)/300;plot(f(49950:50050),abs(c2); 一 什么是小波yanshi(1,1/10,3) 一 什么是小波yanshi(1,1/2,3) yanshi(1,1,3) 為了更有效確定信號(hào)的時(shí)-頻信息，引入具有下列特性的窗函數(shù)：小波的概念是J.Morlet(莫

4、萊)在進(jìn)行地震數(shù)據(jù)分析中提出的。它是窗函數(shù)因?yàn)?稱 為（母）小波函數(shù)所以， 一 什么是小波 對(duì)母小波函數(shù)進(jìn)行伸縮（伸縮因子為a0 ）和平移變換（平移因子為 ），得下列函數(shù)族稱為分析小波(小波基函數(shù))。 一 什么是小波 母小波函數(shù) 是單窗函數(shù)！設(shè)其時(shí)窗中心 為 ，頻窗中心為 ，時(shí)窗半徑 ，頻窗半徑為 。則窗為:低頻高頻b代表時(shí)間！代表頻率，（1），窗面積恒定?。?），窄窗看高頻；，寬窗看低頻。（3） 分析小波 具有變焦性能： 一 什么是小波稱為函數(shù) 的連續(xù)小波變換。 這里，當(dāng)（母）小波函數(shù) 滿足條件時(shí)， 有： 二 連續(xù)小波變換段 上，包含中心頻率 ，帶寬 的窗內(nèi)的頻率分量的大小。 的連續(xù)小波變換

5、系數(shù) 表示：位置b處，時(shí)間MATLAB實(shí)現(xiàn)： COEFS = CWT(S,SCALES,wname,PLOTMODE，XLIM) PLOTMODE = plot (By scale) or PLOTMODE = lvl (By scale) or PLOTMODE = glb (All scales) or PLOTMODE = abslvl or lvlabs (Absolute value and By scale) or PLOTMODE = absglb or glbabs (Absolute value and All scales) XLIM:圖的范圍 二 連續(xù)小波變換例1：連續(xù)小

6、波變換fc=1000;Tp=1;T=1/2048;N=1/T;t=0:T:(N-1)*T;s=sin(2*pi*500*t)+sin(2*pi*1000*t);ss=s;ss(165)=3;ss(207)=3;subplot(2,1,1),plot(s);subplot(2,1,2),plot(ss) 二 連續(xù)小波變換c=cwt(ss,0.01:5,db4,plot); 二 連續(xù)小波變換 c=cwt(ss,0.01:5,coif5,plot); 二 連續(xù)小波變換 三 正交小波將a、b離散化，取 ，則有或連續(xù)小波基函數(shù) 具有很大的相關(guān)性，因此 之間是有冗余的。大部分工作是希望分解系數(shù) 之間沒(méi)有冗

7、余。我們知道，函數(shù)按照線性無(wú)關(guān)的基展開(kāi)的系數(shù)沒(méi)有冗余，特別是標(biāo)準(zhǔn)正交基更好！如果 構(gòu)成 的標(biāo)準(zhǔn)正交基，則 其中，稱為正交小波變換 表示：在位置 處，時(shí)間段 上，包含中心頻率 ，帶寬 的窗內(nèi)的頻率分量的大小. 三 正交小波多尺度分析（MRA）：1. 選取一個(gè)具有緊支集的函數(shù) （稱為尺度函數(shù)），且 關(guān)于k是兩兩標(biāo)準(zhǔn)正交的。定義線性空間2. 令 , 關(guān)于k是兩兩標(biāo)準(zhǔn)正交的。建立線性空間 三 正交小波（尺度空間 ）具有下列性質(zhì)：（2）伸縮規(guī)則性：（1）一致單調(diào)性：（3）漸進(jìn)完全性： （4）平移不變性：易證： 是 的標(biāo)準(zhǔn)正交基。但是，不能構(gòu)成 的標(biāo)準(zhǔn)正交基。 三 正交小波令容易證明：3.尋找 的標(biāo)準(zhǔn)正交

8、基其中，稱為低通濾波器系數(shù)。由 ，可得兩尺度方程： 三 正交小波聯(lián)立式（1）和式（2）,并結(jié)合 正交性，推出：定為正交小波函數(shù)! ，有其中， 稱為高通濾波器系數(shù)。 三 正交小波構(gòu)成 的標(biāo)準(zhǔn)正交基；構(gòu)成 的標(biāo)準(zhǔn)正交基。 三 正交小波綜上所述，尋找正交小波的方法為：（1）由正交尺度函數(shù) ，計(jì)算（2）計(jì)算（3）確定正交小波函數(shù) 三 正交小波例：構(gòu)造shannon小波取正交尺度函數(shù) ， 三 正交小波 三 正交小波常用小波函數(shù)29四 正交小波分析一層分解： 設(shè) 的實(shí)測(cè)信號(hào)為 ，則 其中m層分解：二層分解：四 正交小波分析 mallat算法分解：重構(gòu):四 正交小波分析一維離散小波變換的MATLAB實(shí)現(xiàn)重構(gòu)

9、命令:s=idwt(ca,cd,wname) ;% 單層重構(gòu)s=waverec(c,l,wname) ;%多層重構(gòu)s=wrcoef(type,c,l,wname,N) ;%重構(gòu)指定層（N）的系數(shù)，type=a或d.s=appcoef(c,l,wname,N) ;%提取N層的近似系數(shù)s=detcoef(c,l,N) ;%提取N層的細(xì)節(jié)系數(shù)分解命令：ca,cd=dwt(s,wname,mode,MODE) ,單層分解c,l=wavedec(s,n,wname) , n層分解ca,cd=dwt(s,Lo_D,Hi_D) ,單層分解function yanshi2k=1:2048;s=sin(2*pi

10、*500*(k-1)/2048)+sin(2*pi*1000*(k-1)/2048);ss=s;ss(165)=3;ss(207)=3;c,l=wavedec(ss,4,sym1);cd1,cd2,cd3,cd4=detcoef(c,l,1,2,3,4);ca4=appcoef(c,l,sym1,4);subplot(321),plot(ss);title(ss);subplot(322),plot(ca4);title(ca4);subplot(323),plot(cd4);title(cd4);subplot(324),plot(cd3);title(cd3);subplot(325),p

11、lot(cd2);title(cd2);subplot(326),plot(cd1);title(cd1);3435function yanshi3k=1:2048;s=sin(2*pi*500*(k-1)/2048)+sin(2*pi*1000*(k-1)/2048);ss=s;ss(165)=3;ss(207)=3;c,l=wavedec(ss,4,sym1);cd1,cd2,cd3,cd4=detcoef(c,l,1,2,3,4);ca4=appcoef(c,l,sym1,4);ca3=idwt(ca4,cd4,sym1);ca2=idwt(ca3,cd3,sym1);ca1=idwt(

12、ca2,cd2,sym1);x=idwt(ca1,cd1,sym1);subplot(321),plot(ca4),title(ca4);subplot(322),plot(ca3),title(ca3);subplot(323),plot(ca2),title(ca2);subplot(324),plot(ca1),title(ca1);subplot(325),plot(x),title(x);subplot(326),plot(ss),title(ss);3637 給定二維信號(hào) , 其實(shí)測(cè)信號(hào)為:四 正交小波分析一層分解:四 正交小波分析二層分解：一直分解下去 二維離散小波變換matla

13、b實(shí)現(xiàn)分解命令: CA,CH,CV,CD = dwt2(X,wname)；%單層分解C,S = wavedec2(X,N,wname)；%N層分解S(1,:);%第 N層近似系數(shù)的長(zhǎng)度S(i,:),i=2:N+1;%第 (N-i+2)層細(xì)節(jié)系數(shù)的長(zhǎng)度S(N+2,:),%原始信號(hào)的長(zhǎng)度x=idwt2(ca,ch,cv,cd,wname); %本層系數(shù)重建上一層近似系數(shù)xx=waverec2(c,s,wname) ; %多層重構(gòu)原始信號(hào)x=wrcoef2(type,c,s,wname,N) ; %重構(gòu)N層的類type=a,h,v,d的系數(shù)重構(gòu)命令:41D=detcoef2(o,c,s,wname,

14、N) ; %提取第N層o型的系數(shù)，其中o=h,v,d,compactH,V,D=detcoef2(all,c,s,wname,N) ; %提取第N層所有的細(xì)節(jié)系數(shù)A=appcoef2(c,s,wname,N) ; %提取第N層近似系數(shù)更多重構(gòu)命令四 正交小波分析42 load clown; whos Name Size Bytes Class Attributes X 200 x320 512000 double caption 2x1 4 char map 81x3 1944 double image(X);colormap(map),colorbar;索引表不連續(xù)，不能對(duì)X直接進(jìn)行小波變換

15、！四 正交小波分析% 依索引表提取三個(gè)通道的分量值 R=map(X,1);R=reshape(R,size(X);G=map(X,2);G=reshape(G,size(X);B=map(X,3);B=reshape(B,size(X);%轉(zhuǎn)為灰度矩陣X1=0.2990*R+0.5870*G+0.1140*B;%給定灰度階數(shù)n=256;%從灰度矩陣轉(zhuǎn)化為索引號(hào)矩陣X2=round(X1*(n-1)+1; map2=gray(n); );%建立灰度索引表image(X2);colormap(map2),colorbar;索引表連續(xù)了，能進(jìn)行小波變換，但沒(méi)有顏色信息了！四 正交小波分析load n

16、oiswom;whosName Size Bytes Class Attributes X 96x96 73728 double map 255x3 6120 double image(X);colormap(map),colorbar;索引表連續(xù)，可直接進(jìn)行小波變換。ca1,ch1,cv1,cd1=dwt2(X,db4);A1=upcoef2(a,ca1,db4,1);H1=upcoef2(h,ch1,db4,1);V1=upcoef2(v,cv1,db4,1);D1=upcoef2(d,cd1,db4,1);colormap(map);subplot(2,2,1);image(wcodem

17、at(A1,96); subplot(2,2,2);image(wcodemat(H1,96);subplot(2,2,3);image(wcodemat(V1,96); subplot(2,2,4);image(wcodemat(D1,96); 四 正交小波分析colormap(map);subplot(2,2,1);image(wcodemat(A1,96); subplot(2,2,2);image(wcodemat(H1,96);subplot(2,2,3);image(wcodemat(V1,96); subplot(2,2,4);image(wcodemat(D1,96); 四 正

18、交小波分析四 正交小波分析實(shí)際應(yīng)用時(shí)，選擇小波的原則：1. 自相似 ；2. 支集長(zhǎng)度3. 對(duì)稱性,消失矩階數(shù),正則性對(duì)稱性用于圖像處理可避免相移;消失矩階數(shù)高的小波變換能量更集中;正則性好可使重構(gòu)的信號(hào)更平滑;支集長(zhǎng)度越長(zhǎng)消失矩和正則性越高！五 小波應(yīng)用舉例降噪準(zhǔn)則： 1. 光滑性，即降噪后的信號(hào)至少與原信號(hào)具有相同的光滑性； 2. 相似性，即降噪后的信號(hào) 與原信號(hào)的方差最小，即應(yīng)用1： 小波變換用于信號(hào)降噪(壓縮)sigma=wnoisest(c,l,s); %估算第s層的細(xì)節(jié)系數(shù)的標(biāo)準(zhǔn)差，作為噪聲強(qiáng)度；噪聲模型為：這里， 為信號(hào)， 為被噪聲 污染后的信號(hào) ， 為噪聲強(qiáng)度。 51小波降噪過(guò)程

19、：1.信號(hào)的小波分解，一般需要多層分解；2.對(duì)細(xì)節(jié)系數(shù)進(jìn)行閥值處理（降噪）；3.重構(gòu)。關(guān)鍵是閥值！缺省策略 （全局閾值）：thr,sorh,keepapp=ddencmp(o,wv,x); %x為原始信號(hào)，o=den降噪，o=cmp壓縮; thr閾值，sorh=s軟閾值， sorh=h硬閾值, keepapp保留的近似系數(shù)的層數(shù)。Birge-Massart策略（分層閾值）：thr,nkeep=wdcbm(c,l, alpha)； %Alpha=1.5壓縮，=2，3降噪。thr各層 閾值，nkeep保留的近似系數(shù)的層數(shù)。 penalty策略(用于小波包)：thr=wbmpen(c,l,sigma

20、,alpha)；thr=thselet(y,tptr)；% 用指定的閾值選取方法tptr（tptr=rigsure嚴(yán)格無(wú)偏估計(jì)閾值選擇、sqtwolog對(duì)數(shù)長(zhǎng)度、heusure啟發(fā)式無(wú)偏估計(jì)、minimaxi最大最小方差），得到對(duì)信號(hào)y降噪的閾值thr。應(yīng)用1： 小波變換用于信號(hào)降噪(壓縮)確定閾值的命令53xc,cxc,lxc,perf0,perfl2=wdencmp(gbl,x,wname,N,thr,sorh,k)%全局(gbl)閾值的降噪（壓縮）命令,保留k層以上的所有系數(shù), PERF0、PERFL2降噪或壓縮所保留的能量成分（%）。xc,cxc,lxc,perf0,perfl2=wd

21、encmp(lvd,x,wname,N,thr,sorh)%分層(lvd）閾值的降噪（壓縮）命令。THR各層閾值。xc,cxc,lxc,perf0,perfl2=wdencmp(lvd,c,l,wname,N,thr,sorh)%分層(lvd）閾值的降噪（壓縮）命令。THR各層閾值。Matlab降噪（壓縮）函數(shù)xd,cxd,lxd = wden(x,tptr,sorh,scal,N,wname)；%根據(jù)傳入的參數(shù)對(duì)信號(hào)x降噪，結(jié)果返回xd，cxd,lxd是xd的小波分解系數(shù)。Tptr如上，sorh=s軟閾值、h硬閾值，scal作用閾值時(shí)的乘子，scal=one不使用乘子、sln用對(duì)第一層細(xì)節(jié)系

22、數(shù)的估計(jì)作為乘子、mln用層數(shù)相關(guān)的乘子，N分解層數(shù)。54一維信號(hào)降噪function yanshi10load noisbump; x=noisbump;c,l=wavedec(x,5,sym6);thr,sorh,keepapp=ddencmp(den,wv,x); %全局降噪的閥值thr。thr1,nkeep=wdcbm(c,l,2); %使用birge-massart策略確定各層降噪的閥值。%全局軟閾值降噪信號(hào)xc,保留一層以上的所有分解系數(shù)cxc,lxcxc,cxc,lxc,perf0,perfl2 =wdencmp(gbl,c,l,sym6,5,thr,s,1); %分層軟閥值降噪

23、xc;xc1,cxc1,lxc1,perf01,perfl21 =wdencmp(lvd,c,l,sym6,5,thr1,s);subplot(3,1,1);plot(x);title(原始信號(hào),fontsize,8);subplot(3,1,2);plot(xc);title(全局軟閾值降噪信號(hào)，能量成分perf01,perfl21,fontsize,8);subplot(3,1,3);plot(xc1);title(分層軟閥值降噪信號(hào)，能量成分perf01,perfl21,fontsize,8);二維信號(hào)降噪load womanrandn(seed,2055615866);x=X+15*r

24、andn(size(X); %加噪 thr,sorh,keepapp=ddencmp(den,wv,x);%全局閥值降噪xc,cxc,lxc,perf0,perfl2=wdencmp(gbl,x,sym4,2,thr,sorh,keepapp); %顯示colormap(pink(255),sm=size(map,1);subplot(221),image(wcodemat(X,sm);title(X);subplot(222),image(wcodemat(x,sm);title(x);subplot(223),image(wcodemat(xc,sm);title(xc);c,s=wave

25、dec2(x,3,sym4);%三層二維小波分解 th1,nkeep1=wdcbm2(c,s,1.5,2.7*prod(s(1,:);xc1,cxc1,lxc1,perf01,perfl21=wdencmp(lvd,c,s,sym4,3,th1,s); subplot(224),image(wcodemat(xc1,sm);title(xc1);應(yīng)用2：小波變換用于信號(hào)奇異點(diǎn)檢測(cè) 信號(hào)的奇異點(diǎn)（突變點(diǎn)）常包含重要的故障信息（對(duì)于圖像它是邊緣）。提取此突變點(diǎn)的位置及奇異性（光滑度）是信號(hào)分析與處理的重要工作之一函數(shù)奇異性的描述，李氏（Lipschitz）指數(shù)：連續(xù)不可導(dǎo)函數(shù)：稱函數(shù)f(t)在點(diǎn)

26、的李氏指數(shù)為（ ），如果滿足：不連續(xù)函數(shù)：稱函數(shù)f(t)在點(diǎn) 的李氏指數(shù)（ ），如果其原函數(shù)F(t)的李氏指數(shù)為+1。n階連續(xù)導(dǎo)數(shù)，但n+1階不可導(dǎo)函數(shù)：稱函數(shù)f(t)在點(diǎn) 的李氏指數(shù)為（ ），如果滿足：例如：李氏指數(shù)=1李氏指數(shù)=0李氏指數(shù)=-1函數(shù)奇異性在小波變換下的特征： 設(shè) 的李氏指數(shù)為 (0 =1），則其小波變換系數(shù)滿足： 是小波變換模的上界（模極大值），若在某時(shí)刻小波變換模達(dá)到這個(gè)上界，則說(shuō)明f(t)在此處的李氏指數(shù)為。 另外， 時(shí)，小波變換模隨j增而增； 時(shí)，小波變換模隨j減而減； 時(shí)，小波變換模不隨j變！ 綜上，用小波變換模極大值，可以確定信號(hào)奇異點(diǎn)及奇異性。64連續(xù)間斷點(diǎn)的檢

27、測(cè)function yanshi13load freqbrks=freqbrk;c,l=wavedec(s,5,db5);for i=1:5 ca(i,:)=wrcoef(a,c,l,db5,i); cd(i,:)=wrcoef(d,c,l,db5,i);endsubplot(611),plot(s);title(s);subplot(612),plot(ca(1,:);title(ca1);subplot(613),plot(ca(2,:);title(ca2);subplot(614),plot(ca(3,:);title(ca3);subplot(615),plot(ca(4,:);ti

28、tle(ca4);subplot(616),plot(ca(5,:);title(ca5);figure;subplot(612),plot(cd(1,:);title(cd1);subplot(613),plot(cd(2,:);title(cd2);subplot(614),plot(cd(3,:);title(cd3);subplot(615),plot(cd(4,:);title(cd4);subplot(616),plot(cd(5,:);title(cd5);end導(dǎo)數(shù)間斷點(diǎn)的檢測(cè)load scddvbrk;s=scddvbrk;c,l=wavedec(s,2,db4);for i

29、=1:2 ca(i,:)=wrcoef(a,c,l,db4,i); cd(i,:)=wrcoef(d,c,l,db4,i);endt=350:650;subplot(321),plot(t,s(t);title(s);subplot(323),plot(t,ca(1,t);title(ca1);subplot(324),plot(t,cd(1,t);title(cd1)；subplot(325),plot(t,ca(2,t);title(ca2);subplot(326),plot(t,cd(2,t);title(cd2)；應(yīng)用3 識(shí)別信號(hào)發(fā)展趨勢(shì)function yanshi15load c

30、noislop;s=cnoislop;c,l=wavedec(s,6,db3);for i=1:6 ca(i,:)=wrcoef(a,c,l,db3,i); cd(i,:)=wrcoef(d,c,l,db3,i);endsubplot(711),plot(s),title(s);subplot(712),plot(ca(1,:);title(ca1);subplot(713),plot(ca(2,:);title(ca2);subplot(714),plot(ca(3,:);title(ca3);subplot(715),plot(ca(4,:);title(ca4);subplot(716)

31、,plot(ca(5,:);title(ca5);subplot(717),plot(ca(6,:);title(ca6);figure;subplot(712),plot(cd(1,:);title(cd1);subplot(713),plot(cd(2,:);title(cd2);subplot(714),plot(cd(3,:);title(cd3);subplot(715),plot(cd(4,:);title(cd4);subplot(716),plot(cd(5,:);title(cd5);subplot(717),plot(cd(6,:);title(cd6);end在一個(gè)電廠采

32、集的三天負(fù)荷信號(hào)s=leleccum 。具體是：1. 采樣時(shí)間間隔1分鐘；2. 數(shù)據(jù)的采集是通過(guò)幾百個(gè)高度密集的傳感器完成的；3. 負(fù)荷數(shù)據(jù)大體一半是來(lái)自于工業(yè)，一半來(lái)自于大量個(gè)人用戶。來(lái)自于工業(yè)的負(fù)荷成分規(guī)律性更強(qiáng)，體現(xiàn)為緩慢變化的低頻成分；與此相反，個(gè)人用戶的負(fù)荷成分變化更為劇烈，其信號(hào)成分也以高頻為主； 噪聲可能來(lái)源傳感器或系統(tǒng)狀態(tài)！ 應(yīng)用4 電網(wǎng)監(jiān)測(cè)準(zhǔn)備研究以下問(wèn)題：每日中部時(shí)間段（36003700或12：0013：30）信號(hào)分析；深夜時(shí)段（15301730或01：3004：30）的信號(hào)分析；檢測(cè)傳感器故障；抑制噪聲；識(shí)別信號(hào)模式；定位并抑制偏離點(diǎn)；研究丟失信息。function ya

33、nshi23load leleccum; s=leleccum;c,l=wavedec(s,5,db3); for i=1:5 ca(i,:)=wrcoef(a,c,l,db3,i);cd(i,:)=wrcoef(d,c,l,db3,i);endsubplot(311),plot(s),title(s);subplot(312),plot(ca(1,:);title(ca1);subplot(313),plot(cd(1,:);title(cd1);figure;subplot(311),plot(ca(1,:);title(ca1);subplot(312),plot(ca(2,:);tit

34、le(ca2);subplot(313),plot(cd(2,:);title(cd2);figure;subplot(311),plot(ca(2,:);title(ca2);subplot(312),plot(ca(3,:);title(ca3);subplot(313),plot(cd(3,:);title(cd3);figure;subplot(311),plot(ca(3,:);title(ca3);subplot(312),plot(ca(4,:);title(ca4);subplot(313),plot(cd(4,:);title(cd4);figure;subplot(311)

35、,plot(ca(4,:);title(ca4);subplot(223),plot(ca(5,:);title(ca5);subplot(224),plot(cd(5,:);title(cd5);1. 每日中部時(shí)間段負(fù)荷數(shù)據(jù)分析function yanshi19load leleccum; s=leleccum;t=3600:3700;c,l=wavedec(s(t),5,db3); for i=1:5 ca(i,:)=wrcoef(a,c,l,db3,i); cd(i,:)=wrcoef(d,c,l,db3,i);endsubplot(221),plot(t,s(t),title(s(3600:3700);subplot(223),plot(t,ca(1,:);title(ca1);subplot(224),plot(t,cd(1,:);title(cd1);figure;subplot(221),plot(t,ca(2,:);title(c