Wiener Filter in MATLAB

 

MATLAB Code 

% Wiener Filter Based on Wiener-Hopf Equation
% This script demonstrates how to apply the Wiener filter to recover
% a reference signal from a noisy signal using the Wiener-Hopf equation.
% The filter minimizes the mean squared error between the noisy signal and the reference signal.

clear; close all; clc;

% Signal Parameters
fs = 4000; % Sampling frequency (Hz)
T = 1; % Total recording time (seconds)
L = T * fs; % Signal length (samples)
tt = (0:L-1) / fs; % Time vector
ff = (0:L-1) * fs / L; % Frequency vector

% Generate Reference Signal (a sinusoid)
y = sin(2 * pi * 120 * tt); % Reference sinusoidal signal
y = y(:); % Ensure column vector

% Create Noisy Signal by Adding Gaussian Noise
x = 0.50 * randn(L, 1) + y; % Noisy signal
x = x(:); % Ensure column vector

% Define Filter Order (Number of Coefficients)
N = 200;

% Apply Wiener Filter using custom function
[xest, b, MSE] = wienerFilt(x, y, N);

% Plot Results
figure;
subplot(411);
plot(tt, x, 'k'), hold on, plot(tt, y, 'r');
title('Wiener Filtering Example');
legend('Noisy Signal', 'Reference Signal');

subplot(412);
plot(tt(N+1:end), xest, 'k');
legend('Estimated Signal');

subplot(413);
plot(tt(N+1:end), (x(N+1:end) - xest), 'k');
legend('Residue Signal');

subplot(414);
stem(0:N-1, b, 'k'); % Stem plot for Wiener filter coefficients
title('Wiener Filter Coefficients');
xlabel('Coefficient Index');
ylabel('Magnitude');
grid on;

% Function to Compute Wiener Filter Based on Wiener-Hopf Equations
function [xest, B, MSE] = wienerFilt(x, y, N)
% Inputs:
% x - Noisy signal
% y - Reference signal
% N - Filter order (number of filter coefficients)

% Outputs:
% xest - Estimated signal after filtering
% B - Wiener filter coefficients
% MSE - Mean squared error (performance metric)

% Perform FFT on the first N samples of x (noisy signal) and y (reference)
X = 1/N .* fft(x(1:N)); % FFT of noisy signal
Y = 1/N .* fft(y(1:N)); % FFT of reference signal

% Reshape to column vectors
X = X(:);
Y = Y(:);

% Compute the Autocorrelation Function (Rxx) and Cross-correlation Function (Rxy)
Rxx = N * real(ifft(X .* conj(X))); % Autocorrelation
Rxy = N * real(ifft(X .* conj(Y))); % Cross-correlation

% Convert Rxx into a Toeplitz matrix
Rxx = toeplitz(Rxx);

% Transpose Rxy for the equation
Rxy = Rxy';

% Solve the Wiener-Hopf equation to get filter coefficients: B = inv(Rxx) * Rxy
B = Rxy / Rxx;
B = B(:); % Flatten the filter coefficients vector

% Apply Wiener filter to estimate the signal from the noisy signal
xest = fftfilt(B, x);
xest = xest(N+1:end); % Remove the first N samples to avoid distortion

% Calculate the Mean Squared Error (MSE)
MSE = mean((y(N+1:end) - xest).^2); % Calculate MSE
end

Output

 

 

Further Reading

1.  Wiener Filter (Theory)