Bohrium
robot
新建

空间站广场

论文
Notebooks
比赛
课程
Apps
我的主页
我的Notebooks
我的论文库
我的足迹

我的工作空间

任务
节点
文件
数据集
镜像
项目
数据库
公开
Digital-Signal-Processing_04_Random Signals and LTI-Systems_The Wiener Filter
English
notebook
Signal Processing
EnglishnotebookSignal Processing
喇叭花
发布于 2023-08-02
推荐镜像 :digital-signal-processing:Digital-Signal-Processing-Lecture
推荐机型 :c2_m4_cpu
Random Signals and LTI-Systems
The Wiener Filter
Signal Model
Transfer Function of the Wiener Filter
Wiener Deconvolution
Interpretation

Random Signals and LTI-Systems

🏃🏻 Getting Started Guide
This document can be executed directly on the Bohrium Notebook. To begin, click the Connect button located at the top of the interface, then select the digital-signal-processing:Digital-Signal-Processing-Lecture Image and choose your desired machine configuration to proceed.

📖 Source
This Notebook is from a famous Signal Processing Lecture. The notebooks constitute the lecture notes to the masters course Digital Signal Processing read by Sascha Spors, Institute of Communications Engineering, Universität Rostock.

This jupyter notebook is part of a collection of notebooks on various topics of Digital Signal Processing. Please direct questions and suggestions to Sascha.Spors@uni-rostock.de.

代码
文本

The Wiener Filter

The Wiener filter, named after Norbert Wiener, aims at estimating an unknown random signal by filtering a noisy distorted observation of the signal. It has a wide range of applications in signal processing. For instance, noise reduction, system identification, deconvolution and signal detection. The Wiener filter is frequently used as basic building block for algorithms that denoise audio signals, like speech, or remove noise from a image.

代码
文本

Signal Model

The following signal model is underlying the Wiener filter

model_wiener_filter.png

The random signal is subject to distortion by the linear time-invariant (LTI) system and additive noise , resulting in the observed signal . The additive noise is assumed to be uncorrelated from . It is furthermore assumed that all random signals are wide-sense stationary (WSS). This distortion model holds for many practical problems, like e.g. the measurement of a physical quantity by a sensor.

The basic concept of the Wiener filter is to apply the LTI system to the observed signal such that the output signal matches as close as possible. In order to quantify the deviation between and , the error signal

is introduced. The error equals zero, if the output signal perfectly matches . In general, this goal cannot be achieved in a strict sense. As alternative one could aim at minimizing the linear average of the error . However, this quantity is not very well suited for optimization since it is signed. Instead, the quadratic average of the error is used in the Wiener filter and other techniques. This measure is known as mean squared error (MSE). It is defined as

the equalities are deduced from the properties of the power spectral density (PSD) and the auto-correlation function (ACF). Above equation is referred to as cost function. We aim at minimizing the cost function, hence minimizing the MSE between the original signal and its estimate . The solution is referred to as minimum mean squared error (MMSE) solution.

代码
文本

Transfer Function of the Wiener Filter

At first, the Wiener filter shall only have access to the observed signal and selected statistical measures. It is assumed that the cross-power spectral density between the observed signal and the original signal , and the PSD of the observed signal are known. This knowledge can either be gained by estimating both from measurements taken at an actual system or by using suitable statistical models.

The optimal filter is found by minimizing the MSE with respect to the transfer function . The solution of this optimization problem goes beyond the scope of this notebook and can be found in the literature, e.g. [Girod et. al]. The transfer function of the Wiener filter is given as

for . No knowledge on the actual distortion process is required besides the assumption of an LTI system and additive noise. Only the PSDs and have to be known in order to estimate from in the MMSE sense. Care has to be taken that the filter is causal and stable in practical applications.

代码
文本

Example - Denoising of an Deterministic Signal

The following example considers the estimation of the original signal from a distorted observation. It is assumed that the original signal is a deterministic signal which is distorted by an LTI system and additive wide-sense ergodic normal distributed zero-mean white noise with PSD . The PSDs and are estimated from and using the Welch technique. The Wiener filter is applied to the observation in order to compute the estimate of . Note that the discrete signals have been illustrated by continuous lines for ease of illustration.

代码
文本
[1]
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sig

N = 2**14 # number of samples
M = 256 # length of Wiener filter
Om0 = 0.1*np.pi # frequency of original signal
N0 = 0.1 # PSD of additive white noise

