傅里叶变换、理想低通、巴特沃斯低通、高斯低通滤波器的matlab
2019-11-26 本文已影响0人
glider312
1D傅里叶变换
1.画出一个有噪声的一维信号：
clc
clear all
close all
fs = 1000; % Sampling frequency
T = 1/fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(x)); % Sinusoids plus noise
plot(fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
2.用fft()函数对信号进行傅里叶变换，画出频域图像
figure(2);
f1_index=(0:L-1)*(fs/L);
y1=abs(fft(x)).*2/L;
plot(f1_index,y1);
xlabel('频率(Hz)');ylabel('幅度(V)');title('fft sin');
最终结果：
image.png
image.png
2D傅里叶变换
1.构造一个有噪声的二维信号：
clc
clear all
close all
mRows = 1024;
x = linspace(-1, 1, mRows); %生成 -1 .... 1 共1024个数
[X, Y] = meshgrid(x, x); % X Y
m = size(X,1);
n = peaks(m);
OPD = 0*X +0*Y + 1*n;
a = 1;
b = 1;
I = a + b.*cos(2*pi*OPD) + (rand(mRows)-0.5);
figure(1);
imshow(I, []);
title('spatial domain')
2.通过fft2()函数进行傅里叶变换，通过fftshift()函数将频谱图移到中央：
figure(2);
F = fft2(I);
Fc = fftshift(F);
imshow(abs(Fc),[]);
title('frequency domain')
figure(3);
imshow(abs(log2(Fc)),[]);
title('frequency domain(log2 form)')
画出了原图、频域幅值图和以2为底的幅值的对数图：
image.png
image.png
image.png
三种低通滤波器
低通滤波器实现步骤：
- 给定一幅大小为m*n的图像f(x,y)。选择适当的填充参数P和Q，一般令P = 2m，Q = 2n。
- 对图像f(x, y)填充0，填充后得到图像大小为P*Q的图像fp(x, y)。
- 用(-1)^(x+y)乘以fp(x,y)将其移到变换中心（中心化）。
- 计算fp(x, y)的DFT，得到F(u,v)。
- 生成一个实的，对称的滤波函数H(u, v)，大小为P*Q，中心在(P/2, Q/2)处。然后相乘（矩阵点乘）得到G(u,v) = H(u,v)F(u,v)。
- 对G(u, v)反傅里叶变换，然后取实部，再乘以(-1)^(x+y)进行反中心变换最后得到gp(x,y)。
- 提取gp(x,y)左上角的m*n区域，对提取的部分进行标准化处理，得到最终的结果图像g(x,y)。
#步骤1
clear all;
f = imread('<>') #<>为图片路径
f = double(f);
[m,n]=size(f);
P=2*m;
Q=2*n;
#步骤2
A=zeros(P,Q);
for i=1:m
for j=1:n
A(i,j)=f(i,j);
end
end
#步骤3
for i=1:P
for j=1:Q
A(i,j)=A(i,j)*(-1)^(i+j);
end
end
#步骤4
f1 = fft2(A,P,Q);
步骤5：设计三个低通滤波器的频域函数（D=10/30/60/160/460）
理想低通(步骤5)
H_0 = zeros(P,Q);
H_1 = zeros(P,Q);
H_2 = zeros(P,Q);
H_3 = zeros(P,Q);
H_4 = zeros(P,Q);
for x = (-P/2):1:(P/2)-1
for y = (-Q/2):1:(Q/2)-1
D = sqrt(x^2+y^2);
if(D <= 10) H_0(x+(P/2)+1,y+(Q/2)+1) = 1; end
if(D <= 30) H_1(x+(P/2)+1,y+(Q/2)+1) = 1; end
if(D <= 60) H_2(x+(P/2)+1,y+(Q/2)+1) = 1; end
if(D <= 160) H_3(x+(P/2)+1,y+(Q/2)+1) = 1; end
if(D <= 460) H_4(x+(P/2)+1,y+(Q/2)+1) = 1; end
end
end
butterworth低通滤波器(步骤5)
H_0 = zeros(P,Q);
H_1 = zeros(P,Q);
H_2 = zeros(P,Q);
H_3 = zeros(P,Q);
H_4 = zeros(P,Q);
for x=-(P/2):1:(P/2)-1
for y=(-Q/2):1:(Q/2)-1
D = sqrt(x^2+y^2);
H_0(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/10)^4);
H_1(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/30)^4);
H_2(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/60)^4);
H_3(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/160)^4);
H_4(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D/460)^4);
end
end
高斯低通滤波器(步骤5)
H_0 = zeros(P,Q);
H_1 = zeros(P,Q);
H_2 = zeros(P,Q);
H_3 = zeros(P,Q);
H_4 = zeros(P,Q);
for x=-(P/2):1:(P/2)-1
for y=(-Q/2):1:(Q/2)-1
D = sqrt(x^2+y^2);
D_0 = 10;
H_0(x+(P/2)+1,y+(Q/2)+1) = exp((-D^2)/(2*(D_0^2)));
D_0 = 30;
H_1(x+(P/2)+1,y+(Q/2)+1) = exp((-D^2)/(2*(D_0^2)));
D_0 = 60;
H_2(x+(P/2)+1,y+(Q/2)+1) = exp((-D^2)/(2*(D_0^2)));
D_0 = 160;
H_3(x+(P/2)+1,y+(Q/2)+1) = exp((-D^2)/(2*(D_0^2)));
D_0 = 460;
H_4(x+(P/2)+1,y+(Q/2)+1) = exp((-D^2)/(2*(D_0^2)));
end
end
下面是步骤6、步骤7：
#步骤6
G_0 = H_0 .* f1;
G_1 = H_1 .* f1;
G_2 = H_2 .* f1;
G_3 = H_3.*f1;
G_4 = H_4.*f1;
g_0 = real(ifft2(G_0));
g_0 = g_0(1:1:m,1:1:n);
g_1 = real(ifft2(G_1));
g_1 = g_1(1:1:m,1:1:n);
g_2 = real(ifft2(G_3));
g_2 = g_2(1:1:m,1:1:n);
g_3 = real(ifft2(G_3));
g_3 = g_3(1:1:m,1:1:n);
g_4 = real(ifft2(G_4));
g_4 = g_4(1:1:m,1:1:n);
#步骤7
for x = 1:1:m
for y = 1:1:n
g_0(x,y) = g_0(x,y) * (-1)^(x+y);
g_1(x,y) = g_1(x,y) * (-1)^(x+y);
g_2(x,y) = g_3(x,y) * (-1)^(x+y);
g_3(x,y) = g_3(x,y) * (-1)^(x+y);
g_4(x,y) = g_4(x,y) * (-1)^(x+y);
end
end
最后输出图像
figure();
subplot(2,3,1);
imshow(f,[ ]);
xlabel('1).原图');
subplot(2,3,2);
imshow(g_0,[ ]);
xlabel('2).效果图(D=10)');
subplot(2,3,3);
imshow(g_1,[ ]);
xlabel('3).效果图(D=30)');
subplot(2,3,4);
imshow(g_2,[ ]);
xlabel('4).效果图(D=60)');
subplot(2,3,5);
imshow(g_3,[ ]);
xlabel('5).效果图(D=160)');
subplot(2,3,6);
imshow(g_4,[ ]);
xlabel('6).效果图(D=460)');
