wt.rar_小波变换wt_小波变换去噪_时频去噪

function [tdst_at_certain_point,tdst_at_certain_freq]=tdst(image) %2013年1月2号完成，在2G内存情况下，可以处理84*84的图像，更大的图像无法处理，会出现out of memory的情况 %说明，由于kx=0，ky=0处的值远大于其他点的值，因此默认将其抑制，以方便图像的显示 %用到了fft2,用到了元胞矩阵，也可以用结构体来储存，不过貌似不如用元胞矩阵 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ %输入参数：image：待变换图像，不能是彩色图像，也就是说必须是一个二维矩阵 %输出参数：tdst_at_certain_point：各点处的S变换结果。是一个元胞矩阵，行列数与图像矩阵相同 % tdst_at_certain_freq：各频率处对整幅图像进行S变换的矩阵，行列数比图像矩阵多1 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ image=imread('lenagray.bmp'); tic; TRUE=1; FALSE=0; factor_x=1; %x方向上分辨率控制因子，通常取1，其值增加可以提高分辨率 factor_y=1; %y方向上分辨率控制因子，通常取1，其值增加可以提高分辨率 suppress_kx0ky0=TRUE; %抑制kx=0，ky=0处的值 [N,M]=size(image); fft2_image=fft2(image); %对图像进行傅里叶变换，为后面移位做准备 tdst_at_certain_freq=cell(N+1,M+1); %初始化元胞矩阵，每个元素是一个矩阵，代表着在某频率点分析整幅图像的结果 kx0ky0=mean(mean(image)); %kx=0,ky=0处二维S变换 %计算二维S变换矩阵 for kx=0:N for ky=0:M tdst_at_certain_freq{kx+1,ky+1}=ones(N,M); %初始化该频率点处的变换矩阵 if (kx==0&&ky==0) %零频处的值为图像特点均值 if suppress_kx0ky0==1 %检查是否选择抑制kx=0，ky=0处的值 tdst_at_certain_freq{kx+1,ky+1}=kx0ky0*zeros(N,M); %抑制该点的值 elseif suppress_kx0ky0==0 tdst_at_certain_freq{kx+1,ky+1}=kx0ky0*ones(N,M); %正确生成该点的值 else disp(sprintf('Error.Please choose whether to suppress the value at kx=0,ky=0')); end else tdst_at_certain_freq{kx+1,ky+1}=ifft2(circshift(fft2_image,[-kx-1,-ky-1]).*gauss(N,M,kx+1,ky+1,factor_x,factor_y));%移位是用circshift完成的，x是竖直向下的方向，y是水平向右的方向。移位后与高斯函数窗相乘，再进行傅里叶反变换 end end end %上面得到的变换结果是在各频率点处的整幅图像的变换值，下面将其转换为在各左边点处整个频率平面上的变换值 tdst_at_certain_point=cell(N,M); %初始化矩阵 for x=1:N for y=1:M tdst_at_certain_point{x,y}=extract(tdst_at_certain_freq,x,y); %这里无法避免使用循环，对某个坐标点，从tdst_at_certain_freq中抽取相应元素，得到该点在各频率处的变换值 end end time=toc; disp(sprintf('Estimated time is %f',time)); function gauss_w=gauss(length_x,length_y,freq_x,freq_y,factor_x,factor_y) %二维高斯函数窗 %输入参数：length_x：x方向上长度，即图像矩阵的行数 % length_y：y方向上长度，即图像矩阵的列数 % freq_x：需要计算高斯函数点处x方向上的频率值 % freq_y：需要计算高斯函数点处y方向上的频率值 % factor_x:x方向上分辨率控制因子 % factor_y:y方向上分辨率控制因子 %输出参数：gauss_w：二维高斯函数，是一个矩阵 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ if freq_x~=0 %频率为0时，高斯函数需要单独处理 vector_x=[0:length_x-1;-length_x:-1]; %定义行,这里使用两行来计算高斯函数，是参照Stockwell给出的一维S变换程序，具体原因不详，但是这样可以获得更好的分辨率 vector_x=factor_x*vector_x.^2/freq_x^2; vector_x=sum(vector_x); %两行相加 vector_x=vector_x'; %得到列 matrix_x=vector_x*ones(1,length_y); %把列变为矩阵 else matrix_x=zeros(length_x,length_y); %频率为0时，没有该部分 end if freq_y~=0 %频率为0时，高斯函数需要单独处理 vector_y=[0:length_y-1;-length_y:-1]; %定义行 vector_y=factor_y*vector_y.^2/freq_y^2; vector_y=sum(vector_y); %两行相加 matrix_y=ones(length_x,1)*vector_y; %把行变成矩阵 else matrix_y=zeros(length_x,length_y); %频率为0时，没有该部分 end matrix=matrix_x+matrix_y; %得到高斯函数系数雏形 gauss_w=matrix*(-2)*pi^2; %整个二维高斯函数窗的指数 gauss_w=exp(gauss_w); %二维高斯函数窗 function matrix=extract(matrix_to_be_extracted,x,y) %将tdst_at_certain_freq变成tdst_at_certain_point用到的抽取函数，从tdst_at_certain_freq各矩阵元素的同一位置抽取元素按次序组成tdst_at_certain_point的某一矩阵元素 %输入参数：matrix_to_be_extracted:待抽取的元胞矩阵 % x：需要抽取的位置的坐标x % y：需要抽取的位置的坐标y %输出参数：matrix：由输入矩阵抽取相应元素得到的新矩阵 %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ [N M]=size(matrix_to_be_extracted); %确定输入矩阵的行列数，即确定在x，y方向上共有多少频率点 matrix=ones(N,M); %初始化矩阵 for j=1:N for k=1:M temp=matrix_to_be_extracted{j,k}; %找出输入矩阵某一位置的矩阵 matrix(j,k)=temp(x,y); %将该矩阵的元素按输入的要求取出并赋给matrix end end