信号处理，想到就写

从一个实际的小问题出发，回顾一些基础的知识点：

平稳噪声x[n]激励系统h[n]，输出y[n]能量需要和x[n]一致，系统的形式：

A

H(z) = --------------------------

1 - sum(i=1:m, a[i]h[i])

求A。

1. 首先看系统的响应：

平稳随机过程x通过LSI系统h的输出。

输出的自相关是输入自相关、系统和系统反演的卷积。输出功率谱是输入功率谱和系统函数绝对值平方的乘积。

signal: x[n]
system: h[n]
output: y[n]

y[n] = x[n] * h[n]

R{n,r,y} = E{y[n] y[n-r]'}
= E{ (x[n]*h[n]) (x[n-r]'*h[n-r]') }
= E{ sum(i; x[i] h[n-i]) sum(j; x[j]' h[n-r-j]') }
= E{ sum(i,j; x[i] h[n-i] x[j]' h[n-r-j]') }

= E{ sum(i,j; x[i] x[j]' h[n-i] h[n-r-j]') }
since x and h are not correlated and h is deterministic:

= sum(i,j; E{x[i] x[j]'} h[n-i] h[n-r-j]' ) }
let j = i - s
= sum(i,s; E{x[i] x[i-s]'} h[n-i] h[n-r+s-i]' ) }
if x is stationary
(note, if x were not stationary, Rx and h would not even be set apart)
= sum(i,s; R{s,x} h[n-i] h[-(r-s - (n-i))]' ) }
= sum(i,s; R{s,x} h[i] h[-(r-s - i)]' ) }
= sum(s; R{s,x} h[r-s] * h-[r-s]' ) }
= R{r,x} * h[r] * h-[r]'
y is also stationary

in conclusion:
Py(w) = Px(w) | H(w) |^2

如果噪声是白的，就要求 intg(w, |H(w)|^2) = 1

严格的讲H(w)是h[n]的Fourier变换： H(w) = sum(k, h[k] e^(-j*w*k))

而此处H(w)是IIR滤波器。IIR的特点是极点(poles)都在单位圆内部。同时它是一个因果(casual)系统，因此ROC是一个z平面去掉一个圆域，这个圆域应该是包含极点的最小圆。如果ROC包含单位圆那么系统是BIBO的（因此极点都在单位圆内部），这也是此处暗含的要求。另外为什么要对离散时间系统用z变换讨论？实际上z变换只是Laplace变换自变量的一个映射，而用z变换就能包含所有信息（离散变换，即周期冲击的连续信号在Laplace变换域具有纵向的周期性），Laplace变换虚轴左侧都映射到单位圆内部。

根据Parseval定理（实际上就是一个正交变换模长不变性）：

intg(w, |H(w)|^2) = 2 pi sum(k, |h[k]|^2)

因此对于BIBO的IIR情况，忽略冲击响应中远端的较小项，用h[k]的前几项通过上式计算可以估计出这个能量，于是就能求出A进行调整。

另外也有人提出用对偶的FIR进行近似和用FFT转换到频域进行估算的方法。

2. 输入输出信号的互相关/互功率谱

Pxy(w) 为 Rxy[r] = E{x[n] y[n-r]'} = E{ x[n+r]' y[n] }' 的Fourier变换。
y[n] = x[n]*h[n] = sigma(k; x[n-k] h[k])
E{ x[n+r]' y[n] } = E{ x[n+r]' sigma(k; x[n-k] h[k]) } = E{ sigma(k; x[n+r]' x[n-k] h[k]) }
= sigma(k; E{x[n+r]' x[n-k]} h[k]) = sigma(k; Rxx[-r-k] h[k]) = Rxx[-r] * h[-r]
E{ x[n+r]' y[n] }' = Rxx[-r]' * h[-r]' = Rxx[r] * h[r]
因此Pxy(w) = Pxx(w) H(w)