hw1.zip_dct_dct idct_iff

function [] = Untitled2() I1=imread('I1.jpg'); I2=imread('I2.jpg'); I3=imread('I3.jpg'); I1=rgb2gray(I1); I2=rgb2gray(I2); I3=rgb2gray(I3); %对图像2,3进行fft变换，然后交换二者相位值求逆fft变换 i2=my_fft2(I2); i3=my_fft2(I3); i2_amp=abs(i2); i3_amp=abs(i3); i2_phase=angle(i2); i3_phase=angle(i3); I_2=my_ifft2(i2_amp.*exp(i*i3_phase)); I_3=my_ifft2(i3_amp.*exp(i*i2_phase)); %对图像1进行dct变换，并分别求DCT系数保留前五个，前十个和前二十时其逆dct变换的结果 I_dct_inital=my_dct2(I1); [I_dct_5,I_dct_10,I_dct_20]=deal(zeros(size(I1))); I_dct_5=dct_filter(I_dct_inital,3); I_dct_10=dct_filter(I_dct_inital,4); I_dct_20=dct_filter(I_dct_inital,6); I_inital=my_idct2(I_dct_inital)./255; I_5=my_idct2(I_dct_5)./255; I_10=my_idct2(I_dct_10)./255; I_20=my_idct2(I_dct_20)./255; X=dct2(I1); x=idct2(X)./255; imwrite(I_inital,'I_inital.jpg'); imwrite(I_5,'I_5.jpg'); imwrite(I_10,'I_10.jpg'); imwrite(I_20,'I_20.jpg'); %交换2和3的相位后恢复出的图像 figure(1); subplot(2,2,1);imshow(I2);title('人像'); subplot(2,2,2);imshow(I3);title('飞机'); subplot(2,2,3);colormap(gray);imagesc(abs(I_2));title('交换相位后人像'); subplot(2,2,4);colormap(gray);imagesc(abs(I_3));title('交换相位后飞机'); %对比自定义fft2函数所得频谱和matlab内置fft2函数所得频谱 figure(2); i1=my_fft2(I1); %自定义fft2 sfftI1=fftshift(i1); %直流分量移到频谱中心 A1=abs(sfftI1); A1=(A1-min(min(A1)))/(max(max(A1))-min(min(A1)))*225;%归一化 i1_matlab=fft2(I1); %MATLAB内置fft2 s_matlab=fftshift(i1_matlab); A1_matlab=abs(s_matlab); A1_matlab=(A1_matlab-min(min(A1_matlab)))/(max(max(A1_matlab))-min(min(A1_matlab)))*225;%归一化 subplot(1,2,1);imshow(A1);title('自定义fft所得频谱'); subplot(1,2,2);imshow(A1_matlab);title('内置fft所得频谱'); %自定义dct和内置dct恢复图像对比 figure(3); subplot(1,2,1);imshow(I_inital);title('自定义dct和idct'); subplot(1,2,2);imshow(x);title('内置dct和idct'); %保留不同比例的dct系数后恢复出的图像 figure(4); subplot(2,2,1);imshow(I1);title('原图'); subplot(2,2,2);imshow(I_5);title('保留前5个dct系数'); subplot(2,2,3);imshow(I_10);title('保留前10个dct系数'); subplot(2,2,4);imshow(I_20);title('保留前20个dct系数'); %计算峰值信噪比 MES5=sum(sum((double(I1)-double(I_5)).^2))/(512*512) PSNR5=20*log10(255/sqrt(MES5)) MES10=sum(sum((double(I1)-double(I_10)).^2))/(512*512) PSNR10=20*log10(255/sqrt(MES10)) MES20=sum(sum((double(I1)-double(I_20)).^2))/(512*512) PSNR20=20*log10(255/sqrt(MES20)) end function I_fft=my_fft2(I) %2d-fft sizex=size(I); row=sizex(1);col=sizex(2); [U,V]=deal(zeros(row,col)); for p=0:row-1 for q=0:row-1 U(p+1,q+1)=exp(-2*i*pi*p*q/row); end end for p=0:col-1 for q=0:col-1 V(p+1,q+1)=exp(-2*i*pi*p*q/col); end end I=double(I); I_fft=U'*I*V; end function I_ifft=my_ifft2(I) %2d-ifft sizex=size(I); row=sizex(1);col=sizex(2); [U,V]=deal(zeros(row,col)); for p=0:row-1 for q=0:row-1 U(p+1,q+1)=exp(-2*i*pi*p*q/row); end end for p=0:col-1 for q=0:col-1 V(p+1,q+1)=exp(-2*i*pi*p*q/col); end end I=double(I); I_ifft=inv(U')*I*inv(V); end function I_dct=my_dct2(I) %2d-dct sizex=size(I); row=sizex(1);col=sizex(2); nr=row/8;nc=col/8; I_dct=zeros(row,col); for r=1:nr for c=1:nc I_dct(8*r-7:8*r,8*c-7:8*c)=my_dct8x8(I(8*r-7:8*r,8*c-7:8*c)); end end end function I_dct=my_dct8x8(I) %8x8块大小的2d-dct C=zeros(8,8); for p=0:7 for q=0:7 if p==0 a=sqrt(1/8); else a=sqrt(1/4); end C(p+1,q+1)=a*cos(pi*(q+0.5)*p/8); end end I=double(I); I_dct=C*I*C'; end function I_idct=my_idct2(I) %2d-idct sizex=size(I); row=sizex(1);col=sizex(2); nr=row/8;nc=col/8; I_idct=zeros(row,col); for r=1:nr for c=1:nc I_idct(8*r-7:8*r,8*c-7:8*c)=my_idct8x8(I(8*r-7:8*r,8*c-7:8*c)); end end end function I_idct=my_idct8x8(I) %8x8块大小的2d-idct C=zeros(8,8); for p=0:7 for q=0:7 if p==0 a=sqrt(1/8); else a=sqrt(1/4); end C(p+1,q+1)=a*cos(pi*(q+0.5)*p/8); end end I=double(I); I_idct=C'*I*C; end function I_x=dct_filter(I,number) Isize=size(I); I_x=zeros(Isize); tmp=zeros(8,8); for i=1:number-1 for j=1:number-1 if j+i<=number tmp(i,j)=1; end end end for r=0:Isize(1)/8-1 for c=0:Isize(2)/8-1 I_x(8*r+1:8*r+8,8*c+1:8*c+8)=I(8*r+1:8*r+8,8*c+1:8*c+8).*tmp; end end end