小波分析的matlab程式

飄香的城堡發表於2011-06-23

完整版

[分享]個人收集的一些關於小波分析的matlab程式
http://www.chinavib.com/thread-9141-1-1.html

中國振動聯盟's Archiver

simon21 發表於 2006-3-31 08:34

[分享]個人收集的一些關於小波分析的matlab程式

都是從網上收集來的,由於時間比較久,處處都忘記了,如果是誰的原創請和我聯絡,我在帖子中標出來的
內容比較多,將會逐步貼出來

提升法97經典程式 (二樓)
2代小波示意程式 (三樓)
二代小波漫談 (四樓)
小波濾波器構造和消噪程式(2個) (五樓)
小波譜分析mallat演算法經典程式 (六樓)
2維小波變換經典程式 (七樓)
基於LeventCodes平臺的小波去噪程式包 (十一樓)
連續小波和離散小波分析的應用例項(十二樓)
小波插值與小波構造(3個程式)(十三樓)
採用多孔trous演算法(undecimated wavelet transform)實現小波變換(十四樓)
Daubechies小波基的構造(十五樓)
消失矩作用的程式(二十三樓)
平移變換平移法(cycle_spinning)消除gibbs效應 (二十四樓)
使用小波包變換分析訊號的MATLAB程式(五十四樓)
基於小波消噪的雷達回波檢測,可以檢測雷達回波的有無及其準確的位置(五十五樓)
二維小波變換(正和逆變換)(五十六樓)
第二代小波變換原始碼(五十七樓)
利用小波變換實現對電能質量檢測的演算法實現(五十八樓)
基於小波變換的圖象去噪 Normalshrink演算法(五十九樓)
基於小波變換模極大的多尺度影象邊緣檢測(六十樓)
利用小波變根據二進位制數(水印)來改變圖片,提取其中一個子帶的直方圖(六十一樓)
用小波函式構建神經網路的源程式(六十二樓)
利用小波和霍夫曼對單聲道檔案進行壓縮編碼 並解碼輸出(六十三樓)
用小波神經網路來對時間序列進行預測(六十四樓)
基於小波特徵提取方法的圖象匹配演算法(六十五樓)

今天(2007年4月4日)先貼到這裡

[[i] 本帖最後由 simon21 於 2007-4-4 07:39 編輯 [/i]]

simon21 發表於 2006-3-31 08:36

提升法97經典程式

[code]
%% 本程式實現任意偶數大小影象第二代雙正交97提升小波變換 %% 注1: 採用標準正交方法,對行列採用不同矩陣(和matlab裡不同)
%% 注2: 為了保證正交,所有邊界處理,全部採用迴圈處理
%% 注3: 正交性驗證,將單位陣帶入函式,輸出仍是單位陣(matlab不具有此性質)
%% 注4: 此程式是矩陣實現,所以影象水平分量和垂直分量估計被交換位置
%% 注5: 此程式實現的是類小波(wavelet-like)變換,是介於小波包變換與小波變換之間的變換
%% 注6: 此程式每層變換相對原影象矩陣,產生的矩陣都是正交陣,這和小波包一致
%% 注7: 但小波變換每層產生的矩陣,是相對每個待分解子塊的正交矩陣,而不是原影象的正交矩陣
%% 注8: 且小波變換產生的正交矩陣維數,隨分解層數2分減少
%% 注9: 提升係數可以在MATLAB7.0以上版本,用liftwave('9.7')獲取,這裡直接給出,考慮相容性
%% 注10:由於MATLAB陣列下標從1開始,所以注意奇偶序列的變化
%% 注11:d為對偶上升,即預測;p為原上升,即更新 %% 程式設計人 沙威 安徽大學
%% 程式設計時間 2004年12月18日 %% x輸入影象,y輸出影象
%% flag_trans為正變換或反變換標誌,0執行正變換,1執行反變換
%% flag_max,是否最大層數變換標誌,0執行使用者設定層數,1執行最大層數變換
%% layer,使用者層數設定(小於最大層) function y=db97(x,flag_trans,flag_max,layer); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 1.輸入引數檢查 % 矩陣維數判斷
[sa,sb]=size(x); if (sa~=sb) % 防止非影象資料
errordlg('非影象資料!');
error('非影象資料!');
end; % 變換標誌判斷
[sa,sb]=size(flag_trans);
if ((sa~=1) | (sb~=1)) % 變換標誌錯誤
errordlg('變換標誌錯誤!');
error('變換標誌錯誤!');
end; if ((flag_trans~=1) & (flag_trans~=0)) % 變換標誌錯誤
errordlg('變換標誌錯誤!');
error('變換標誌錯誤!');
end; % 最大層數標誌判斷
[sa,sb]=size(flag_max);
if ((sa~=1) | (sb~=1)) % 最大層數標誌錯誤
errordlg('最大層數標誌錯誤!');
error('最大層數標誌錯誤!');
end; if ((flag_max~=1) & (flag_max~=0)) % 最大層數標誌錯誤
errordlg('最大層數標誌錯誤!');
error('最大層數標誌錯誤!');
end; % 使用者設定層數判斷
if (flag_max~=1) [sa,sb]=size(layer);
if ((sa~=1) | (sb~=1)) % 層數設定錯誤
errordlg('層數設定錯誤!');
error('層數設定錯誤!');
end; if (flag_max<0) % 層數設定錯誤
errordlg('層數設定錯誤!');
error('層數設定錯誤!');
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 2.提升係數確定
% t1=liftwave('9.7'); % 獲取提升係數(MATLAB7.0以後) d1=[-1.586100000000000e+000,-1.586134342069360e+000];
p1=[1.079600000000000e+000,-5.298011857188560e-002];
d2=[-8.829110755411875e-001,-8.829110755411875e-001];
p2=[4.435068520511142e-001,1.576123746148364e+000];
d3=-8.698644516247808e-001;
p3=-1.149604398860242e+000;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 3.分解層數確定
% 採用使用者輸入和自動給出最大層數兩種方法 N=length(x); % 矩陣大小
S=N; % 變數
s=log2(N); % 最大迴圈次數
n1=N/2; % 初始一半矩陣大小
n2=N; % 初始矩陣大小
u=0; % 初始值 % 對非2的整數冪大小影象確定最大分解層數
for ss=1:s
if (mod(S,2)==0)
u=u+1;
S=S/2;
end;
end;
u=u-1; % 分解最大層數減1(後面的邊界處理造成) % 最大層數確定
if (flag_max==0) % 手動輸入
T=layer; % 使用者輸入值
else % 自動確定最大層數
T=u; % 分解最大層數
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 4.最大層數和影象大小檢查 if (T>u) % 防止使用者層數越界
errordlg('已超過最大分解層數!或者非偶數大小影象!');
error('已超過最大分解層數!或者非偶數大小影象!');
end; if (mod(N,2)~=0) % 防止影象大小錯誤
errordlg('非偶數大小影象!');
error('非偶數大小影象!');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 5.提升法正變換 if (flag_trans==0)
for time=1:T; % 行正變換

