MATLAB 灰度图像的二维傅里叶变换-深圳市維司達科技有限公司

基本傅里叶变换

基础实现代码

functionbasic_fft2_demo()% 读取灰度图像I=imread('cameraman.tif');% 使用MATLAB自带图像或替换为你的图像路径% 转换为double类型以便计算I_double=im2double(I);% 执行二维傅里叶变换F=fft2(I_double);% 将零频率分量移到中心F_shift=fftshift(F);% 计算幅度谱和相位谱magnitude_spectrum=abs(F_shift);phase_spectrum=angle(F_shift);% 对幅度谱取对数以便更好地显示log_magnitude=log(1+magnitude_spectrum);% 显示结果figure('Position',[100,100,1200,300]);subplot(1,4,1);imshow(I);title('原始图像');subplot(1,4,2);imshow(log_magnitude,[]);title('对数幅度谱');subplot(1,4,3);imshow(phase_spectrum,[]);title('相位谱');% 逆变换重建图像F_ishift=ifftshift(F_shift);I_reconstructed=real(ifft2(F_ishift));subplot(1,4,4);imshow(I_reconstructed,[]);title('重建图像');% 计算重建误差reconstruction_error=norm(I_double-I_reconstructed,'fro');fprintf('重建误差: %e\n',reconstruction_error);end

频域滤波应用

频域滤波器实现

functionfrequency_filtering_demo()% 读取图像I=imread('cameraman.tif');I=im2double(I);[M,N]=size(I);% 傅里叶变换F=fft2(I);F_shift=fftshift(F);% 创建频率坐标网格u=(0:M-1)-floor(M/2);v=(0:N-1)-floor(N/2);[U,V]=meshgrid(v,u);D=sqrt(U.^2+V.^2);% 到中心的距离% 设计不同的滤波器cutoff_low=30;% 低通截止频率cutoff_high=30;% 高通截止频率cutoff_band=50;% 带通中心频率% 低通滤波器 (理想低通)H_low=double(D<=cutoff_low);% 高通滤波器 (理想高通)H_high=double(D>=cutoff_high);% 高斯低通滤波器D0=30;% 截止频率H_gaussian=exp(-(D.^2)/(2*D0^2));% 应用滤波器F_low=F_shift.*H_low;F_high=F_shift.*H_high;F_gaussian=F_shift.*H_gaussian;% 逆变换I_low=real(ifft2(ifftshift(F_low)));I_high=real(ifft2(ifftshift(F_high)));I_gaussian=real(ifft2(ifftshift(F_gaussian)));% 显示结果figure('Position',[100,100,1200,600]);subplot(2,4,1);imshow(I);title('原始图像');subplot(2,4,2);imshow(H_low,[]);title('理想低通滤波器');subplot(2,4,3);imshow(H_high,[]);title('理想高通滤波器');subplot(2,4,4);imshow(H_gaussian,[]);title('高斯低通滤波器');subplot(2,4,5);imshow(I,[]);title('原始图像');subplot(2,4,6);imshow(I_low,[]);title('理想低通滤波结果');subplot(2,4,7);imshow(I_high,[]);title('理想高通滤波结果');subplot(2,4,8);imshow(I_gaussian,[]);title('高斯低通滤波结果');end

图像增强应用

高频增强和同态滤波

functionimage_enhancement_demo()% 读取图像I=imread('pout.tif');% 使用低对比度图像I=im2double(I);[M,N]=size(I);% 傅里叶变换F=fft2(I);F_shift=fftshift(F);% 高频增强滤波器u=(0:M-1)-floor(M/2);v=(0:N-1)-floor(N/2);[U,V]=meshgrid(v,u);D=sqrt(U.^2+V.^2);D_max=max(D(:));% 高频增强传递函数alpha=1.5;% 增强系数beta=0.5;H_enhance=alpha+beta*(D/D_max);% 应用高频增强F_enhanced=F_shift.*H_enhance;I_enhanced=real(ifft2(ifftshift(F_enhanced)));% 同态滤波 (用于同时增强对比度和压缩动态范围)% 对图像取对数I_log=log(I+0.01);% 加小值避免log(0)% 傅里叶变换F_log=fft2(I_log);F_log_shift=fftshift(F_log);% 同态滤波函数 (高频增强，低频抑制)gammaL=0.5;% 低频增益gammaH=2.0;% 高频增益c=1;% 锐化常数D0=30;% 截止频率H_homomorphic=(gammaH-gammaL)*(1-exp(-c*(D.^2/D0^2)))+gammaL;% 应用同态滤波F_homomorphic=F_log_shift.*H_homomorphic;I_homomorphic=real(ifft2(ifftshift(F_homomorphic)));I_homomorphic=exp(I_homomorphic);% 指数变换% 显示结果figure('Position',[100,100,1200,400]);subplot(1,4,1);imshow(I);title('原始图像');subplot(1,4,2);imshow(H_enhance,[]);title('高频增强滤波器');subplot(1,4,3);imshow(I_enhanced,[]);title('高频增强结果');subplot(1,4,4);imshow(I_homomorphic,[]);title('同态滤波结果');end

纹理分析和特征提取

频谱分析和特征计算

