DCT及FFT2变换解释

A = [73 82 73 91 164 156 74 83 74 164 165 165 165 94 91 91 166 166 29 21 101 91 165 102 111 93 95 22 102 30 29 94 83 101 175 84 92 82 103 82 10 82 82 12 21 102 92 74 0 82 82 82 10 83 73 73 91 82 82 91 82 91 90 91 ];

B = dct2(A);

% High pass filter

hpfilter = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 1 1 1 1 ];

C = B.*hpfilter;

D = idct2(C);

figure(1);

subplot(1,2,1); imshow(A, [0 255]); subplot(1,2,2); imshow(D, [0 255]);

figure(2); subplot(1,2,1); imshow(abs(B), [0 255]); figure(2); subplot(1,2,2); imshow(abs(C), [0 255]);   

使用matlab自带dct2函数进行DCT变换，然后对频谱使用高通滤波，然后将滤波结果进行DCT反变换 

Orignal DCTFiltered DCT Orignal imageFiltered image200100 0  

1 

2001000 

clc; clear;

A = [73 82 73 91 164 156 74 83 74 164 165 165 165 94 91 91 166 166 29 21 101 91 165 102 111 93 95 22 102 30 29 94 83 101 175 84 92 82 103 82 10 82 82 12 21 102 92 74 0 82 82 82 10 83 73 73 91 82 82 91 82 91 90 91 ];

Const = [1/sqrt(2) 1 1 1 1 1 1 1]; M = 8; for u=1:M for v=1:M temp = 0; for x=1:M

for y=1:M % DCT formula use x,y,u,v=0...N-1, so cos() use (x-1)

temp = temp + A(x,y)*Const(u)*Const(v)*cos( (2*(x-1)+1)*(u-1)*pi/(2*M) )

*cos( (2*(y-1)+1)*(v-1)*pi/(2*M) );

end end

F(u,v) = temp*2/sqrt(M*M); end end

B = dct2(A);   

使用DCT的定义式计算DCT 

710.50 121.82 31.31 -45.22 8.00 -102.02 -21.09 44.83

5.10 30.27 -39.37 -47.28 18.38 -27.05 -44.25 -11.02 2.71 -29.77 -77.28 -76.72 -26.83 44.88 85.50 22.50 -9.47 55.88 19.29 -32.76 14.52 37.10 97.59 26.52 -59.75 47.41 20.74 -53.39 59.25 13.89 -15.90 -21.37 -54.44 -48.44 6.12 0.71 -5.62 -20.95 -4.03 33.97 -9.78 -30.43 -5.25 50.91 49.54 8.09 -16.47 -22.43 21.15 32.09 -56.29 -18.54 39.33 -19.36 3.65 -6.55

2 

RGB = imread('autumn.tif'); I = rgb2gray(RGB); J = dct2(I);

imshow(log(abs(J)),[]), colormap(jet(64)), colorbar

1050-5  

RGB = imread('autumn.tif'); I = rgb2gray(RGB); J = dct2(I);

imshow(log(abs(J)),[]), colormap(gray), colorbar

1050-5  

RGB = imread('autumn.tif'); I = rgb2gray(RGB); J = dct2(I);

imshow((abs(J)),[]), colormap(jet(64)), colorbar

x 102.521.510.5 4 

不加log由于数据太大，动态范围大，所以不能完全显示。   

3 

clc; clear;

I = imread('cameraman.tif'); D = fft2(I);

figure(1);imshow(log(abs(D)), []), colorbar title('FFT2 without shift') DF = fftshift(D);

figure(2);subplot(1,2,1);

imshow(log(abs(DF)), []), colorbar title('FFT2 with shift')

D = dct2(I); subplot(1,2,2);

imshow(log(abs(D)), []), colorbar title('DCT2')

FFT2 without shift 25020015010050 510 15FFT2 with shift 15105 DCT2 1050-5 

4 

clc; clear;

x = zeros(256,1); n=100; nw=1; for i=n:n+nw x(i) = 1; end

f = fft(x);

figure(1); subplot(2,1,1); plot(x); axis([0 300 -1 2]); title(['Signal, width=', num2str(nw)]);

figure(1); subplot(2,1,2); plot(abs(f(1:128))) title('Spectrum');

Signal, width=121.510.50-0.5-1050100150Spectrum21.510.500200250300Signal, width=1021.510.50-0.5-1050100150Spectrum12108642200250300204060801001201400020406080100120140 

Signal, width=3021.510.50-0.5-1050100150Spectrum403020200250300Signal, width=8021.510.50-0.5-1050100150Spectrum100806040200250300100020002040608010012014020406080100120140矩形脉冲宽度窄时，高频分量多，宽度宽时，高频分量减少 

5 

clc; clear;

A = zeros(256, 256); A = im2uint8(A); m=100; n=100; mw=50; nw=1; for i=m:m+mw for j=n:n+nw

A(i,j) = 255; end end

figure(1); subplot(2,1,1); imshow(A), colorbar

title(['Image, pulse width=', num2str(nw), ', pulse height=', num2str(mw)]); B = fft2(A); C = fftshift(B);

figure(1);subplot(2,1,2); imshow((abs(C)), []), colorbar title('Image spectrum');  

Image, pulse width=1, pulse height=30 25020015010050 Image spectrum 15000100000Image, pulse width=10, pulse height=30 25020015010050 Image spectrum 04x 108645000 02  

Image, pulse width=30, pulse height=30 25020015010050 Image spectrum 05Image, pulse width=80, pulse height=30 25020015010050 Image spectrum 05x 1021.510.5x 10642 只是考察了水平方向的变换，脉冲水平宽度窄时，水平高频分量多，水平宽度宽时，水平高频分量减少 垂直方向可以类推。    

6 

目前还不明白斜对角的频谱表示什么含义？ clc; clear;

A = zeros(256, 256); A = im2uint8(A); m=100; n=100; r=1;

for i=1:256 for j=1:256

if( sqrt( (i-m)^2+(j-n)^2 ) <= r ) A(i,j) = 255; end end end

figure(1); subplot(2,1,1); imshow(A), colorbar

title(['Image, pulse center=', num2str(n), ', pulse radius=', num2str(r)]); B = fft2(A); C = fftshift(B);

figure(1);subplot(2,1,2); imshow((abs(C)), []), colorbar title('Image spectrum');  

Image, pulse center=100, pulse radius=1 200Image, pulse center=100, pulse radius=3 200100100 Image spectrum 0 Image spectrum012001000800600400200 600040002000 Image, pulse center=100, pulse radius=5 200Image, pulse center=100, pulse radius=8 200100100 Image spectrum 04 Image spectrum 04x 1021.510.5x 1054321  

可以看出，FFT2虽然只是在水平和垂直方向的频谱分析，但是也显示了斜方向的频谱情况，如果使用一个椭圆脉冲，应该可以更清楚的看出在斜方向的频谱情况。 

7 

Image, pulse center=100, pulse radius=10 Image, pulse center=100, pulse radius=30200100 0Image spectrum x48 10642 8 

200100 0Image spectrumx 105 642