% d;
x1(n1,:)=x(n2,:)+d1(2)*x(n2-1,:)+d1(1)*x(1,:);
x1([1:n1-1],:)=x([2:2:n2-2],:)+d1(2)*x([1:2:n2-3],:)+d1(1)*x([3:2:n2-1],:);

% p;
x(1,:)=x(1,:)+p1(2)*x1(n1,:)+p1(1)*x1(1,:);
x([2:n1],:)=x([3:2:n2-1],:)+p1(2)*x1([1:n1-1],:)+p1(1)*x1([2:n1],:);
x([n1+1:n2],:)=x1([1:n1],:);

% d;
x(n1+1,:)=x(n1+1,:)+d2(2)*x(n1,:)+d2(1)*x(1,:);
x([n1+2:n2],:)=x([n1+2:n2],:)+d2(2)*x([1:n1-1],:)+d2(1)*x([2:n1],:);

% p;
x(n1,:)=x(n1,:)+p2(2)*x(n1+1,:)+p2(1)*x(n1+2,:);
x(n1-1,:)=x(n1-1,:)+p2(2)*x(n2,:)+p2(1)*x(n1+1,:);
x([1:n1-2],:)=x([1:n1-2],:)+p2(2)*x([n1+2:n2-1],:)+p2(1)*x([n1+3:n2],:);

% 歸一
x([1:n1],:)=p3*x([1:n1],:);
x([n1+1:n2],:)=d3*x([n1+1:n2],:); clear x1;

% 列正變換

% d;
x1(:,[1:n1])=x(:,[2:2:n2]);

% p;
x(:,1)=x(:,1)-d1(1)*x1(:,n1)-d1(2)*x1(:,1);
x(:,[2:n1])=x(:,[3:2:n2-1])-d1(1)*x1(:,[1:n1-1])-d1(2)*x1(:,[2:n1]);
x(:,[n1+1:n2])=x1(:,[1:n1]);

% d;
x(:,n2)=x(:,n2)-p1(1)*x(:,n1)-p1(2)*x(:,1);
x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])-p1(1)*x(:,[1:n1-1])-p1(2)*x(:,[2:n1]);

% p;
x(:,n1,:)=x(:,n1)-d2(1)*x(:,n2)-d2(2)*x(:,n1+1);
x(:,[1:n1-1])=x(:,[1:n1-1])-d2(1)*x(:,[n1+1:n2-1])-d2(2)*x(:,[n1+2:n2]);

% d;
x(:,n1+1)=x(:,n1+1)-p2(1)*x(:,n1-1)-p2(2)*x(:,n1);
x(:,n1+2)=x(:,n1+2)-p2(1)*x(:,n1)-p2(2)*x(:,1);
x(:,[n1+3:n2])=x(:,[n1+3:n2])-p2(1)*x(:,[1:n1-2])-p2(2)*x(:,[2:n1-1]);

% 歸一
x(:,[1:n1])=d3*x(:,[1:n1]);
x(:,[n1+1:n2])=p3*x(:,[n1+1:n2]); clear x1;

n2=n2/2; % 原大小
n1=n2/2; % 一半大小
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 6.提升法反變換 else
n2=N/(2.^(T-1)); % 分解最小子塊維數
n1=n2/2;
for time=1:T; % 行反變換

% 去歸一
x([1:n1],:)=x([1:n1],:)/p3;
x([n1+1:n2],:)=x([n1+1:n2],:)/d3; % 反p;
x(n1,:)=x(n1,:)-p2(2)*x(n1+1,:)-p2(1)*x(n1+2,:);
x(n1-1,:)=x(n1-1,:)-p2(2)*x(n2,:)-p2(1)*x(n1+1,:);
x([1:n1-2],:)=x([1:n1-2],:)-p2(2)*x([n1+2:n2-1],:)-p2(1)*x([n1+3:n2],:);

% 反d;
x(n1+1,:)=x(n1+1,:)-d2(2)*x(n1,:)-d2(1)*x(1,:);
x([n1+2:n2],:)=x([n1+2:n2],:)-d2(2)*x([1:n1-1],:)-d2(1)*x([2:n1],:);

% 反p;
x1(1,:)=x(1,:)-p1(2)*x(n2,:)-p1(1)*x(n1+1,:);
x1([2:n1],:)=x([2:n1],:)-p1(2)*x([n1+1:n2-1],:)-p1(1)*x([n1+2:n2],:);

% 反d;
x(n2,:)=x(n2,:)-d1(2)*x1(n1,:)-d1(1)*x1(1,:);
x([2:2:n2-2],:)=x([n1+1:n2-1],:)-d1(2)*x1([1:n1-1],:)-d1(1)*x1([2:n1],:);

