用Matlab画一下Hanning和Hamming滤波函数的频域响应图
时间: 2024-05-10 15:20:27 浏览: 22
以下是Matlab代码和图像:
```
% 生成Hanning滤波函数
N = 51; % 滤波器长度
hanning = hann(N);
hann_resp = fftshift(fft(hanning, 1024)); % 计算频域响应
freq = linspace(-pi, pi, length(hann_resp));
figure;
subplot(2,1,1);
plot(freq, abs(hann_resp));
title('Hanning滤波函数频域响应');
xlabel('频率(rad)');
ylabel('幅值');
% 生成Hamming滤波函数
hamming = hamming(N);
hamm_resp = fftshift(fft(hamming, 1024)); % 计算频域响应
subplot(2,1,2);
plot(freq, abs(hamm_resp));
title('Hamming滤波函数频域响应');
xlabel('频率(rad)');
ylabel('幅值');
```
![Hanning和Hamming滤波函数的频域响应图](https://img-blog.csdn.net/20180528151121664?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3F1aW5fYmFzZTY0/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70/q/70)
相关问题
使用RL,SL,Hamming,Hanning,Cosine滤波函数的滤波反投影重建图像的matlab代码
由于没有提供具体的数据集和重建算法,以下代码仅提供基本的滤波函数实现:
1. RL滤波函数:
```matlab
function h = RL_filter(N, d0, a)
% N: filter size
% d0: cutoff frequency
% a: parameter
h = zeros(N,1);
for i = 1:N
u = (i - N/2 - 1)/(N/2);
if abs(u) <= d0/a
h(i) = 1;
elseif abs(u) > d0/a && abs(u) <= d0
h(i) = 0.5 * (1 + cos(pi*a*(abs(u)-d0/a)/(d0-d0/a)));
else
h(i) = 0;
end
end
end
```
2. SL滤波函数:
```matlab
function h = SL_filter(N, d0, a)
% N: filter size
% d0: cutoff frequency
% a: parameter
h = zeros(N,1);
for i = 1:N
u = (i - N/2 - 1)/(N/2);
if abs(u) <= d0/a
h(i) = 1;
elseif abs(u) > d0/a && abs(u) <= d0
h(i) = (d0 - abs(u))/(d0 - d0/a);
else
h(i) = 0;
end
end
end
```
3. Hamming窗函数:
```matlab
function w = Hamming_window(N, alpha)
% N: window size
% alpha: parameter
w = zeros(N,1);
for i = 1:N
w(i) = alpha - (1-alpha)*cos((2*pi*(i-1))/(N-1));
end
end
```
4. Hanning窗函数:
```matlab
function w = Hanning_window(N)
% N: window size
w = zeros(N,1);
for i = 1:N
w(i) = 0.5 - 0.5*cos((2*pi*(i-1))/(N-1));
end
end
```
5. Cosine窗函数:
```matlab
function w = Cosine_window(N)
% N: window size
w = zeros(N,1);
for i = 1:N
w(i) = sin((pi*(i-0.5))/N);
end
end
```
6. 滤波反投影重建代码:
```matlab
function img = FBPR_reconstruction(proj, theta, N, d, filter, window)
% proj: projection data (N x n_theta)
% theta: projection angles (n_theta x 1)
% N: image size (N x N)
% d: pixel size
% filter: filter function (handle)
% window: window function (N x N)
n_theta = length(theta);
img = zeros(N,N);
for i = 1:n_theta
% filtering
proj(:,i) = filter_projection(proj(:,i), filter, N, d);
% backprojection
img = img + backprojection(proj(:,i), theta(i), N, d);
end
% windowing
img = img .* window;
% normalization
img = img / (pi*n_theta/2);
end
function proj_f = filter_projection(proj, filter, N, d)
% proj: projection data
% filter: filter function (handle)
% N: image size (N x N)
% d: pixel size
% FFT
proj_f = fft(proj);
% frequency axis
freq = (-N/2:N/2-1)/(N*d);
% filtering
filter_f = filter(N, freq);
proj_f = proj_f .* filter_f.';
% inverse FFT
proj_f = ifft(proj_f);
end
function img_bp = backprojection(proj, theta, N, d)
% proj: projection data
% theta: projection angle
% N: image size (N x N)
% d: pixel size
img_bp = zeros(N,N);
% coordinates of image center
xc = (N+1)/2;
yc = (N+1)/2;
% coordinates of image pixels
[x,y] = meshgrid(1:N,1:N);
x = (x - xc)*d;
y = (y - yc)*d;
% backprojection
for i = 1:length(proj)
t = theta(i);
p = proj(i);
img_bp = img_bp + p * interp1([cos(t-pi/2),cos(t+pi/2)],[x(:),y(:)],cos(t)*x(:)+sin(t)*y(:),'linear',0);
end
end
```
为什么RL、SL、Hamming和Hanning滤波函数都存在一定的频谱泄漏问题
RL、SL、Hamming和Hanning滤波函数都是窗函数,它们的作用是在时域上加窗,使得频域上的滤波响应更加平滑。然而,这些窗函数都存在一定的频谱泄漏问题,主要原因是窗函数的频率响应不为零。
在频域上,窗函数的频率响应类似于低通滤波器,它们会对低频信号进行增益,对高频信号进行衰减。但是,由于窗函数的频率响应不为零,它们会对一定范围内的中心频率进行增益,这就导致了频谱泄漏问题。
具体来说,当使用窗函数对信号进行滤波时,窗函数会在时域上对信号进行加窗,然后再在频域上对加窗后的信号进行滤波。由于窗函数的频率响应不为零,当滤波信号的频率和窗函数的中心频率相差较小时,窗函数会将这些频率成分进行增益,导致频谱泄漏。这种泄漏会在频率相差越大时越小,所以通常在滤波时会将窗函数的长度设置得越长,以减小频谱泄漏的影响。