反问题: Wiener滤波器进行反卷积
编辑目的
在信号处理中，真实信号x 经过某个 滤波器H 之后，我们获得了观测信号y，但通常y 会受到 噪声n 的影响，也就是信噪混合（要么脉冲相近产生混淆，要么被淹没在强噪声里），使得我们无法直接反卷积获得高质量的x。
因此维纳滤波器作为最优线性滤波器，它通过频域处理来逆转卷积的影响并去除噪音，主要基于信号和噪声的功率谱密度（DSP）来调整增益，使得输出信号尽可能接近原始信号。
背景
假设系统模型为
其中 h(t) 表示系统冲激响应（或点扩散函数）， x(t) 是原始信号， e(t) 是加性噪声。
在频域中可写为
我们希望设计一个滤波器 W(\omega) ，使得对 Y(\omega) 进行滤波后，能够得到对原始信号 X(\omega) 的最佳还原 \hat{X}(\omega) ，数学形式表达为:
我们最简单的想法是直接使用 H^{-1}(\omega)=1/H(\omega) ，也就是直接除以原滤波器，但是这样有一个问题就是，如果 H(\omega) 的某些频段幅度很小或接近零，那么在一除就会导致噪声 E(\omega) 大幅度放大，并且如果 H(\omega) 不可逆，那更没办法求了。因此，在有噪声的情况下，单纯的逆滤波不行，这就需要Wiener滤波器。
Wiener 滤波器的目标是最小化均方误差（MSE），即最小化 \mathbb{E}\{|X(\omega)-\hat{X}(\omega)|^2\} 。推导过程基于以下假设：
x(t) 和 e(t) 是平稳随机过程
x(t) 与 e(t) 相互统计独立
系统 h(t) 是线性、时不变（LTI系统）
在频域中，Wiener滤波器 W(\omega) 的最优解如下（具体推导过程见 随机信号理论）
S_x(\omega) ：信号 x(t) 的功率谱，可以根据对信号特性的先验知识来进行估计
S_e(\omega) ：噪声 e(t) 的功率谱，可以根据系统在无信号输入的时候的输出来估计噪声分布
H(\omega) ：系统冲激响应 h(t) 的频域表示
H^*(\omega) 是 H(\omega) 的复共轭， |H(\omega)|^2 = H(\omega)H^*(\omega)
当信噪比较高（ S_X(\nu) \gg S_B(\nu) ）时，分母中\frac{S_B(\nu)}{S_X(\nu)} 这一项接近 0，此时W(\nu) \approx \frac{1}{H(\nu)}，接近逆滤波，恢复原始信号。
当信噪比较低（ S_X(\nu) \ll S_B(\nu) ）时，噪声占主导， W(\nu) 变小，滤波器会抑制该频率分量，避免噪声放大。
在实际应用中，我们需要实现知道 S_x(\omega) 和 S_e(\omega) ，如果不知道这些统计量，那么维纳滤波器就玩不了，这是其一个不好的地方。
功率谱密度的选择
短平稳序列的情况
若观测信号y 与输入信号x 具有很高的相关性，可以使用 自回归（AR）模型 近似输入信号的行为。AR 模型是基于当前值由前面的值加上一个随机误差来决定的，其数学形式如下:
复杂信号的情况
如果信号是较为复杂的 图像、重叠频率、密集功率谱分布，功率谱密度（DSP）可以用来反映这些纹理的频率特征，因此研究DSP就可以研究图像纹理细节。
综上所述，对于一个相对平滑、连续且正相关(当前值与之前的值有很强的正相关性)的信号，我们可以依赖于一个自回归系数接近1的信号的功率谱密度。自回归系数接近1意味着 当前值与之前值关系非常强，信号变化非常平缓。
于此同时，我们希望对所有频率成分一视同仁，不额外惩罚某些频率，因此选定单位 DSP (即每个频率的强度都是相同的，其功率谱密度是平坦的)，这就是白噪音假设的来源
利用先验信息对功率谱密度的估计思路
我们之前提到用先验信息来估计信号的功率谱S_X(\nu)，具体有以下几种思路：
1. 直接利用历史数据（已有大量历史数据情况）
可通过周期图来计算其功率谱密度:
其中: X(\nu) 是 x(t) 的傅里叶变换
2. 自回归（AR）模型
信号如果可以用一个自回归过程（Auto-Regressive） 来描述：
那么 AR 模型的功率谱密度可以表示为：
适用于平稳信号
3. 马尔可夫模型
如果信号是非平稳的，可以使用 马尔可夫模型（Hidden Markov Model, HMM） 估计不同状态下的功率谱：
每个隐藏状态 s_i 对应一个功率谱 S_X^{(i)}(\nu)
最终功率谱是所有状态的加权和：
4. 最大后验估计（MAP）
最大后验估计（MAP）是一种基于贝叶斯思想的方法，用于在结合观测数据和先验知识的基础上估计参数。在估计功率谱密度 S_X(\nu) 的问题中，我们可以将 S_X(\nu) 看作待估计的参数，然后利用 MAP 方法求解其后验分布的最大值。
首先，我们需要构造似然函数，即在给定谱密度 S_X(\nu) 的条件下观测数据 X 出现的概率。对于很多实际问题（例如高斯过程模型），假设观测数据服从高斯分布，其协方差矩阵与 S_X(\nu) 有关，那么似然函数 p(X|S_X(\nu)) 就可以根据高斯分布给出。
其次，我们引入先验分布 p(S_X(\nu)) ，以表达对谱密度的先验认识。例如，我们可能期望谱密度具有平滑性或者满足某些物理约束，通过先验可以将这些信息融入估计中。
根据贝叶斯公式，后验分布为
MAP 估计即为使后验概率最大的 S_X(\nu) ，可以写作
通常，为了简化计算，我们取对数将乘积转化为求和，得到
通过这种方式，MAP 估计不仅依赖于观测数据（通过似然函数），还融入了先验知识，使得在数据有限或噪声较大的情况下能够得到更稳健的谱密度估计。实际求解时，往往需要利用数值优化方法（例如梯度下降或 EM 算法）来求解上述最大化问题。
综上所述
自回归模型和线性回归模型的区别
由于这里有一篇 线性回归模型 的内容，为了避免混淆简单做一下区分，首先自回归模型主要用于时间序列数据，其核心思想是利用序列过去的值来预测当前或未来的值。其形式为:
而线性回归模型则主要用于描述因变量与一个或多个自变量之间的线性关系，模型通常写成
实现
一维情况
加载数据文件 Observations.mat，并查看滤波器在时域和傅里叶域中的特性
clear all
clc
close all
%% 加载数据并可视化滤波器
% 加载数据文件 Observations.mat
load('Observations.mat'); % 这将加载变量 RI(滤波器) 和 Donnees(观测信号)
% 在时域中可视化滤波器
figure;
subplot(2,1,1);
plot(RI);
title('时域中的滤波器');
xlabel('时间样本点');
ylabel('幅度');
% 计算滤波器的傅里叶变换
H = MyFFT(RI); % 使用 MyFFT 得到中心化频域的 H
N = length(RI);
% 频率轴设置
freq = (-N/2:N/2-1)/N; % 从 -0.5 到 0.5 的归一化频率轴
% 在频域中可视化滤波器的幅度
subplot(2,1,2);
plot(freq, abs(H));
title('频域中的滤波器(幅度)');
xlabel('归一化频率');
ylabel('幅度');
在时间域，滤波器的特性通过其冲激响应展现出来。冲激响应是指滤波器对单位冲激信号（Dirac delta 函数）的响应，完全描述了该滤波器的行为。
在频率域，滤波器的特性通过其频率响应的幅值图来表示。幅值图表明了滤波器如何处理不同的频率分量，是高通还是低通。这里明显低通。
根据时频分析的基本原理，脉冲越窄，频谱越宽，符合我们的结果图
然后绘制观测信号
%% 可视化观测信号
figure;
plot(Donnees);
title('观测信号');
xlabel('时间样本点');
ylabel('幅度');
信号经过滤波器处理时，卷积（convolution）使信号变得模糊，也就是我们观测到的信号是不纯的（不完美尖锐），再加上噪声的干扰，非常不好，使得单脉冲变得不易区分，甚至有可能出现邻近脉冲的重叠问题，这里的观测信号出现了尿分叉的现象，所以要用反卷积思想，这里我们用维纳滤波器来做
维纳滤波器
我们先固定一个 λ 的值来 绘制维纳滤波器的增益，然后将该增益曲线与滤波器 H 的特性放在一起看一下。 λ 是噪声功率谱密度与信号功率谱密度之间的比率:
这里实验中做了简化，不再估计计算 S_X (ν) ，而是直接假设 λ =0.9
%% 计算并可视化维纳滤波增益
% 固定 lambda 的值 (lambda = S_B / S_X)
lambda = 0.9;
% 计算维纳滤波器增益
W = conj(H) ./ (abs(H).^2 + lambda);
% 绘制维纳滤波增益并与滤波器 H 对比
figure;
plot(freq, abs(W), 'LineWidth', 2);
hold on;
plot(freq, abs(H), '--', 'LineWidth', 2);
title(['维纳滤波增益与滤波器 H (\lambda = ', num2str(lambda), ')']);
xlabel('归一化频率');
ylabel('幅度');
legend('维纳滤波增益', '滤波器 H');
grid on;
这时在频域（横轴是归一化频率）下，将维纳滤波器增益（Wiener Filter Gain）与系统原滤波器 H 的增益进行对比，直观看到维纳滤波器如何对噪声部分进行调节和削弱，这里结果显示维纳滤波器和原滤波器增益曲线基本重叠，因为我们把 λ 设为0.9 ，代表噪声影响大，信号影响小，信噪比SNR 约为 1.1，不大不小，所以维纳滤波器基本没做明显的增益衰减来抑制噪声，所以我们把 λ 设为10，这样信噪比SNR就是0.1，噪声强度大于信号强度，维纳滤波器倾向于抑制噪声，因此会大幅度衰减维纳滤波器增益，如下图所示。
我们再取不同 λ 来看其他情况
%% 不同 lambda 对维纳滤波器的影响
% 定义一组 lambda 值
lambda_values = [0.01, 0.1, 0.9];
figure;
hold on;
for i = 1:length(lambda_values)
lambda_val = lambda_values(i);
W_temp = conj(H) ./ (abs(H).^2 + lambda_val);
plot(freq, abs(W_temp), 'LineWidth', 1.5);
end
title('不同 \lambda 值下的维纳滤波增益');
xlabel('归一化频率');
ylabel('幅度');
legend('\lambda=0.01', '\lambda=0.1', '\lambda=0.9');
grid on;
增益是指滤波器对信号的增强或衰减程度，λ 的本质其实就是 信噪比的倒数
综上所述我们得到结论:
如果 \lambda \to 0 ，λ 小 → 噪声较弱 → 信噪比大，维纳滤波器趋于完全恢复信号，即 理想逆滤波形态，在主要频率段不做太多衰减（甚至有所放大），因为此时噪声影响小，可以放心补偿。
但是如果 λ 大 → 噪声较强 → 信噪比小，Wiener 滤波器在低信噪比的情形下会大幅衰减增益，以防止过度放大噪声。
%% 使用维纳滤波器对信号进行重建
lambda = 0.1; % 为重建选择合适的 lambda
W = conj(H) ./ (abs(H).^2 + lambda);
% 计算观测信号的傅里叶变换
Y = MyFFT(Donnees); % 频域观测信号(已居中)
% 在频域中应用维纳滤波器
X_hat = Y .* W; % 不要加abs，会丢相位信息
% 在时域中重建信号
x_hat = MyIFFT(X_hat);
% 对重建的时域信号进行循环移位以校正相位
x_hat = circshift(x_hat, N/2);
% 可视化观测与重建信号
figure;
subplot(2,1,1);
plot(Donnees);
title('观测信号');
xlabel('时间样本点');
ylabel('幅度');
subplot(2,1,2);
plot(real(x_hat));
title('使用维纳滤波器重建的信号(相位已校正)');
xlabel('时间样本点');
ylabel('幅度');
根据结果可见，维纳滤波器很好的从被模糊或被卷积后的信号中还原原始信号形状，让脉冲特征更加易于辨别 ，但是由于我们的原始信号不可知，所以没办法直接比较真实的原始信号来评估性能，非常遗憾
二维情况
clear all
clc
close all
%% 加载数据
% 假设 DataOne.mat 中包含变量 Data、IR、TrueImage
Data1 = load('DataOne.mat');
Y_obs = Data1.Data; % 观测图像（已退化）
RI = Data1.IR; % 2D滤波器内核
X_true = Data1.TrueImage; % 与 Y_obs 尺寸相同的真实清晰图像
%% 显示观测图像
figure;
imagesc(Y_obs);
colormap gray;
axis image off;
title('Observed Degraded Image');
drawnow;
%% 计算滤波器的频域表示
[rows, cols] = size(Y_obs);
Long = max(rows, cols); % 根据需要设置扩展后的大小，可与图像大小匹配
% 计算滤波器内核 RI 的傅里叶变换表示 H
H = MyFFT2RI(RI, Long);
% 频率轴（归一化频率，仅用于绘图，可选）
freq_x = (-cols/2:cols/2-1)/cols;
freq_y = (-rows/2:rows/2-1)/rows;
%% 定义维纳滤波参数
% 这里设定 lambda 为噪声功率与信号功率的比值
lambda = 0.05;
% 构建二维维纳滤波器增益
W = conj(H) ./ (abs(H).^2 + lambda);
%% 对观测图像进行频域变换
Y = MyFFT2(Y_obs);
%% 应用维纳滤波器恢复图像
X_hat_freq = Y .* W;
x_hat = MyIFFT2(X_hat_freq);
%% 计算 MSE 和 PSNR
if ~exist('X_true','var')
error('真实图像 X_true 缺失，请后续补充提供。');
end
mse_value = mean((double(X_true(:)) - double(x_hat(:))).^2);
max_val = max(X_true(:));
psnr_value = 10*log10((max_val^2)/mse_value);
fprintf('MSE: %.4f\n', mse_value);
fprintf('PSNR: %.4f dB\n', psnr_value);
%% 可视化恢复结果
figure;
subplot(1,3,1);
imagesc(Y_obs);
colormap gray;
axis image off;
title('Observed (Degraded) Image');
subplot(1,3,2);
imagesc(real(x_hat));
colormap gray;
axis image off;
title(sprintf('Restored Image using Wiener Filter\nMSE=%.4f, PSNR=%.2f dB', mse_value, psnr_value));
subplot(1,3,3);
imagesc(X_true);
colormap gray;
axis image off;
title('True Image');
drawnow;
%% 分图绘制像素值随位置变化的对比图
% 取图像中间行
row_idx = floor(rows/2);
true_line = double(X_true(row_idx, :)); % 真实图像中该行像素强度
obs_line = double(Y_obs(row_idx, :)); % 观测图像中该行像素强度
est_line = double(x_hat(row_idx, :)); % Wiener 恢复图像中该行像素强度
% 图1：真实图像与观测图像对比
figure;
plot(true_line, 'k-', 'LineWidth', 2); % 黑色实线：真实图像
hold on;
plot(obs_line, 'r-.', 'LineWidth', 2); % 红色点划线：观测图像
hold off;
title('Row Intensity Profile: True vs. Observed');
xlabel('Pixel Index');
ylabel('Intensity');
legend('True','Observed','Location','best');
grid on;
% 图2：真实图像与 Wiener 恢复图像对比
figure;
plot(true_line, 'k-', 'LineWidth', 2); % 黑色实线：真实图像
hold on;
plot(est_line, 'b--', 'LineWidth', 2); % 蓝色虚线：Wiener 恢复图像
hold off;
title('Row Intensity Profile: True vs. Wiener');
xlabel('Pixel Index');
ylabel('Intensity');
legend('True','Wiener','Location','best');
grid on;
第一个图对比 观测图像 与 真实图像 的像素强度，第二个图对比 Wiener 恢复图像 与 真实图像的像素强度。
部分自定义函数如下
%% MyIFFT
function Temporel = MyIFFT(Frequentiel)
% MyIFFT
% 对一维频域数据执行逆FFT操作，包括 fftshift 和 sqrt(N) 归一化，恢复时域信号。
%
% 输入参数：
% Frequentiel : 一维频域数据(已 fftshift)
%
% 输出参数：
% Temporel : 逆FFT后的时域信号数组
Temporel = ifft( fftshift(Frequentiel) ) * sqrt(length(Frequentiel));
end
%% MyFFTRI
function Frequentiel = MyFFTRI(Temporel)
% MyFFTRI
% 对一维时域数据进行：先 fftshift，再 fft，然后再 fftshift，
% 得到频域中心化的频谱数据(不归一化)。
%
% 输入参数：
% Temporel : 一维时域数据
%
% 输出参数：
% Frequentiel : 中心化后的频域数据(不做归一化)
Frequentiel = fftshift( fft(fftshift(Temporel)) );
end
%% MyFFT
function Frequentiel = MyFFT(Temporel)
% MyFFT
% 对一维时域数据进行 fft，然后 fftshift，并通过 sqrt(N) 归一化结果，
% 生成零频居中的一维频域数据。
%
% 输入参数：
% Temporel : 一维时域数据
%
% 输出参数：
% Frequentiel : 经 fftshift 和 sqrt(N) 归一化的一维频域数据
Frequentiel = fftshift( fft(Temporel) ) / sqrt(length(Temporel));
end
%% MyFFT2
function Frequentiel = MyFFT2(Spatial)
% MyFFT2
% 对二维空间域数据 Spatial 进行 2D FFT，并对结果进行 fftshift 和归一化处理。
%
% 输入参数：
% Spatial : 二维时域(或空间域)数据矩阵
%
% 输出参数：
% Frequentiel : 已经经过 fftshift 和归一化的二维频域数据矩阵
Frequentiel = fftshift( fft2(Spatial) ) / length(Spatial);
end
%% MyFFT2RI
function Frequentiel = MyFFT2RI(Spatial, Long)
% MyFFT2RI
% 对输入的二维数据 Spatial 进行补零扩展至尺寸 Long x Long，以便在频域处理时
% 保持良好的居中对齐。补零后对数据进行 fftshift、2D FFT 再 fftshift。
%
% 输入参数：
% Spatial : 原始二维数据矩阵
% Long : 补零后的目标矩阵尺寸(长宽相等)
%
% 输出参数：
% Frequentiel : 已补零扩展并居中对齐的二维频域数据矩阵
Taille = length(Spatial);
Ou = 1 + Long/2 - (Taille-1)/2 : 1 + Long/2 + (Taille-1)/2;
SpatialComplet = zeros(Long, Long);
SpatialComplet(Ou, Ou) = Spatial;
Frequentiel = fftshift( fft2( fftshift(SpatialComplet) ) );
end
%% MyIFFT2
function Spatial = MyIFFT2(Frequentiel)
% MyIFFT2
% 对二维频域数据执行逆FFT操作并进行 fftshift 和归一化，以还原回空间/时域信号。
%
% 输入参数：
% Frequentiel : 二维频域数据矩阵(已 fftshift)
%
% 输出参数：
% Spatial : 逆变换并归一化后的二维时域(空间域)数据矩阵
Spatial = ifft2( fftshift(Frequentiel) ) * length(Frequentiel);
end
总结
总而言之我们想做的就是一件事，反卷积，直接采用逆滤波器 1/H(\nu) 往往会在 H(\nu) 较小时放大噪声，为此我们使用维纳滤波器，它在频域中加入噪声和信号的先验知识（功率谱密度之比）来实现均方误差最小化，根据信噪比大小不同来着重恢复信号或者抑制噪音，在实验中我们通过调节 λ 的取值（简化）来调节信噪比，展示了维纳滤波器在两种情况（大 小）下的增益变化。