% 偶數
x([1:2:n2-1],:)=x1([1:n1],:);

clear x1;

% 列反變換

% 歸一
x(:,[1:n1])=x(:,[1:n1])/d3;
x(:,[n1+1:n2])=x(:,[n1+1:n2])/p3; % 反d;
x(:,n1+1)=x(:,n1+1)+p2(1)*x(:,n1-1)+p2(2)*x(:,n1);
x(:,n1+2)=x(:,n1+2)+p2(1)*x(:,n1)+p2(2)*x(:,1);
x(:,[n1+3:n2])=x(:,[n1+3:n2])+p2(1)*x(:,[1:n1-2])+p2(2)*x(:,[2:n1-1]);

% 反p;
x(:,n1,:)=x(:,n1)+d2(1)*x(:,n2)+d2(2)*x(:,n1+1);
x(:,[1:n1-1])=x(:,[1:n1-1])+d2(1)*x(:,[n1+1:n2-1])+d2(2)*x(:,[n1+2:n2]);

% 反d;
x(:,n2)=x(:,n2)+p1(1)*x(:,n1)+p1(2)*x(:,1);
x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])+p1(1)*x(:,[1:n1-1])+p1(2)*x(:,[2:n1]);

% 反p;
x1(:,1)=x(:,1)+d1(1)*x(:,n2)+d1(2)*x(:,n1+1);
x1(:,[2:n1])=x(:,[2:n1])+d1(1)*x(:,[n1+1:n2-1])+d1(2)*x(:,[n1+2:n2]); % 奇偶
x(:,[2:2:n2])=x(:,[n1+1:n2]);
x(:,[1:2:n2-1])=x1(:,[1:n1]); clear x1;

n2=n2*2; % 原大小
n1=n2/2; % 一半大小 end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 7.結果輸出 y=x;
% 傳輸最後結果 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 8.記憶體清理 clear x;
clear flag_max;
clear layer;
clear flag_trans;
clear N;
clear n1;
clear n2;
clear s;
clear ss;
clear u;
clear d1;
clear d2;
clear d3;
clear p1;
clear p2;
clear p3;
clear sa;
clear sb;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[/code]

[[i] 本帖最後由 yejet 於 2006-8-31 20:32 編輯 [/i]]

simon21 發表於 2006-3-31 08:38

2代小波示意程式

[code]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 此程式用提升法實現第二代小波變換
%% 我用的是非整數階小波變換
%% 採用時域實現,步驟先列後行
%% 正變換:分裂,預測,更新;
%% 反變換:更新,預測,合併
%% 只做一層(可以多層,而且每層的預測和更新方程不同) clear;clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%

% 1.調原始影象矩陣 load wbarb; % 下載影象
f=X; % 原始影象
% f=[0 0 0 0 0 0 0 0 ;...
% 0 0 0 1 1 0 0 0 ;...
% 0 0 2 4 4 2 0 0 ;...
% 0 1 4 8 8 4 1 0 ;...
% 0 1 4 8 8 4 1 0 ;...
% 0 0 2 4 4 2 0 0 ;...
% 0 0 0 1 1 0 0 0 ;...
% 0 0 0 0 0 0 0 0 ;]; % 原始影象矩陣 N=length(f); % 影象維數
T=N/2;

% 子影象維數 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%正變換%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 1.列變換

% A.分裂(奇偶分開) f1=f([1:2:N-1],:); % 奇數
f2=f([2:2:N],:); % 偶數 % f1(:,T+1)=f1(:,1); % 補列
% f2(T+1,:)=f2(1,:); % 補行 % B.預測 for i_hc=1:T;
high_frequency_column(i_hc,:)=f1(i_hc,:)-f2(i_hc,:);
end; % high_frequency_column(T+1,:)=high_frequency_column(1,:); % 補行 % C.更新 for i_lc=1:T;
low_frequency_column(i_lc,:)=f2(i_lc,:)+1/2*high_frequency_column(i_lc,:);
end; % D.合併
f_column([1:1:T],:)=low_frequency_column([1:T],:);
f_column([T+1:1:N],:)=high_frequency_column([1:T],:);


figure(1)
colormap(map);
image(f); figure(2)
colormap(map);
image(f_column);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 2.行變換

% A.分裂(奇偶分開) f1=f_column(:,[1:2:N-1]); % 奇數
f2=f_column(:,[2:2:N]); % 偶數
% f2(:,T+1)=f2(:,1); % 補行 % B.預測 for i_hr=1:T;
high_frequency_row(:,i_hr)=f1(:,i_hr)-f2(:,i_hr);
end; % high_frequency_row(:,T+1)=high_frequency_row(:,1); % 補行 % C.更新 for i_lr=1:T;
low_frequency_row(:,i_lr)=f2(:,i_lr)+1/2*high_frequency_row(:,i_lr);
end; % D.合併
f_row(:,[1:1:T])=low_frequency_row(:,[1:T]);
f_row(:,[T+1:1:N])=high_frequency_row(:,[1:T]);
figure(3)
colormap(map);
image(f_row);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%反變換%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% 1.行變換
% A.提取(低頻高頻分開) f1=f_row(:,[T+1:1:N]); % 奇數
f2=f_row(:,[1:1:T]); % 偶數
% f2(:,T+1)=f2(:,1); % 補行 % B.更新 for i_lr=1:T;
low_frequency_row(:,i_lr)=f2(:,i_lr)-1/2*f1(:,i_lr);
end; % C.預測 for i_hr=1:T;
high_frequency_row(:,i_hr)=f1(:,i_hr)+low_frequency_row(:,i_hr);
end; % high_frequency_row(:,T+1)=high_frequency_row(:,1); % 補行
% D.合併(奇偶分開合併)
f_row(:,[2:2:N])=low_frequency_row(:,[1:T]);
f_row(:,[1:2:N-1])=high_frequency_row(:,[1:T]);
figure(4)
colormap(map);
image(f_row);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% 2.列變換

