基於奇異值分解的影象壓縮與除噪
一、本報告所用的一些基本原理
第一，慮噪過程。
巴特沃斯濾波器的特點是通頻帶內的頻率響應曲線最大限度平坦，沒有起伏，而在阻頻帶則逐漸下降為零。在振幅的對數對角頻率的波得圖上，從某一邊界角頻率開始，振幅隨著角頻率的增加而逐步減少，趨向負無窮大。
本程式通過不斷調整通帶截止頻率、阻帶下限截止頻率、通帶衰減及阻帶衰減，力圖使濾噪後的影象逼近未加噪影象。
第二，影象壓縮思路。
影象資訊的資料量非常的大，隨著各種成像裝置的解析度不斷提高，單幅影象所包含的資料量也越來越大，大資料量的影象資訊會給儲存器的儲存容量、通訊通道的頻寬以及計算機的處理速度增加極大的壓力。為了解決這個問題，必須對影象進行壓縮處理
在本次數值分析的影象壓縮課程模擬實踐當中，原圖是一個2000*3000的一個矩陣，這裡就是一個10^6的一個很大的資料量，我們就要想辦法在保證影象質量的前提下，儘可能的壓縮影象，在實踐當中分別選擇了1/10、1/100、1/1000、1/10000進行練習。
我們小組在壓縮實踐中主要運用的思路是基於SVD演算法的影象壓縮，在進行1/10的壓縮時，我們對2000*3000的矩陣進行20*30的分塊，每一小塊為一個100*100的矩陣，進行奇異值分解時，我們只選取了前5個奇異值，這樣在傳輸過程中的資料量僅為K=100*5+5*5+5*100=1025，K/（100*100）近似等於1/10，這樣在影象的傳輸過程中，就只需對U（:,1:5)、S(1:5,1:5)、V(:,1:5)'三個矩陣的資料進行傳輸，而不必在傳輸以前的矩陣資料，資料接收方只需對這三個矩陣進行處理，即可得到壓縮了1/10後的影象。
在進行1/100的壓縮時，對濾波後的影象矩陣整體進行奇異值分解，這一次取前12個奇異值，這樣的資料量K=2000*12+12*12+12*3000=60144，K/(2000*3000)近似於1/100。
這樣的思路在進行1/1000及更高比例的壓縮時遇到了問題，因為要進行1/1000的壓縮，只能取1個奇異值左後才能夠滿足壓縮要求，但是這樣壓縮之後的影象卻嚴重的失真了，原因當然是清楚的，對於2000個奇異值只取最大的1個奇異值，影象的很多東西就被丟棄了，這樣影象壓縮遇到了瓶頸，不得不停滯下來思考其他方法；經過小組成員之間的討論和與其他小組的交流中，有了新的思路，先將2000*3000的濾波後矩陣進行200*300的分塊，每一小塊有100個數據，通過c=mean（CB（:))函式對每一小塊的矩陣進行均值化處理，這樣一個10*10的矩陣變成了一個均值c，再通過新矩陣賦值CD(i,j)=c,最後一個2000*3000的矩陣通過均值化處理再賦值變成了200*300的矩陣，在保持了影象的基本質量中影象已經壓縮了1/100，後面只需對這200*300的矩陣壓縮1/10即可達到壓縮1/1000的要求，在顯示回影象的過程中，因為2000*3000的矩陣變成了200*300的矩陣，因此x、y座標的長度的縮小十分之一（這個問題一直沒有注意到，導致壓縮程式老是出現錯誤，顯示不出影象，後再小組成員的討論中終於尋找出了錯誤）
最後是進行1/10000的壓縮，也就是對上一步中的那個200*300的矩陣壓縮1/100，進行奇異值分解，只取前2個奇異值，這樣K=200*2+2*2+2*300=1004，K/（200*300）近似於1/60，已和1/100處在同一個數量級，能夠滿足要求。
二、關於影象噪點過濾、壓縮及誤差所編寫程式
clear all ;
clc
%% 生成2000*3000資料
f = @(x,y) exp(x-2*x.^2-y.^2).*sin(6*(x + y +x.*y.^2));
x =linspace(-1,1,3e3);
y =linspace(-1.2,1.2,2e3); y = y(:) ;
A = f(repmat(x,2e3,1),repmat(y,1,3e3)) ;
AA=A;
r = -0.2 + (0.2+0.2).*rand(2000,3000); %增加噪聲項
A=A+r;
%% 原始資料圖
figure,mesh(x,y,A); %三維網格圖
xlabel('原始加噪點圖');
view([-40,80]); %調整視角
%% Butterworth對影象進行濾波
disp('１．Butterworth Filtering');
%低通濾波器技術要求，經上分析，擬定
%通帶截止頻率為1.5Hz，阻帶下限截止頻率為4Hz
%通帶衰減為0.7dB，阻帶衰減為1.3dB
fs=3000;Wp=2*pi*1.5/fs;Ws=2*pi*4/fs;Rp=0.7;Rs=1.3;
Omip=Wp/pi;Omis=Ws/pi; %歸一化技術要求
[N,Wn]=buttord(Omip,Omis,Rp,Rs); %確定濾波器的階數
disp(['The order of Butterworth Filtering is',num2str(N)]);
[b,a]=butter(N,Wn); %確定Butterworth濾波器轉移函式係數向量
R=filter(b,a,A);
figure,mesh(x,y,R); %三維網格圖,濾波後圖像
xlabel('濾波後的影象');
view([-40,80]); %調整視角
VC=AA(:); %將矩陣變成列向量
VC0=R(:); %將矩陣變成列向量
RMSE0=sqrt(sum((VC-VC0).^2)/6000000); %求均方根誤差
%% 分塊奇異值壓縮1/10
BB=R;
for i=1:1:20
for j=1:1:30
BC=(BB(((i-1)*100+1):i*100,((j-1)*100+1):j*100));
[U,S,V]=svd(BC); %對每個100*100塊做奇異值分解
BD=U(:,1:5)*S(1:5,1:5)*V(:,1:5)'; %取前5個特徵值
BK(((i-1)*100+1):i*100,((j-1)*100+1):j*100)=BD;
end
end
%% 分塊奇異值分解後的圖
figure,mesh(x,y,BK);
xlabel('壓縮1/10後的影象');
view([-40,80]);
VC1=BK(:); %將矩陣變成列向量
RMSE1=sqrt(sum((VC-VC1).^2)/6000000); %求均方根誤差
%% 奇異值分解進行壓縮1/100
[U,S,V]=svd(R); %奇異值分解出U,S,V三個矩陣
C = U(:,1:12)*S(1:12,1:12)*V(:,1:12)'; %取U的前12列，S的前12個奇異值，V的前12列
U2= U(:,1:12);
S2= S(1:12,1:12);
V2= V(:,1:12)';
%% 奇異值分解後的圖
figure,mesh(x,y,C);
xlabel('壓縮1/100後的影象');
view([-40,80]);
VC2=C(:); %將矩陣變成列向量
RMSE2=sqrt(sum((VC-VC2).^2)/6000000); %求均方根誤差
%% 分塊奇異值10*10均值後壓縮1/1000
CA=R;
for i=1:1:200
for j=1:1:300
CB=(CA(((i-1)*10+1):i*10,((j-1)*10+1):j*10));
c=mean(CB(:));%全部平均
CD(i,j)=c;
end
end
%% 分塊奇異值壓縮1/10
CE=CD;
for i=1:1:2
for j=1:1:3
CF=(CE(((i-1)*100+1):i*100,((j-1)*100+1):j*100));
[U,S,V]=svd(CF); %對每個10*10塊做奇異值分解
CG=U(:,1:5)*S(1:5,1:5)*V(:,1:5)'; %取前5個特徵值
CH(((i-1)*100+1):i*100,((j-1)*100+1):j*100)=CG;
end
end
%% 分塊奇異值壓縮後的圖
x =linspace(-1,1,3e2);
y =linspace(-1.2,1.2,2e2); y = y(:) ;
figure,mesh(x,y,CH);
xlabel('壓縮1/1000後的影象');
view([-40,80]);
%% 壓縮後的CD奇異值分解進行壓縮1/10000
E=CD;
[U,S,V]=svd(E); %奇異值分解出U,S,V三個矩陣
D = U(:,1:2)*S(1:2,1:2)*V(:,1:2)'; %取U的前2列，S的前2個奇異值，V的前2列
U4= U(:,1:2);
S4= S(1:2,1:2);
V4= V(:,1:2)';
%% 奇異值分解後的圖
figure,mesh(x,y,D);
xlabel('壓縮1/10000後的影象');
view([-40,80]);
三、程式執行後的結果見附件二
評價：對影象進行壓縮，不可避免的要引入失真。我們要做的就是在影象訊號的終端使用者察覺不出或能夠忍受這些是真的前提下，進一步提高壓縮比，以換取更高的編碼效率，這就需要一些標準來對影象壓縮的質量進行評價。在我們的程式中，主要運用的是均方根誤差來對影象壓縮的質量進行評價。
四、參考文獻及資料
【1】ShoichiroNakamura .科學計算引論.電子工業出版社，2002.
【2】低通濾波技術在數字影象處理中的應用以及模擬實現（百度資料）.
【3】語音訊號的濾波處理（課程設計）.
【4】胡鄉峰，衛金茂. 基於奇異值分解的影象壓縮.東北師範大學學報，2006.
【5】曾超，張衛東.基於奇異值分解的影象壓縮及其Matlab實現.高校理科研究，2010.