前言
本文集中前面主要介紹了離散數(shù)據(jù)的傅里葉變換狼忱，并且得到了較好的效果扭仁！那既然有了傅里葉變換這個(gè)工具讼油，為什么還需要小波變換呢晓折？因?yàn)椋焊道锶~變換只能告訴你原始信號(hào)中有哪些頻率惑朦，但不能告訴你這些頻率的信號(hào)出現(xiàn)在什么時(shí)間！也就說(shuō)明：如果信號(hào)是"時(shí)變"的(頻率隨著時(shí)間是改變的)漓概，那么單純用傅里葉變換所能反映的信息就十分有限了漾月！因此，針對(duì)時(shí)變信號(hào)垛耳，我們使用小波變換栅屏。圖1展示"時(shí)變信號(hào)"與"時(shí)不變信號(hào)"區(qū)別：
時(shí)不變與時(shí)變的區(qū)別飘千，看下面的實(shí)現(xiàn)的代碼就很輕易理解：
x = 0:0.001:1;
% 4個(gè)頻率:
f1 = 50; f2 = 80; f3 = 110; f4 = 20;
% 時(shí)不變信號(hào): 多頻
y1 = sin(2*pi*f1*x) + sin(2*pi*f2*x) + sin(2*pi*f3*x) + sin(2*pi*f4*x);
% 時(shí)變信號(hào): 多頻
y2 = sin(2*pi*f1*x).*(x<=0.3) + sin(2*pi*f2*x).*(x>0.3 & x<=0.6)+...
sin(2*pi*f3*x).*(x>0.6 & x<=0.8) + sin(2*pi*f4*x).*(x>0.8);
figure(1);
subplot(2,1,1);
plot(y1);
grid on;
xlabel('采樣點(diǎn)'); ylabel('振幅');
axis([0 1000 -inf inf]);
title('時(shí)不變信號(hào)');
subplot(2,1,2);
plot(y2);
grid on;
xlabel('采樣點(diǎn)'); ylabel('振幅');
axis([0 1000 -inf inf]);
title('時(shí)變信號(hào)')
本文同樣考慮的是離散數(shù)據(jù)的小波變換使用。通過(guò)手動(dòng)matlab編程實(shí)現(xiàn)小波變換"塔式分解"與"重構(gòu)"來(lái)深刻了解小波變換實(shí)現(xiàn)的內(nèi)在含義栈雳。之后护奈，借助matlab自帶的一系列相關(guān)小波變換程序來(lái)實(shí)現(xiàn)"時(shí)頻分析"和"小波去噪"。
說(shuō)明：本文更加側(cè)重詳細(xì)介紹matlab自帶各種小波功能函數(shù)的使用哥纫！除了小波分解與重構(gòu)的程序我們手動(dòng)實(shí)現(xiàn)外霉旗，其他的各種操作都建議用自帶函數(shù)實(shí)現(xiàn)。
小波分解：
小波分解的流程總結(jié)為：先將信號(hào)對(duì)半分解成"低頻近似"與"高頻細(xì)節(jié)"2個(gè)部分蛀骇；同樣的操作每次將上一次的"低頻近似"部分再分成低頻近似和高頻細(xì)節(jié)部分厌秒，逐次細(xì)分(最多分解到每個(gè)部分只有1個(gè)點(diǎn))。每次分出的高頻細(xì)節(jié)部分不做分解擅憔。因此：每次分出低頻近似部分相當(dāng)于對(duì)本次信號(hào)做"低通濾波"鸵闪，分出的高頻細(xì)節(jié)部分相當(dāng)于對(duì)本次信號(hào)做"高通濾波"。所以：每次小波分解就是用1個(gè)低通濾波器和1個(gè)高通濾波器對(duì)本次信號(hào)做1次低通濾波和1次高通濾波而已暑诸。
由上述說(shuō)明可得：小波分解的關(guān)鍵在于2個(gè)(一組)濾波器蚌讼。對(duì)于現(xiàn)實(shí)的離散數(shù)據(jù)而言，濾波器看上去很高大上其實(shí)就是很簡(jiǎn)單的數(shù)字而已个榕，濾波聽(tīng)起來(lái)很難篡石，其實(shí)就是做"點(diǎn)乘相加"而已。離散小波分解中最簡(jiǎn)單的一組濾波器為：
低通濾波器：西采；高通濾波器：
凰萨。
下面舉例說(shuō)明如何用上面這一組最簡(jiǎn)單濾波器對(duì)離散數(shù)據(jù)進(jìn)行小波分解：
假設(shè)我們的離散數(shù)據(jù)為：
(1) 第一級(jí)分解：
低通濾波：
得到低頻近似部分：
高通濾波：
得到高頻細(xì)節(jié)部分：
（2）第二級(jí)分解：
分解上一次分解得到的低頻近似部分，高頻細(xì)節(jié)不動(dòng)械馆。即待分解信號(hào)為：
低通濾波：
得到的低頻近似部分為：
高通濾波：
得到的高頻細(xì)節(jié)部分為：
（3）第3級(jí)分解：
同理胖眷，待分解信號(hào)是上一次分解得到的低頻近似部分：
低通濾波：
得到的低頻近似部分為：
高通濾波：
得到的高頻細(xì)節(jié)部分為：
到此為止，例子中給出的原始離散信號(hào)就達(dá)到了其最大分解級(jí)數(shù)(分解到元素只有1個(gè))狱杰。整個(gè)過(guò)程很簡(jiǎn)單瘦材，只需注意每次點(diǎn)乘求和元素是無(wú)重合的。整個(gè)的多級(jí)分解過(guò)程如圖2所示：
注意：不同組的高通和低通濾波中都有這樣的一個(gè)規(guī)律：兩者的區(qū)別只是高通濾波器中第2個(gè)值是負(fù)數(shù)而已仿畸；數(shù)都是一樣的食棕。
小波多級(jí)分解清楚了，那怎么"重構(gòu)/恢復(fù)"回去呢错沽？塔式分解的逆向合成而已簿晓。根據(jù)濾波器的規(guī)律，我們可以設(shè)：
以2級(jí)分解到3級(jí)為例千埃，我們知道：
那么逆過(guò)程就是憔儿，我們知道了3級(jí)低通近似和高通細(xì)節(jié)的值，我們還知道濾波器的數(shù)值(a已知)放可，然后反推2級(jí)低通近似和高通細(xì)節(jié)數(shù)值谒臼，即：
所以：
這樣我們就得到了2級(jí)低通近似和2解高通細(xì)節(jié)朝刊，然后逐級(jí)往上遞推即可實(shí)現(xiàn)重構(gòu)/恢復(fù)~ so easy.
整個(gè)分解過(guò)程我們清楚了，現(xiàn)在我們引入一些專(zhuān)業(yè)的名詞：在離散數(shù)據(jù)中蜈缤，一組低通高通濾波器拾氓，其實(shí)就是"小波基函數(shù)"！取不同的小波基函數(shù)其實(shí)就是濾波器里面的數(shù)值不同而已底哥。最常用的"haar小波基"咙鞍。下面我們就利用haar小波基，在matlab里手動(dòng)實(shí)現(xiàn)小波分解與重構(gòu)：
matlab手動(dòng)實(shí)現(xiàn)小波分解程序：
clc ; clear;
% 每次修改這里的原始數(shù)據(jù), 個(gè)數(shù)最好是2^n
% x = [9 7 3 5];
x = [2 5 8 9 7 4 -1 1];
% x = [2 5 8 9 7 4 -1 1 2 1 8 3 8 0 3 1];
order_max = log(length(x))/log(2);
fprintf('當(dāng)前數(shù)據(jù)最多分解%d階\n',order_max);
order = double(input('自定義分解階數(shù)( order<order_max ):'));
% matlab默認(rèn)的haar小波基函數(shù):
low = [1/sqrt(2) 1/sqrt(2)];
high = [1/sqrt(2) -1/sqrt(2)];
new = zeros(1,length(x)); % 最后結(jié)果 —— 均值 + 差值
ave = zeros(1,length(x)/2); % 均值(低頻)記錄
dec = zeros(1,length(x)/2); % 差值(高頻)記錄
% 小波循環(huán)分解部分: 其實(shí)就是低通和高通2種情況的卷積計(jì)算
m = 1;
xtmp = x;
for norder = 1:order
for n = 1:2:length(xtmp)
ave(m) = sum(xtmp(n:n+1).*low);
dec(m) = sum(xtmp(n:n+1).*high);
m = m + 1;
end
% 下面2句的賦值過(guò)程, 總體是從后往前賦值的~
% 進(jìn)入過(guò)的高頻就不會(huì)動(dòng)了; 進(jìn)入的低頻再下次就會(huì)被自己分解出的高頻和低頻取代——總體從后往前√
new( length(xtmp)/2+1:length(xtmp) ) = dec( 1:2^(order_max-norder) ); % 記錄每次更新的高頻內(nèi)容;
new( 1:length(xtmp)/2 ) = ave( 1:2^(order_max-norder) ); % 記錄每次更新的低頻內(nèi)容趾徽；
xtmp = ave( 1:2^(order_max-norder) );
m = 1;
end
fprintf('手寫(xiě)%d級(jí)分解結(jié)果為:\n',order);
new
fprintf('matlab自帶%d級(jí)分解結(jié)果為:\n',order);
wavedec(x,order,'haar')
matlab手動(dòng)實(shí)現(xiàn)小波重構(gòu)/恢復(fù)程序：
clc ; clear;
% 每次修改這里的原始數(shù)據(jù), 個(gè)數(shù)最好是2^n
% x = [9 7 3 5];
x = [2 5 8 9 7 4 -1 1];
% x = [2 5 8 9 7 4 -1 1 2 1 8 3 8 0 3 1];
order_max = log(length(x))/log(2);
fprintf('當(dāng)前數(shù)據(jù)最多分解%d階\n',order_max);
order = double(input('自定義分解階數(shù)( order<order_max ):'));
% 直接用知自帶的分解函數(shù)：
new = wavedec(x,order,'haar');
% 小波重構(gòu)——任意haar基函數(shù)全能重構(gòu)回去
% 計(jì)算的系數(shù): 與基函數(shù)有關(guān)
tmp1 = 1/( low(1) + high(1) );
% 迭代重構(gòu)開(kāi)始:
xrec = zeros(1,length(new)); % 記錄復(fù)原的數(shù)據(jù)
new_rec = new;
for norder = order_max+1-order : order_max % 這個(gè)規(guī)律讓任意低階都可以適用
m = 1; % 專(zhuān)門(mén)用來(lái)記錄"前半段"位置——給奇數(shù)位用的
for n = 1:2^norder
half = 2^norder/2; % 當(dāng)前重構(gòu)區(qū)間的一半——也是用來(lái)給奇數(shù)位計(jì)算用的
if mod(n,2) ~= 0
% 奇數(shù)時(shí)操作:
xrec(n) = tmp1*( new_rec(m) + new_rec(m+half) );
m = m + 1;
else
% 偶數(shù)時(shí)操作:
xrec(n) = (new_rec(m-1) - low(1)*xrec(n-1))/low(2);
end
end
new_rec(1:n) = xrec(1:n); % 每次要更新的(后一次重構(gòu)基于前一次結(jié)果): 從前往后
end
fprintf('%d級(jí)分解后重構(gòu):\n',order);
new_rec
fprintf('自帶的%d級(jí)重構(gòu)結(jié)構(gòu):\n',order);
[C,S] = wavedec(x,order,'haar');
waverec(C,S,'haar')