% A.提取(低頻高頻分開) f1=f_row([T+1:1:N],:); % 奇數
f2=f_row([1:1:T],:); % 偶數 % f1(:,T+1)=f1(:,1); % 補列
% f2(T+1,:)=f2(1,:); % 補行 % B.更新 for i_lc=1:T;
low_frequency_column(i_lc,:)=f2(i_lc,:)-1/2*f1(i_lc,:);
end; % C.預測 for i_hc=1:T;
high_frequency_column(i_hc,:)=f1(i_hc,:)+low_frequency_column(i_hc,:);
end; % high_frequency_column(T+1,:)=high_frequency_column(1,:); % 補行 % D.合併(奇偶分開合併)
f_column([2:2:N],:)=low_frequency_column([1:T],:);
f_column([1:2:N-1],:)=high_frequency_column([1:T],:);


figure(5)
colormap(map);
image(f_column);
[/code]

[[i] 本帖最後由 yejet 於 2006-8-31 20:32 編輯 [/i]]

simon21 發表於 2006-3-31 08:40

二代小波漫談

現在我就舉例,對一個8點序列,怎樣實現第二代小波變換。

1. 奇偶分開。
非常簡單,就是[2,4,6,8]組成一列向量,[1,3,5,7]組成一列向量。

2. 預測。
用[2,4,6,8]來預測[1,3,5,7]。比如說1,3估計2; 3,5估計4; 5,7估計6; 7,1估計8。(邊緣處理,我採用迴圈方法)。估計公式可以用別人的,也可以自己做。舉一個線性的例子:2=1*a+3*b,4=3*a+ 5*b,...,其他的都一樣。這樣我們就可找到最優的a,b,使得(2-(1*a+3*b)).^2+(4-(3*a+5*b)).^2+...最小化。就是最小均方準則。若正好為零,說明偶可以完全預測奇,也就是我們只要儲存偶數列向量,和a,b就可以了,壓縮也就是實現了。對於訊號很長序列,就等於壓縮了一半。當然,我們可以採用更復雜的立方差值預測,多項式預測,或其它的準則,來使其最小,這樣我們的壓縮也就得到了最優。

3. 提升。
我們總希望,均方為零,但可望不可及。於是,提升就需要了。我們經過預測後,要儲存的是偶數序列[2,4,6,8],新的奇數序列[n1,n3,n5, n7]=[2-(1*a+3*b),4-(3*a+5*b),...]和線性變換系數(a,b)。這裡新的奇數序列就是高頻分量。但偶數序列是不能完全代表訊號的性質的,有所差距。所以我們要對偶數序列進行修正。即所謂的提升。我們這次用個簡單的提升吧。[n2,n4,n6,n8]=[2,4,6,8]+ k*[n1,n3,n5,n7]。[n2,n4,n6,n8],就是要分解的低頻分量。那k怎麼求呢?因為要保持n2,n4,n6,n8和原始訊號 [1,2,3,4,5,6,7,8]一樣的性質。一般就是均值和高階矩。這裡就一個未知數k,所以用均值相等,就行了。1/8*(1+2+3+..8)= 1/4*(n2+n4+n6+n8)。k很容易就求出來了。我們最終儲存的就是[n1,n3,n5,n7]和[n2,n4,n6,n8]以及a,b,k。

現在,所謂的第二代就完了。再說幾句。
1.反變換,就是3->2->1。

2.二維。先行提升,再列提升。(我置頂的貼子裡有harr二維提升的原始碼)。

3.整數階。就是加一個取整。

4.多層或小波包提升,就是在對序列[n1,n3,n5,n7]或[n2,n4,n6,n8],再做1->2->3。

5.靈活。不一定是a,b,也可能就一個a,或a,b,c;不一定是一個k,也可能是k1,k2。但越多計算量太大。最好是用大師們做好的CDF,5/3,7/9等。

6.最重要的,任何一代小波,總可以通過一次或多次提升實現。它和一代小波沒有本質區別。

7.優勢。文獻都有,我隨便談談。時域實現,最優壓縮,無邊緣效應,靈活多變,無失真壓縮,程式設計方便,速度快。

文章寫完了,希望對大家有幫助。最主要的,動手編,不要依賴MATLABM,這樣才有所體會。希望和大家多交流。

給 simon21 加一點人氣

[[i] 本帖最後由 simon21 於 2007-4-4 07:05 編輯 [/i]]

simon21 發表於 2006-3-31 08:45

小波濾波器構造和消噪程式(2個)

1.重構

[code]% mallet_wavelet.m

% 此函式用於研究Mallet演算法及濾波器設計

% 此函式僅用於消噪

a=pi/8; %角度賦初值

b=pi/8;

%低通重構FIR濾波器h0(n)衝激響應賦值

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; %確定h1(n)=(-1)^n,高通重建濾波器

if(mod(i_high,2)==0);

coefficient=-1;

else

coefficient=1;

end

high_construct(1,i_high)=low_decompose(1,i_high)*coefficient;

end

high_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);

title('原訊號');

check1=sum(high_decompose); %h0(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); %抽取,得低頻細節

h_fre=conv(y,high_decompose);

h_fre_down=dyaddown(h_fre); %訊號高頻細節

figure(2);

subplot(2,1,1)

plot(l_fre_down);

title('小波分解的低頻係數');

subplot(2,1,2);

plot(h_fre_down);

title('小波分解的高頻係數');

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); %取結果的中心部分,消除卷積影響

h_fre_keep=wkeep(h_fre_denoise,L_signal);

sig_denoise=l_fre_keep+h_fre_keep; %訊號重構

compare=sig_denoise-y; %與原訊號比較

