matlab求图像频谱

clc;clear;close all;
ori = imread('A.bmp');
% [h, w,ch] = size(ori);
I = im2double(ori);

%% 用matlab函数实现求频谱
F = fft2(I(:,:,1)); % 针对R通道做fft变换
F1 = log(abs(F)+1); % abs取绝对值让复数变实数，log压缩数据范围，一定要压缩，不然数据太大了！！
F2 = log(abs(fftshift(F))+1);
figure; imshow(F1,[]);  % 没有经过中心化处理得到的频谱图
figure; imshow(F2,[]);  % 经过中心化处理得到的频谱图

%% 根据原理自行实现求频谱过程
function [Wu,Wv] = trans(M,N,flag)
Wu=zeros(M,M);Wv=zeros(N,N);
for u = 0:M-1
    for x = 0:M-1
        if flag ==1
            Wu(u+1,x+1) = exp(-1j*2*pi*u*x/M);%FFT
        else
            Wu(u+1,x+1) = exp(1j*2*pi*u*x/M);%IFFT
        end
    end
endfor v = 0:N-1
    for y = 0:N-1
        if flag ==1
            Wv(v+1,y+1) = exp(-1j*2*pi*v*y/N);%FFT
        else
            Wv(v+1,y+1) = exp(1j*2*pi*v*y/N);%IFFT
        end
    end
end

%原理部分自行补充
[Wu,Wv] = trans(h,w,1); % 调用trans函数，见前文。当flag==1，做FFT变换，当flag==0，做IFFT变换I_trans=zeros(h,w,1); % 针对一通道做的处理，多通道时将1改成ch，并且下面的循环中要多加一层循环k=1：ch，I_trans(i,j,1)改为I_trans(i,j,k)
for i = 1:h
	for j = 1:w
		I_trans(i,j,1) = I(i,j,1)*(-1)^(i+j); % 使得变换的频谱中心化，如果没有这一步，得到的是没有中心化的频谱
	end
end
I_trans1 = Wu*I_trans(:,:,1)*Wv;
I_trans2 = log(abs(I_trans1)+1);
figure; imshow(I_trans2,[]); %中心化频谱

两种方法求得的频谱图有细微出入，原因不明，但是都可以使用。