functiontexture_analysis_demo()% 读取纹理图像I=imread('texture.png');% 替换为你的纹理图像I=im2double(I);ifsize(I,3)==3I=rgb2gray(I);end[M,N]=size(I);% 傅里叶变换F=fft2(I);F_shift=fftshift(F);magnitude_spectrum=abs(F_shift);% 计算径向功率谱center_x=floor(N/2)+1;center_y=floor(M/2)+1;[X,Y]=meshgrid(1:N,1:M);R=sqrt((X-center_x).^2+(Y-center_y).^2);R=round(R);% 计算径向平均功率谱r_max=min(floor(M/2),floor(N/2));radial_profile=zeros(1,r_max);forr=0:r_max-1mask=(R==r);ifsum(mask(:))>0radial_profile(r+1)=mean(magnitude_spectrum(mask));endend% 计算纹理特征total_power=sum(magnitude_spectrum(:).^2);% 低频能量比例 (中心区域)low_freq_mask=(R<=min(M,N)/8);low_freq_energy=sum(magnitude_spectrum(low_freq_mask).^2)/total_power;% 高频能量比例 (外围区域)high_freq_mask=(R>=min(M,N)/4);high_freq_energy=sum(magnitude_spectrum(high_freq_mask).^2)/total_power;% 显示结果figure('Position',[100,100,1000,400]);subplot(1,3,1);imshow(I);title('原始纹理图像');subplot(1,3,2);imshow(log(1+magnitude_spectrum),[]);title('对数幅度谱');subplot(1,3,3);plot(0:r_max-1,radial_profile);xlabel('径向频率');ylabel('平均幅度');title('径向功率谱');grid on;% 输出纹理特征fprintf('纹理特征分析:\n');fprintf('总功率: %.4f\n',total_power);fprintf('低频能量比例: %.4f\n',low_freq_energy);fprintf('高频能量比例: %.4f\n',high_freq_energy);fprintf('低频/高频能量比: %.4f\n',low_freq_energy/high_freq_energy);end

噪声检测和图像恢复

周期性噪声检测和去除

functionnoise_detection_removal()% 创建带有周期性噪声的图像I_clean=im2double(imread('cameraman.tif'));[M,N]=size(I_clean);% 添加周期性噪声[X,Y]=meshgrid(1:N,1:M);noise=0.1*sin(2*pi*X/20)+0.1*sin(2*pi*Y/15);I_noisy=I_clean+noise;I_noisy=min(max(I_noisy,0),1);% 限制在[0,1]范围内% 傅里叶变换F_noisy=fft2(I_noisy);F_shift=fftshift(F_noisy);magnitude_spectrum=abs(F_shift);log_magnitude=log(1+magnitude_spectrum);% 检测噪声频率 (寻找明显的峰值)threshold=mean(log_magnitude(:))+2*std(log_magnitude(:));noise_mask=(log_magnitude>threshold)&...~((X>=N/2-5&X<=N/2+5)&(Y>=M/2-5&Y<=M/2+5));% 排除DC分量% 创建陷波滤波器H_notch=ones(M,N);H_notch(noise_mask)=0;% 应用陷波滤波器F_filtered=F_shift.*H_notch;I_filtered=real(ifft2(ifftshift(F_filtered)));% 显示结果figure('Position',[100,100,1200,400]);subplot(1,4,1);imshow(I_noisy);title('带噪声图像');subplot(1,4,2);imshow(log_magnitude,[]);title('噪声图像频谱');hold on;[noise_y,noise_x]=find(noise_mask);plot(noise_x,noise_y,'r.','MarkerSize',10);title('频谱 (红点:检测到的噪声)');subplot(1,4,3);imshow(H_notch,[]);title('陷波滤波器');subplot(1,4,4);imshow(I_filtered,[]);title('去噪后图像');% 计算信噪比改善snr_original=10*log10(var(I_clean(:))/var(I_noisy(:)-I_clean(:)));snr_filtered=10*log10(var(I_clean(:))/var(I_filtered(:)-I_clean(:)));fprintf('信噪比分析:\n');fprintf('原始图像SNR: %.2f dB\n',snr_original);fprintf('滤波后SNR: %.2f dB\n',snr_filtered);fprintf('SNR改善: %.2f dB\n',snr_filtered-snr_original);end

实用工具函数

通用的傅里叶变换工具包

functionfft_toolkit()% 这是一个综合的傅里叶变换工具包% 1. 图像傅里叶变换可视化functionvisualize_fft2(I,title_str)I=im2double(I);F=fft2(I);F_shift=fftshift(F);magnitude=abs(F_shift);phase=angle(F_shift);figure('Position',[100,100,800,200]);subplot(1,4,1);imshow(I);title([title_str,' - 原图']);subplot(1,4,2);imshow(log(1+magnitude),[]);title('幅度谱');subplot(1,4,3);imshow(phase,[]);title('相位谱');% 重建验证I_recon=real(ifft2(ifftshift(F_shift)));subplot(1,4,4);imshow(I_recon,[]);title('重建图像');end% 2. 频域滤波器设计functionH=create_filter(type,M,N,params)u=(0:M-1)-floor(M/2);v=(0:N-1)-floor(N/2);[U,V]=meshgrid(v,u);D=sqrt(U.^2+V.^2);switchtypecase'lowpass'D0=params.D0;H=double(D<=D0);case'highpass'D0=params.D0;H=double(D>=D0);case'gaussian_lowpass'D0=params.D0;H=exp(-(D.^2)/(2*D0^2));case'butterworth_lowpass'D0=params.D0;n=params.n;H=1./(1+(D./D0).^(2*n));otherwiseH=ones(M,N);endend% 3. 频域卷积functionI_filtered=frequency_convolution(I,H)I=im2double(I);F=fft2(I);F_shift=fftshift(F);F_filtered=F_shift.*H;I_filtered=real(ifft2(ifftshift(F_filtered)));end% 使用示例I=imread('cameraman.tif');visualize_fft2(I,'测试图像');[M,N]=size(I);params.D0=30;H=create_filter('gaussian_lowpass',M,N,params);I_filtered=frequency_convolution(I,H);figure;subplot(1,2,1);imshow(I);title('原图');subplot(1,2,2);imshow(I_filtered,[]);title('高斯低通滤波');end