figure(3);

subplot(3,1,1)

plot(y);

ylabel('y'); %原訊號

subplot(3,1,2);

plot(sig_denoise);

ylabel('sig/_denoise'); %重構訊號

subplot(3,1,3);

plot(compare);

ylabel('compare'); %原訊號與消噪後訊號的比較[/code]


2.消噪

[quote]% mallet_wavelet.m

% 此函式用於研究Mallet演算法及濾波器設計

% 此函式用於消噪處理

%角度賦值

%此處賦值使濾波器係數恰為db9

%分解的高頻係數採用db9較好,即它的消失矩較大

%分解的有用訊號小波高頻係數基本趨於零

%對於噪聲訊號高頻分解係數很大,便於閾值消噪處理

[l,h]=wfilters('db10','d');

low_construct=l;

L_fre=20; %濾波器長度

low_decompose=low_construct(end:-1:1); %確定h0(-n),低通分解濾波器

for i_high=1:L_fre; %確定h1(n)=(-1)^n,高通重建濾波器

if(mod(i_high,2)==0);

coefficient=-1;

else

coefficient=1;

end

high_construct(1,i_high)=low_decompose(1,i_high)*coefficient;

end

high_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);

zero1=zeros(1,60); %訊號加噪聲訊號產生

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),性質校驗

check2=sum(low_decompose);

check3=norm(high_decompose);

check4=norm(low_decompose);

l_fre=conv(y_noise,low_decompose); %卷積

l_fre_down=dyaddown(l_fre); %抽取,得低頻細節

h_fre=conv(y_noise,high_decompose);

h_fre_down=dyaddown(h_fre); %訊號高頻細節

figure(2);

subplot(2,1,1)

plot(l_fre_down);

title('小波分解的低頻係數');

subplot(2,1,2);

plot(h_fre_down);

title('小波分解的高頻係數');

% 消噪處理

for i_decrease=31:44;

if abs(h_fre_down(1,i_decrease))>=0.000001

h_fre_down(1,i_decrease)=(10^-7);

end

end

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); %取結果的中心部分,消除卷積影響

h_fre_keep=wkeep(h_fre_denoise,L_signal);

sig_denoise=l_fre_keep+h_fre_keep; %消噪後訊號重構

%平滑處理

for j=1:2

for 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(sig_denoise);

ylabel('sig/_denoise'); %消噪後訊號

subplot(3,1,3);

plot(compare);

ylabel('compare'); %原訊號與消噪後訊號的比較 [/quote]

[[i] 本帖最後由 simon21 於 2007-4-4 07:04 編輯 [/i]]

simon21 發表於 2006-3-31 08:46

小波譜分析mallat演算法經典程式