# generate original signal
s = np.cos(Om0 * np.arange(N))
# generate observed signal
g = 1/20*np.asarray([1, 2, 3, 4, 5, 4, 3, 2, 1])
np.random.seed(1)
n = np.random.normal(size=N, scale=np.sqrt(N0))
x = np.convolve(s, g, mode='same') + n
# estimate (cross) PSDs using Welch technique
f, Pxx = sig.csd(x, x, nperseg=M)
f, Psx = sig.csd(s, x, nperseg=M)
# compute Wiener filter
H = Psx/Pxx
H = H * np.exp(-1j*2*np.pi/len(H)*np.arange(len(H)) *
(len(H)//2)) # shift for causal filter
h = np.fft.irfft(H)
# apply Wiener filter to observation
y = np.convolve(x, h, mode='same')

# plot (cross) PSDs
Om = np.linspace(0, np.pi, num=len(H))

plt.figure(figsize=(10, 4))
plt.subplot(121)
plt.plot(Om, 20*np.log10(np.abs(.5*Pxx)),
label=r'$| \Phi_{xx}(e^{j \Omega}) |$ in dB')
plt.plot(Om, 20*np.log10(np.abs(.5*Psx)),
label=r'$| \Phi_{sx}(e^{j \Omega}) |$ in dB')
plt.title('(Cross) PSDs')
plt.xlabel(r'$\Omega$')
plt.legend()
plt.axis([0, np.pi, -60, 40])
plt.grid()

# plot transfer function of Wiener filter
plt.subplot(122)
plt.plot(Om, 20*np.log10(np.abs(H)))
plt.title('Transfer function of Wiener filter')
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$| H(e^{j \Omega}) |$ in dB')
plt.axis([0, np.pi, -150, 3])
plt.grid()
plt.tight_layout()

# plot signals
idx = np.arange(500, 600)
plt.figure(figsize=(10, 4))
plt.plot(idx, x[idx], label=r'observed signal $x[k]$')
plt.plot(idx, s[idx], label=r'original signal $s[k]$')
plt.plot(idx, y[idx], label=r'estimated signal $y[k]$')
plt.title('Signals')
plt.xlabel(r'$k$')
plt.axis([idx[0], idx[-1], -1.5, 1.5])
plt.legend()
plt.grid()
代码
文本

Let's listen to the observed signal and the output of the Wiener filter

代码
文本
[2]
from scipy.io import wavfile

fs = 8000
wavfile.write('Wiener_observed_signal.wav', fs, np.int16(x*16384))
wavfile.write('Wiener_output.wav', fs, np.int16(y*16384))
代码
文本

Observed signal

./Wiener_observed_signal.wav

Output of Wiener filter

./Wiener_output.wav

代码
文本

Exercise

  • Take a look at the PSDs and the resulting transfer function of the Wiener filter. How does the Wiener filter remove the noise from the observed signal?
  • Change the frequency Om0 of the original signal and the noise power N0 of the additive noise. What changes?

Solution: The cross-PSD has its maximum at the frequency of the original signal , besides this frequency the PSD has approximately the value of the PSD of the additive noise. The resulting Wiener filter has an attenuation of 0 dB at which increases off this frequency. It has the characteristics of a band-pass filter centered at . By filtering the observed signal with the Wiener filter, the major portion of the additive white noise is removed. However, around both the original signal and noise pass the filter. This explains the remaining small deviations between the estimated signal and the original signal . The adaption of the Wiener filter to the original signal and additive noise changes when changing its frequency and PSD.

代码
文本

Wiener Deconvolution

As discussed above, the general formulation of the Wiener filter is based on the knowledge of the PSDs and characterizing the distortion process and the observed signal, respectively. These PSDs can be derived from the PSDs of the original signal and the noise , and the transfer function of the distorting system.

Under the assumption that is uncorrelated from , the PSD can be derived as

and

Introducing these results into the general formulation of the Wiener filter yields

This specialization is also known as Wiener deconvolution filter. The filter can be derived from the PSDs of the original signal and the noise , and the transfer function of the distorting system. This form is especially useful when the PSDs can be modeled reasonably well by analytic probabilty density functions (PSDs). In many cases, the additive noise can be modeled as white noise .

代码
文本

Interpretation

The Wiener deconvolution filter derived above is also useful for an interpretation of the operation of the Wiener filter. It can be rewritten by introducing the frequency dependent signal-to-noise ratio between the original signal and the noise as

Now two special cases are discussed:

  1. If there is no additive noise , the bracketed expression is equal to 1. Hence, the Wiener filter is simply given as the inverse system to the distorting system

  2. If the distorting system is just a pass-through , the Wiener filter is given as

    Two extreme cases are considered

    • for (high ) at a given frequency the transfer function approaches 1
    • for (low ) at a given frequency the transfer function approaches 0

The discussed cases explain the operation of the Wiener filter for specific situations, the general operation is a combination of these.

代码
文本

Example - Denoising of a deterministic signal with the Wiener deconvolution filter

The preceding example of the general Wiener filter will now be reevaluated with the Wiener deconvolution filter.

代码
文本
[3]
N = 2**14 # number of samples
M = 256 # length of Wiener filter
Om0 = 0.1*np.pi # frequency of original signal
N0 = .1 # PSD of additive white noise

# generate original signal
s = np.cos(Om0 * np.arange(N))
# generate observed signal
g = 1/20*np.asarray([1, 2, 3, 4, 5, 4, 3, 2, 1])
np.random.seed(1)
n = np.random.normal(size=N, scale=np.sqrt(N0))
x = np.convolve(s, g, mode='same') + n
# estimate PSD
f, Pss = sig.csd(s, s, nperseg=M)
f, Pnn = sig.csd(n, n, nperseg=M)
# compute Wiener filter
G = np.fft.rfft(g, M)
H = 1/G * (np.abs(G)**2 / (np.abs(G)**2 + N0/Pss))
H = H * np.exp(-1j*2*np.pi/len(H)*np.arange(len(H)) *
(len(H)//2-8)) # shift for causal filter
h = np.fft.irfft(H)
# apply Wiener filter to observation
y = np.convolve(x, h, mode='same')

# plot (cross) PSDs
Om = np.linspace(0, np.pi, num=len(H))

plt.figure(figsize=(10, 4))
plt.subplot(121)
plt.plot(Om, 20*np.log10(np.abs(.5*Pss)),
label=r'$| \Phi_{ss}(e^{j \Omega}) |$ in dB')
plt.plot(Om, 20*np.log10(np.abs(.5*Pnn)),
label=r'$| \Phi_{nn}(e^{j \Omega}) |$ in dB')
plt.title('PSDs')
plt.xlabel(r'$\Omega$')
plt.legend()
plt.axis([0, np.pi, -60, 40])
plt.grid()

# plot transfer function of Wiener filter
plt.subplot(122)
plt.plot(Om, 20*np.log10(np.abs(H)))
plt.title('Transfer function of Wiener filter')
plt.xlabel(r'$\Omega$')
plt.ylabel(r'$| H(e^{j \Omega}) |$ in dB')
plt.axis([0, np.pi, -150, 10])
plt.grid()
plt.tight_layout()

# plot signals
idx = np.arange(500, 600)
plt.figure(figsize=(10, 4))
plt.plot(idx, x[idx], label=r'observed signal $x[k]$')
plt.plot(idx, s[idx], label=r'original signal $s[k]$')
plt.plot(idx, y[idx], label=r'estimated signal $y[k]$')
plt.title('Signals')
plt.xlabel(r'$k$')
plt.axis([idx[0], idx[-1], -1.5, 1.5])
plt.legend()
plt.grid()
代码
文本

Exercise

  • Compare the transfer function of the Wiener deconvolution filter to the transfer function of the general Wiener filter from the preceding example.
  • What is different compared to the general Wiener filter? Why?

Solution: The transfer function of the Wiener deconvolution filter is much smoother than the transfer function of the general Wiener filter from the preceding example. This is due to the knowledge of the transfer function of the disturbing system. The general Wiener filter inherently estimates this from the noisy observations.

代码
文本

Copyright

This notebook is provided as Open Educational Resource. Feel free to use the notebook for your own purposes. The text is licensed under Creative Commons Attribution 4.0, the code of the IPython examples under the MIT license. Please attribute the work as follows: Sascha Spors, Digital Signal Processing - Lecture notes featuring computational examples.

代码
文本
English
notebook
Signal Processing
EnglishnotebookSignal Processing
点个赞吧
推荐阅读
公开
Digital-Signal-Processing 04Random Signals and LTI-Systems_Measurement of Acoustic Impulse Responses
EnglishnotebookSignal Processing
EnglishnotebookSignal Processing
喇叭花
发布于 2023-08-02
公开
Digital-Signal-Processing_04_Random Signals and LTI-Systems_Auto-Correlation Function
EnglishnotebookSignal Processing
EnglishnotebookSignal Processing
喇叭花
发布于 2023-08-02