[code]clc;clear;
%% 1.正弦波定義
f1=50; % 頻率1
f2=100; % 頻率2
fs=2*(f1+f2); % 取樣頻率
Ts=1/fs; % 取樣間隔
N=120; % 取樣點數
n=1:N;
y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts); % 正弦波混合
figure(1)
plot(y);
title('兩個正弦訊號')
figure(2)
stem(abs(fft(y)));
title('兩訊號頻譜')
%% 2.小波濾波器譜分析
h=wfilters('db30','l'); % 低通
g=wfilters('db30','h'); % 高通
h=[h,zeros(1,N-length(h))]; % 補零(圓周卷積,且增大解析度變於觀察)
g=[g,zeros(1,N-length(g))]; % 補零(圓周卷積,且增大解析度變於觀察)
figure(3);
stem(abs(fft(h)));
title('低通濾波器圖')
figure(4);
stem(abs(fft(g)));
title('高通濾波器圖')
%% 3.MALLET分解演算法(圓周卷積的快速傅立葉變換實現)
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); % 頻譜圖
subplot(2,1,1)
stem(abs(fft(sig1)));
title('分解訊號1頻譜')
subplot(2,1,2)
stem(abs(fft(sig2)));
title('分解訊號2頻譜')
%% 4.MALLET重構演算法
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); % 重構低通
gr=g(end:-1:1); % 重構高通
hr=circshift(hr',1)'; % 位置調整圓周右移一位
gr=circshift(gr',1)'; % 位置調整圓周右移一位
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('重構低頻訊號');
subplot(2,1,2)
plot(real(sig2));
title('重構高頻訊號');
figure(8);
subplot(2,1,1)
stem(abs(fft(sig1)));
title('重構低頻訊號頻譜');
subplot(2,1,2)
stem(abs(fft(sig2)));
title('重構高頻訊號頻譜');
figure(9)
plot(real(sig),'r','linewidth',2);
hold on;
plot(y);
legend('重構訊號','原始訊號')
title('重構訊號與原始訊號比較') [/code]

[[i] 本帖最後由 simon21 於 2007-4-4 07:04 編輯 [/i]]

simon21 發表於 2006-3-31 08:46

2維小波變換經典程式

[code]
% FWT_DB.M;
% 此示意程式用DWT實現二維小波變換
% 程式設計時間2004-4-10,程式設計人沙威
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;clc;
T=256; % 影象維數
SUB_T=T/2; % 子圖維數
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1.調原始影象矩陣
load wbarb; % 下載影象
f=X; % 原始影象
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2.進行二維小波分解
l=wfilters('db10','l'); % db10(消失矩為10)低通分解濾波器衝擊響應(長度為20)
L=T-length(l);
l_zeros=[l,zeros(1,L)]; % 矩陣行數與輸入影象一致,為2的整數冪
h=wfilters('db10','h'); % db10(消失矩為10)高通分解濾波器衝擊響應(長度為20)
h_zeros=[h,zeros(1,L)]; % 矩陣行數與輸入影象一致,為2的整數冪
for i=1:T; % 列變換
row(1:SUB_T,i)=dyaddown( ifft( fft(l_zeros).*fft(f(:,i)') ) ).'; % 圓周卷積<->FFT
row(SUB_T+1:T,i)=dyaddown( ifft( fft(h_zeros).*fft(f(:,i)') ) ).'; % 圓周卷積<->FFT
end;
for j=1:T; % 行變換
line(j,1:SUB_T)=dyaddown( ifft( fft(l_zeros).*fft(row(j,:)) ) ); % 圓周卷積<->FFT
line(j,SUB_T+1:T)=dyaddown( ifft( fft(h_zeros).*fft(row(j,:)) ) ); % 圓周卷積<->FFT
end;
decompose_pic=line; % 分解矩陣
% 影象分為四塊
lt_pic=decompose_pic(1:SUB_T,1:SUB_T); % 在矩陣左上方為低頻分量--fi(x)*fi(y)
rt_pic=decompose_pic(1:SUB_T,SUB_T+1:T); % 矩陣右上為--fi(x)*psi(y)
lb_pic=decompose_pic(SUB_T+1:T,1:SUB_T); % 矩陣左下為--psi(x)*fi(y)
rb_pic=decompose_pic(SUB_T+1:T,SUB_T+1:T); % 右下方為高頻分量--psi(x)*psi(y)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3.分解結果顯示
figure(1);
colormap(map);
subplot(2,1,1);
image(f); % 原始影象
title('original pic');
subplot(2,1,2);
image(abs(decompose_pic)); % 分解後影象
title('decomposed pic');
figure(2);
colormap(map);
subplot(2,2,1);
image(abs(lt_pic)); % 左上方為低頻分量--fi(x)*fi(y)
title('/Phi(x)*/Phi(y)');
subplot(2,2,2);
image(abs(rt_pic)); % 矩陣右上為--fi(x)*psi(y)
title('/Phi(x)*/Psi(y)');
subplot(2,2,3);
image(abs(lb_pic)); % 矩陣左下為--psi(x)*fi(y)
title('/Psi(x)*/Phi(y)');
subplot(2,2,4);
image(abs(rb_pic)); % 右下方為高頻分量--psi(x)*psi(y)
title('/Psi(x)*/Psi(y)');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5.重構源影象及結果顯示
% construct_pic=decompose_matrix'*decompose_pic*decompose_matrix;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l_re=l_zeros(end:-1:1); % 重構低通濾波
l_r=circshift(l_re',1)'; % 位置調整
h_re=h_zeros(end:-1:1); % 重構高通濾波
h_r=circshift(h_re',1)'; % 位置調整

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
top_pic=[lt_pic,rt_pic]; % 影象上半部分
t=0;
for i=1:T; % 行插值低頻

if (mod(i,2)==0)
topll(i,:)=top_pic(t,:); % 偶數行保持
else
t=t+1;
topll(i,:)=zeros(1,T); % 奇數行為零
end
end;
for i=1:T; % 列變換
topcl_re(:,i)=ifft( fft(l_r).*fft(topll(:,i)') )'; % 圓周卷積<->FFT
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bottom_pic=[lb_pic,rb_pic]; % 影象下半部分
t=0;
for i=1:T; % 行插值高頻
if (mod(i,2)==0)
bottomlh(i,:)=bottom_pic(t,:); % 偶數行保持
else
bottomlh(i,:)=zeros(1,T); % 奇數行為零
t=t+1;
end
end;
for i=1:T; % 列變換
bottomch_re(:,i)=ifft( fft(h_r).*fft(bottomlh(:,i)') )'; % 圓周卷積<->FFT
end;

construct1=bottomch_re+topcl_re; % 列變換重構完畢

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
left_pic=construct1(:,1:SUB_T); % 影象左半部分
t=0;
for i=1:T; % 列插值低頻

if (mod(i,2)==0)
leftll(:,i)=left_pic(:,t); % 偶數列保持
else
t=t+1;
leftll(:,i)=zeros(T,1); % 奇數列為零
end
end;
for i=1:T; % 行變換
leftcl_re(i,:)=ifft( fft(l_r).*fft(leftll(i,:)) ); % 圓周卷積<->FFT
end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
right_pic=construct1(:,SUB_T+1:T); % 影象右半部分

t=0;
for i=1:T; % 列插值高頻
if (mod(i,2)==0)
rightlh(:,i)=right_pic(:,t); % 偶數列保持
else
rightlh(:,i)=zeros(T,1); % 奇數列為零
t=t+1;
end
end;
for i=1:T; % 行變換
rightch_re(i,:)=ifft( fft(h_r).*fft(rightlh(i,:)) ); % 圓周卷積<->FFT
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
construct_pic=rightch_re+leftcl_re; % 重建全部影象

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 結果顯示
figure(3);
colormap(map);
subplot(2,1,1);
image(f); % 源影象顯示
title('original pic');

subplot(2,1,2);
image(abs(construct_pic)); % 重構源影象顯示
title('reconstructed pic');
error=abs(construct_pic-f); % 重構圖形與原始影象誤值
figure(4);
mesh(error); % 誤差三維影象
title('absolute error display'); [/code]

[[i] 本帖最後由 simon21 於 2007-4-4 07:05 編輯 [/i]]

eros 發表於 2006-4-2 09:40

請問能不能指教一下二次樣條小波呢?matlab裡面如何實現呢?

simon21 發表於 2006-4-6 14:26

回覆:(simon21)[分享]個人收集的一些關於小波分析的...

基於LeventCodes平臺的小波去噪程式

包括以下方法:

BivaShrink方法、模型1、模型2、模型3(TrivaShrink方法)、BayesShrink方法、 LAWMLShrink方法的DWT實現和DT_CWT實現。

simon21 發表於 2006-4-6 14:30

回覆:(simon21)[分享]個人收集的一些關於小波分析的...

連續小波和離散小波分析的應用例項

simon21 發表於 2006-4-6 14:31

小波插值與小波構造(3個程式)

[code]%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 小波構造 function casade
clear;clc;
t=3;
phi=[0,1,0]; h=wfilters('db7','r');
h=h*sqrt(2); h_e=h(1,[2:2:14]);
h_o=h(1,[1:2:13]); for m=1:15;
stem(phi);
drawnow;
pause(1);
ee=conv(h_e,phi);
oo=conv(h_o,phi);
phi(1,[2:2:2*length(ee)])=ee;
phi(1,[1:2:2*length(oo)-1])=oo;


end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cubic_average(立方b樣條)
% 均值插值 % 初始化
s=[0 0 1 0 0] % 正弦波
% f=50;
% ts=1/200;
% n=0:16;
% s=sin(2*pi*f*n*ts); % 係數
se=[1/8,6/8,1/8];
so=[4/8,4/8] % 迴圈
for p=1:10;
t=length(s)-1;
o(1:t)=s(1:t)*so(1)+s(2:t+1)*so(2);
e(1)=s(t+1)*se(1)+s(1)*se(2)+s(2)*se(3);
e(2:t)=s(1:t-1)*se(1)+s(2:t)*se(2)+s(3:t+1)*se(3);
e(t+1)=s(t)*se(1)+s(t+1)*se(2)+s(1)*se(3);
s([1:2:2*t+1])=e([1:t+1]);
s([2:2:2*t])=o([1:t]);
plot(s);
drawnow;
end; % 抽取
t=length(s); % 總長度
p=128; % 需要點數 % 間隔
d=(t-1)/p; % 最終尺度函式
r=s(2:d:t-1); % 畫圖
figure(2);
plot(r);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cubic_subdivision(立方插值)
% 細分插值 % %% 初始化(尺度函式)
% s=[0,0,1,0,0]; % 正弦函式
n=1:20;
f=50;
ts=1/200;
s=sin(2*pi*f*n*ts); % 指數函式
% n=0:16;
% s=exp(n); % % 係數
a=[-1/16,9/16,9/16,-1/16]; % 迴圈
for p=1:4;
t=length(s)-1;
o(1)=s(4)*a(1)+s(1)*a(2)+s(2)*a(3)+s(3)*a(4);
o(2:t-1)=s(1:t-2)*a(1)+s(2:t-1)*a(2)+s(3:t)*a(3)+s(4:t+1)*a(4);
o(t)=s(t-2)*a(4)+s(t+1)*a(3)+s(t)*a(2)+s(t-1)*a(1);
s([1:2:2*t+1])=s([1:t+1]);
s([2:2:2*t])=o([1:t]);
plot(s);
drawnow;
end; % % 抽取
% t=length(s); % 總長度
% p=128; % 需要點數
%
% % 間隔
% d=(t-1)/p;
%
% % 最終尺度函式
% r=s(2:d:t-1);
%
% % 畫圖
% figure(2);
% plot(r);[/code]

[[i] 本帖最後由 simon21 於 2007-4-4 07:06 編輯 [/i]]

simon21 發表於 2006-4-6 14:33

採用多孔trous演算法(undecimated wavelet transform)實現小波變換

[code]clear;clc;
% 1.生成訊號
f=50; % 頻率
fs=800; % 取樣率
T=128; % 訊號長度
n=1:T;
y=sin(2*pi*f*n/fs)+2*exp(-f*n/(4*fs)); % 訊號
% y=circshift(y.',3).'; %%
2.正變換 l1=wfilters('db4','l')*sqrt(2)/2;
% 參考低通濾波器
l1_zeros=[l1,zeros(1,T-length(l1))]; % 低通濾波器1
h1=wfilters('db4','h')*sqrt(2)/2; % 參考高通濾波器
h1_zeros=[h1,zeros(1,T-length(h1))]; % 高通濾波器1 low1=ifft(fft(y).*fft(l1_zeros)); % 低頻分量1
high1=ifft(fft(y).*fft(h1_zeros)); % 高頻分量1 l2=dyadup(l1); % 原濾波器插值
l2_zeros=[l2,zeros(1,T-length(l2))]; % 低通濾波器2
h2=dyadup(h1); % 原濾波器插值
h2_zeros=[h2,zeros(1,T-length(h2))]; % 高通濾波器2 low2=ifft(fft(low1).*fft(l2_zeros)); % 低頻分量2
high2=ifft(fft(low1).*fft(h2_zeros)); % 高頻分量2
%% 3.反變換 lr2=circshift(l2_zeros(end:-1:1).',1).'; % 重構低通濾波器2
hr2=circshift(h2_zeros(end:-1:1).',1).'; % 重構高通濾波器2
lr1=circshift(l1_zeros(end:-1:1).',1).'; % 重構低通濾波器1
hr1=circshift(h1_zeros(end:-1:1).',1).'; % 重構高通濾波器1 lowr=(ifft(fft(low2).*fft(lr2))+ifft(fft(high2).*fft(hr2))); % 重構低頻分量1(lowr=low1)
r_s=(ifft(fft(lowr).*fft(lr1))+ifft(fft(high1).*fft(hr1))); % 重構源訊號(r_s=y)
%% 4.繪圖 figure(1);
plot(y);
title('源訊號'); figure(2);
plot(low1,'r');
hold on;
plot(low2,'b');
legend('第一層低頻','第二層低頻'); figure(3);
plot(high1,'r');
hold on;
plot(high2,'b');
legend('第一層高頻','第二層高頻'); figure(4);
plot(low1,'r');
hold on;
plot(lowr,'b.');
legend('第一層低頻','重構第一層低頻'); figure(5);
plot(y,'r');
hold on;
plot(r_s,'b.');
legend('源訊號','重構訊號'); disp(norm(low1-lowr))
disp(norm(y-r_s))[/code]

[[i] 本帖最後由 simon21 於 2007-4-4 07:07 編輯 [/i]]

simon21 發表於 2006-4-6 14:35

此程式實現構造小波基

[code]% 此程式實現構造小波基
% periodic_wavelet.m function ss=periodic_wavelet; clear;clc; % global MOMENT; % 消失矩階數
% global LEFT_SCALET; % 尺度函式左支撐區間
% global RIGHT_SCALET; % 尺度函式右支撐區間
% global LEFT_BASIS; % 小波基函式左支撐區間
% global RIGHT_BASIS; % 小波基函式右支撐區間
% global MIN_STEP; % 最小離散步長
% global LEVEL; % 計算需要的層數(離散精度)
% global MAX_LEVEL; % 週期小波最大計算層數
[s2,h]=scale_integer;
[test,h]=scalet_stretch(s2,h);
wave_base=wavelet(test,h);
ss=periodic_waveletbasis(wave_base);
function [s2,h]=scale_integer; % 本函式實現求解小波尺度函式離散整數點的值
% sacle_integer.m MOMENT=10; % 消失矩階數
LEFT_SCALET=0; % 尺度函式左支撐區間
RIGHT_SCALET=2*MOMENT-1; % 尺度函式右支撐區間
LEFT_BASIS=1-MOMENT; % 小波基函式左支撐區間
RIGHT_BASIS=MOMENT; % 小波基函式右支撐區間
MIN_STEP=1/512; % 最小離散步長
LEVEL=-log2(MIN_STEP); % 計算需要的層數(離散精度)
MAX_LEVEL=8; % 週期小波最大計算層數
h=wfilters('db10','r'); % 濾波器係數 h=h*sqrt(2); % FI(T)=SQRT(2)*SUM(H(N)*FI(2T-N)) N=0:2*MOMENT-1; for i=LEFT_SCALET+1:RIGHT_SCALET-1
for j=LEFT_SCALET+1:RIGHT_SCALET-1
k=2*i-j+1;
if (k>=1&k<=RIGHT_SCALET+1)
a(i,j)=h(k); % 矩陣係數矩陣
else
a(i,j)=0;
end
end
end [s,w]=eig(a); % 求特徵向量,解的基
s1=s(:,1);
s2=[0;s1/sum(s1);0]; % 根據條件SUM(FI(T))=1,求解; % 本函式實現尺度函式經伸縮後的離散值
% scalet_stretch.m function [s2,h]=scalet_stretch(s2,h); MOMENT=10; % 消失矩階數
LEFT_SCALET=0; % 尺度函式左支撐區間
RIGHT_SCALET=2*MOMENT-1; % 尺度函式右支撐區間
LEFT_BASIS=1-MOMENT; % 小波基函式左支撐區間
RIGHT_BASIS=MOMENT; % 小波基函式右支撐區間
MIN_STEP=1/512; % 最小離散步長
LEVEL=-log2(MIN_STEP); % 計算需要的層數(離散精度)
MAX_LEVEL=8; % 週期小波最大計算層數
for j=1:LEVEL % 需要計算到尺度函式的層數
t=0;
for i=1:2:2*length(s2)-3 % 需要計算的離散點取值(0,1,2,3 -> 1/2, 3/2, 5/2)
t=t+1;
fi(t)=0;
for n=LEFT_SCALET:RIGHT_SCALET; % 低通濾波器衝擊響應緊支撐判斷
if ((i/2^(j-1)-n)>=LEFT_SCALET&(i/2^(j-1)-n)<=RIGHT_SCALET) % 小波尺度函式緊支撐判斷
fi(t)=fi(t)+h(n+1)*s2(i-n*2^(j-1)+1); % 反覆應用雙尺度方程求解
end
end
end
clear s
n1=length(s2);
n2=length(fi);
for i=1:length(s2)+length(fi) % 變換後的矩陣長度
if (mod(i,2)==1)
s(i)=s2((i+1)/2); % 矩陣奇數下標為小波上一層(0,1,2,3)離散值
else
s(i)=fi(i/2); % 矩陣偶數下標為小波下一層(1/2,3/2,5/2)(經過伸縮變換後)的離散值
end
end
s2=s;
end
% 採用雙尺度方程求解小波基函式 PSI(T)
% wavelet.m function wave_base=wavelet(test,h); MOMENT=10; % 消失矩階數
LEFT_SCALET=0; % 尺度函式左支撐區間
RIGHT_SCALET=2*MOMENT-1; % 尺度函式右支撐區間
LEFT_BASIS=1-MOMENT; % 小波基函式左支撐區間
RIGHT_BASIS=MOMENT; % 小波基函式右支撐區間
MIN_STEP=1/512; % 最小離散步長
LEVEL=-log2(MIN_STEP); % 計算需要的層數(離散精度)
MAX_LEVEL=8; % 週期小波最大計算層數
i=0;
for t=LEFT_BASIS:MIN_STEP:RIGHT_BASIS; % 小波基支撐長度
s=0;
for n=1-RIGHT_SCALET:1-LEFT_SCALET % g(n)取值範圍
if((2*t-n)>=LEFT_SCALET&(2*t-n)<=RIGHT_SCALET) % 尺度函式判斷
s=s+h(1-n+1)*(-1)^(n)*test((2*t-n)/MIN_STEP+1); % 計算任意精度的小波基函式值
end
end
i=i+1;
wave_base(i)=s;
end[/code]

[[i] 本帖最後由 simon21 於 2007-4-4 07:07 編輯 [/i]]

infant 發表於 2006-4-26 16:18

<P>請問樓主有對稱小波的程式嗎?謝謝了</P>[em04]

ruola 發表於 2006-4-30 10:50

請問有沒有小波係數模極大值的程式啊!
頁: [1] 2 3 4

 

Powered by Discuz! X1.5 Archiver   © 2001-2010 Comsenz Inc.

相關文章