精华内容
下载资源
问答
  • 维纳滤波

    2020-07-22 10:46:07
    举例:基于频域维纳滤波的语音增强 还是利用上面的模型: y(n)=x(n)+w(n) 这里y(n)是麦克风接收的带噪语音,x(n)是干净语音信号,w(n)为白噪声。显然相关函数我们无法得知。 利用一种近似的处理思路:利用前面几个分...

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    举例:基于频域维纳滤波的语音增强

    还是利用上面的模型:

    y(n)=x(n)+w(n)

    这里y(n)是麦克风接收的带噪语音,x(n)是干净语音信号,w(n)为白噪声。显然相关函数我们无法得知。

    利用一种近似的处理思路:利用前面几个分帧不带语音,估计噪声,从而得到噪声的功率谱近似,利用带噪语音功率减去噪声功率,得到

    H(ωk) = Pxy(ωk) / Pyy(ωk) = Pxx(ωk) / (Pxx(ωk)+Pnn(ωk)) = (Pyy(ωk)−Pnn(ωk))/ Pyy(ωk)

    利用估计出的维纳滤波器,即可实现信号的频域滤波。这里只是想到的一个实际例子,至于参数估计、迭代方式则是百花齐放了。

    附上主要代码:

    nw = 512; ni = 64;
    NIS = 100;
    Y = fft(enframe(y,hamming(nw),ni)');
    Yab = abs(Y);
    Nest = mean(Yab(:,1:NIS),2);
    Yest = zeros(size(Y));
    for i = 1:size(Y,2)
        Yest(:,i) = Yab(:,i).*((Yab(:,i)-Nest)./Yab(:,i));
    end
    Ye = Yest.*exp(1j*angle(Y));
    result = zeros(1,length(y));%estimation
    for i =1:size(Y,2);
        pos = ((i-1)*ni+1):((i-1)*ni+nw);
       result(pos) = result(pos)+real(ifft(Ye(:,i))).';
    end
    result = result/max(abs(result));
    

    上面思路处理结果:
    在这里插入图片描述

    可以看出维纳降噪多少还是有些效果的,H(ωk) = Pxy(ωk) / Pyy(ωk) = Pxx(ωk) / (Pxx(ωk)+Pnn(ωk))= 11 + 1/SNR可以看出SNR越小,维纳滤波器衰减越大。

    展开全文
  • webrtc中的噪声抑制之一:频域维纳滤波 前言 在开源的噪声抑制算法中,webrtc ns是很有名的,社区里也有很多分享的文章,但要么深要么浅,还有一些误导读者的,所以趁着移植项目的机会,从盲人摸象到庖丁解牛的学习一番...

    webrtc中的噪声抑制之一:频域维纳滤波

    前言

    在开源的噪声抑制算法中,webrtc ns是很有名的,社区里也有很多分享的文章,但要么深要么浅,还有一些误导读者的,所以趁着移植项目的机会,从盲人摸象到庖丁解牛的学习一番这里面的算法原理和工程实现。
    WebRtc Ns模块采用的是频域维纳滤波的方法,结合VAD检测得到前验信噪比和后验信噪比,算出频域维纳滤波器的系数,在频域实现了噪声的滤除。该模块有定点和浮点两部分代码,其中定点代码相对比较复杂,而浮点代码主要分WebRtcNs_AnalyzeCore和WebRtcNs_ProcessCore两部分模块,先分析当前帧,再对当前帧进行滤波,结构清晰,方便阅读和学习,所以本文研究浮点代码的实现。这里吐槽一下有的博主分享的代码,居然只调用WebRtcNs_ProcessCore来实现,不知道这些事例代码能将噪声滤除到什么水平。

    频域维纳滤波器

    使用维纳滤波进行语音降噪的过程,其实是把降噪过程视为一个线性时不变系统,当带噪语音通过这个系统时,在均方误差最小化准则下,使得系统的输出与期望的纯净语音信号最接近的过程。
    实际当中我们能观察到的信号y(n)y(n)是含有噪声的。 假设带噪语音信号可以表述为
    y(n)=s(n)+n(n) y(n) = s(n)+n(n)
    s(n)s(n)为语音信号,d(n)d(n)为加性噪声。维纳滤波方法就是涉及一个数字滤波器h(n)h(n),当输入y(n)y(n)经过滤波器后,滤波器的输出可以最大程度的逼近s(n)s(n)
    s^(n)=y(n)h(n)=k=0K1y(nk)h(k) \hat{s}(n)=y(n)*h(n)=\sum_{k=0}^{K-1} y(n-k)h(k)
    上式进行DFT得到
    S^(ω)=H(ω)Y(ω) \hat{S}(\omega)=H(\omega)Y(\omega)
    在任意频率ωk\omega_k,我们可以计算出误差估计
    E(ωk)=S(ωk)S^(ωk)=S(ωk)H(ωk)Y(ωk) E(\omega_k)=S(\omega_k)-\hat{S}(\omega_k)=S(\omega_k)-H(\omega_k)Y(\omega_k)
    由此定义频域均方误差的代价函数
    J2=E[E(ωk)2]=E{[S(ωk)H(ωk)Y(ωk)][S(ωk)H(ωk)Y(ωk)]}=E[S(ωk)2]H(ωk)E[S(ωk)Y(ωk)]H(ωk)E[S(ωk)Y(ωk)]+H(ωk)2E[Y(ωk)2]=E[S(ωk)2]H(ωk)Pys(ωk)H(ωk)Psy(ωk)+H(ωk)2Pyy(ωk) \begin{aligned} J_2&=E[|E(\omega_k)|^2]\\ &=E\lbrace[S(\omega_k)-H(\omega_k)Y(\omega_k)]^*[S(\omega_k)-H(\omega_k)Y(\omega_k)]\rbrace\\ &=E[|S(\omega_k)|^2]-H(\omega_k)E[S^*(\omega_k)Y(\omega_k)]-H^*(\omega_k)E[S(\omega_k)Y^*(\omega_k)]+|H(\omega_k)|^2E[|Y(\omega_k)|^2]\\ &=E[|S(\omega_k)|^2]-H(\omega_k)P_{ys}(\omega_k)-H^*(\omega_k)P_{sy}(\omega_k)+|H(\omega_k)|^2P_{yy}(\omega_k) \end{aligned}
    公式的最后一行定义了Pyy(ωk)=E[Y(ωk)2]P_{yy}(\omega_k)=E[|Y(\omega_k)|^2],叫做信号的功率谱。Pys(ωk)=E[S(ωk)Y(ωk)]P_{ys}(\omega_k)=E[S^*(\omega_k)Y(\omega_k)] 定义为输入信号y(n)y(n)和期望信号s(n)s(n)的互功率谱,Psy(ωk)=E[S(ωk)Y(ωk)]P_{sy}(\omega_k)=E[S(\omega_k)Y(\omega_k)^*]是互功率谱的另外一种形式,Pys(ωk)=Psy(ωk)P_{ys}(\omega_k)=P_{sy}^*(\omega_k)。对J2J_2求复导数,并令其等于0,得到如下结果:(共轭求偏导这块有点没理解)
    J2H(ωk)=H(ωk)Pyy(ωk)Pys(ωk)=[H(ωk)Pyy(ωk)Psy(ωk)]=0 \frac{\partial J_2}{\partial H(\omega_k)}=H^*(\omega_k)P_{yy}(\omega_k)-P_{ys}(\omega_k)=[H(\omega_k)P_{yy}(\omega_k)-P_{sy}(\omega_k)]^*=0
    由此得到H(ωk)H(\omega_k)的维纳解:
    H(ωk)=Psy(ωk)Pyy(ωk) H(\omega_k)=\frac{P_{sy}(\omega_k)}{P_{yy}(\omega_k)}
    假设噪声为加性噪声,且与语音信号不相关
    Psy(ωk)=Pss(ωk)Pyy(ωk)=Pxx(ωk)+Pnn(ωk) P_{sy}(\omega_k)=P_{ss}(\omega_k)\\ P_{yy}(\omega_k)=P_{xx}(\omega_k)+P_{nn}(\omega_k)
    上式代入维纳解,我们得到:
    H(ωk)=Pss(ωk)Pss(ωk)+Pnn(ωk) H(\omega_k)=\frac{P_{ss}(\omega_k)}{P_{ss}(\omega_k)+P_{nn}(\omega_k)}
    我们定义这种形式为“频域维纳滤波器”。如果得到H(ωk)H(\omega_k),通过频域相乘,很容易得到S^(ω)=H(ω)Y(ω)\hat{S}(\omega)=H(\omega)Y(\omega)。我们观察频域维纳滤波器的计算公式,只涉及到功率谱的计算,看上去也比时域要简单许多,但是语音信号的短时平稳特性,令我们求真实的输入信号功率谱很麻烦,噪声功率谱也不容易得到。所以需要继续推导寻找近似逼近的方法。我们定义ξωk=ξk=Pss(ωk)Pnn(ωk)\xi_{\omega_k}=\xi_{k}=\frac{P_{ss}(\omega_k)}{P_{nn}(\omega_k)}为频率ωk\omega_k处的先验信噪比SNRpreSNR_{pre},即信号没有被噪声干扰的信噪比。γωk=γk=Pyy(ωk)Pnn(ωk)\gamma_{\omega_k}=\gamma_{k}=\frac{P_{yy}(\omega_k)}{P_{nn}(\omega_k)}为频率ωk\omega_k处的后验信噪比SNRpostSNR_{post},即信号引入加性噪声后的信噪比。则H(ωk)H(\omega_k)可以写成着两种信噪比的表达方法。
    H(ωk)=ξk1+ξk=11γk H(\omega_k)=\frac{\xi_{k}}{1+\xi_{k}}=1-\frac{1}{\gamma_{k}}
    至此,对频域维纳滤波器的求解变成了求解信号的先验或者后验信噪比问题。看上去难度似乎大大降低了,然而实际上对两种信噪比的求解也非常困难。所以有的算法希望利用两种信噪比来平滑计算结果,所以引入一个平滑因子 α\alpha,导出 :
    ξi(k)=αξi(k)+(1α)ξi(k)=αξi(k)+(1α)(γi(k)1)αξi1(k)+(1α)(γi(k)1) \begin{aligned} \xi_i(k)&=\alpha\xi_i(k)+(1-\alpha)\xi_i(k)\\ &=\alpha\xi_i(k)+(1-\alpha)(\gamma_{i}(k)-1)\\ &\approx\alpha\xi_{i-1}(k)+(1-\alpha)(\gamma_{i}(k)-1)\\ \end{aligned}
    因为计算当前帧的先验信噪比是一个非因果事件,所以利用上次滤波后的先验信噪比很显然是一种流行的做法,那么最后我们得出先验信噪比的估计值为:
    ξi^(k)=αξi1(k)+(1α)(γi(k)1) \hat{\xi_i}(k)=\alpha\xi_{i-1}(k)+(1-\alpha)(\gamma_{i}(k)-1)
    进而我们得出当前帧的维纳滤波器
    H^(k)=ξi^(k)1+ξi^(k) \hat{H}(k)=\frac{\hat{\xi_i}(k)}{1+\hat{\xi_i}(k)}
    最后利用此公式实现频域的维纳滤波:
    S^(ω)=H(ω)Y(ω) \hat{S}(\omega)=H(\omega)Y(\omega)
    可以看出,这和时域从维纳解导出自适应滤波的过程还是有些差异的。从公式我们观察两点结论:

    1. 通过该公式我们将降噪问题转化为求解信噪比问题,降低了问题复杂度。
    2. 频域维纳滤波需要准确的估算出先验信噪比和后验信噪比,这两个值得准确程度和收敛速度决定了滤波器以及整个降噪算法的性能。

    WebRtc中的WienerFilter

    以下是代码

    // Estimate prior SNR decision-directed and compute DD based Wiener Filter.
    // Input:
    //   * |magn| is the signal magnitude spectrum estimate.
    // Output:
    //   * |theFilter| is the frequency response of the computed Wiener filter.
    static void ComputeDdBasedWienerFilter(const NoiseSuppressionC* self,
                                           const float* magn,
                                           float* theFilter) {
      size_t i;
      float snrPrior, previousEstimateStsa, currentEstimateStsa;
    
      for (i = 0; i < self->magnLen; i++) {
        // Previous estimate: based on previous frame with gain filter.
        previousEstimateStsa = self->magnPrevProcess[i] /
                               (self->noisePrev[i] + 0.0001f) * self->smooth[i];
        // Post and prior SNR.
        currentEstimateStsa = 0.f;
        if (magn[i] > self->noise[i]) {
          currentEstimateStsa = magn[i] / (self->noise[i] + 0.0001f) - 1.f;
        }
        // DD estimate is sum of two terms: current estimate and previous estimate.
        // Directed decision update of |snrPrior|.
        snrPrior = DD_PR_SNR * previousEstimateStsa +
                   (1.f - DD_PR_SNR) * currentEstimateStsa;
        // Gain filter.
        theFilter[i] = snrPrior / (self->overdrive + snrPrior);
      }  // End of loop over frequencies.
    }
    

    从代码的核心行可以看出这个和上一篇算法推算的原理一致,是利用两个信噪比综合得出当前系数的办法。具体细节处理还有些差异,暂时没有对比,当我们完成对整个vad统计下信噪比的计算之后,再回头看看能否理解每一步公式的细节。

    小结

    本篇学习和分析了WebRtc中的频域维纳滤波器,整个Ns都是围绕着这个来工作的,本文参考了引用里提到的书籍和几篇博文,特此感谢。文中有纰漏的地方非常欢迎纠错。

    引用

    1.《MATLAB在语音信号分析与合成中的应用》宋之用 北京航空航天大学出版社
    2.《一个频域语音降噪算法实现及改进方法》 https://www.cnblogs.com/icoolmedia/p/weiner_audio_ns.html
    3.《WebRTC之noise suppression算法》https://blog.csdn.net/shichaog/article/details/52514816
    4.《Webrtc NS模块算法》https://blog.csdn.net/qq_28882043/article/details/80885240
    5.《webrtc 单通道降噪算法(ANS)简析》https://blog.csdn.net/audio_algorithm/article/details/80812408

    展开全文
  • 卡尔曼滤波和维纳滤波的不同点在于:(1)解决最佳滤波的方法不同,维纳滤波是用频域及传递函数的方法,卡尔曼是用时域及状态变量的方法;(2)维纳滤波要求过程的自相关系数和互相关函数的简单知识,而卡尔曼滤波则...

    完整的实验报告下载连接https://download.csdn.net/download/LIsaWinLee/14884356

    一、实验原理

        卡尔曼滤波和维纳滤波都是最小均方误差为准则的线性估计器。卡尔曼滤波和维纳滤波的不同点在于:(1)解决最佳滤波的方法不同,维纳滤波是用频域及传递函数的方法,卡尔曼是用时域及状态变量的方法;(2)维纳滤波要求过程的自相关系数和互相关函数的简单知识,而卡尔曼滤波则要求时域中状态变量及信号产生过程的详细知识;(3)维纳滤波要求平稳,而卡尔曼滤波则不要求。

    1.1维纳滤波器原理

    维纳滤波的原理是根据全部过去的和当前的观测数据xn,xn-1,… 来估计信号的当前值。以均方误差最小条件下求解系统的传递函数H(z) 或单位脉冲响应h(n) 。维纳滤波的信号模型是从信号与噪声的相关函数得到。

    维纳滤波器的一般结构如下:

    有一期望响应d(n) ,滤波器系数的设计准则是使得滤波器的输出y(n) 是均方意义上对期望响应的最优线性关系。

    线性系统输出为:

    yn=sn=mhmx(n-m)

    均方误差为:Ee2n=Esn-m=0∞hmx(n-m)2

    维纳滤波器的设计,实际上就是在最小均方误差的条件下,即∂E[e2n]hj=0 ,确定滤波器的冲激响应h(n) 或系统函数H(z) ,可等效于求解维纳-霍夫方程。

    1.2卡尔曼滤波器原理

    卡尔曼滤波的原理是不需要全部过去的观察数据,只根据前一个估计值和最近一个观察数据来估计信号的当前值。它是用状态空间法描述系统,即由状态方程和量测方程组成。解是以估计值的形式给出的,其算法是递推。卡尔曼滤波的信号模型(一维)如下:

    离散系统的n维状态方程:x(k)=Akxk-1+wk-1

    离散系统的m维测量方差:yk=ckxk+vk

    Ak 表示状态变量之间的增益矩阵,ck 为状态变量与输出信号之间的增益矩阵,不随时间发生变化,动态噪声wk 与观测噪声vk 都是零均值的正态噪声,且两者互不相关,Rk=var[vk] =E[vkvkT] 为量测噪声协方差矩阵,Qk=var[wk] =E[wkwkT] 为动态噪声协方差矩阵。

    系统初始条件为Ex0=μ0

    varx0 =E[( x-μ0)(x-μ0)T ]=p0

    cov[x0 ,wk]=Ex0wkT=0

    cov[x0 ,vk]=Ex0vkT=0

    卡尔曼滤波的基本思想是先不考虑激励噪声wk 和观测噪声vk ,得到的状态估计值xk 和观测数据的估计值yk ,再用观测数据的估计误差yk=yk-yk 去修正状态的估计值xk ,通过选择修正矩阵H使得状态估计误差的均方值Pk 最小。

    卡尔曼滤波的递推公式如下:

    xk=Akxk-1+Hk(yk-CkAkxk-1)

    Hk=Pk'CkT(CkPk'CkT+Rk)-1

    Pk'=AkPk-1AkT+Qk-1

    Pk=(I-HkCk)Pk'

    假设初始条件Ak,Ck,Qk,Rk,yk ,xk-1,Pk-1 已知,其中x0=Ex0 , P0=var[x0] ,则卡尔曼滤波的递推流程如下

    二、实验内容

    随机信号 服从 过程,它是一个宽带过程,参数如下:

    (1)通过观测方程是方差为1的高斯白噪声,要求分别利用Wiener滤波器和Kalman滤波器通过测量信号估计的波形。

    (2)将 的方差改为4,按(1)中要求,重新估计 的波形。

    自行选择Wiener滤波器的阶数,自行选择实验扩展内容,写出实验报告,并进行相关分析。

    三、实验过程

    3.1维纳滤波器

    利用维纳滤波测量信号时,建立模型如下:

           1)将随机信号X(n)看成是由典型白噪声序列源W(n)激励一个线性系统产生,用一个差分方程表示xn-1.352xn-1+1.338xn-2-0.662xn-3+0.24xn-4=w(n) 。进行Z变换得到Hz=X(z)W(z)=11-1.352z-1+1.338z-2-0.662z-3+0.24z-4 ,均值为1的高斯白噪声序列W(n)可以用randn函数产生,再利用函数X=filter(B,A,W)产生随机信号X(n)。

    2)观测方程是Y(n)=X(n)+V(n),V(n)是方差为1的高斯白噪声,产生进入维纳滤波器的信号。

    维纳滤波仿真流程如下:

    3.2卡尔曼滤波器

    卡尔曼滤波实际由两个过程组成:预测与校正。在预测阶段,滤波器使用上一状态的估计,做出对当前状态的预测。在校正阶段,滤波器利用对当前状态的观测值修正在预测阶段获得的预测值,以获得一个更接进真实值的新估计值。采用最陡下降算法,搜索方向为梯度负方向,每一步更新都使梯度函数值最小。

    卡尔曼滤波仿真流程如下:

    四、实验结果

    4.1 维纳滤波器实验结果

    (1)维纳滤波器的阶数取101阶,观测点数取100,噪声方差为1时,通过仿真得到真实轨迹,观测样本及估计轨迹的比较图像,和通过仿真得到的平均误差如图4.1所示。

            (a)真实轨迹、观测样本及估计轨迹的比较                                                            (b)平均误差

                                                                      图4.1 维纳滤波(噪声方差取1)

    从图4.1(a)可以看出,利用维纳滤波器通过测量信号估计真实信号的波形,得到的估计轨迹接近于真实轨迹。当噪声方差为1时,估计轨迹与真实轨迹间的误差很小。

    从图4.1(b)可以看出,利用维纳滤波器通过测量信号估计真实信号的波形,当噪声方差为1时,得到的估计轨迹与真实轨迹的平均误差较小,而且刚开始的时候平均误差相对较大,随着时间的推移,平均误差有逐渐减小的趋势。

    (2)维纳滤波器的阶数取101阶,观测点数取100,噪声方差为4时,通过仿真得到真实轨迹,观测样本及估计轨迹的比较图像,和通过仿真得到的平均误差如图4.2所示。

                   (a)真实轨迹、观测样本及估计轨迹的比较                                                          (b)平均误差

                                                                                   图4.2 维纳滤波(噪声方差取4)

    从图4.2(a)可以看出,当噪声方差为4时,估计轨迹与真实轨迹间的误差变大了。

    从图4.2(b)可以看出,当噪声方差为4时,估计轨迹与真实轨迹间的平均误差变大了。而且刚开始的时候平均误差相对较大,随着时间的推移,平均误差有逐渐见效的趋势。

    4.2 卡尔曼滤波器实验结果

    (1)卡尔曼滤波器的阶数取101阶,观测点数取100,噪声方差为1时,通过仿真得到真实轨迹,观测样本及估计轨迹的比较图像,和通过仿真得到的平均误差如图4.1所示。

       

                       (a)真实轨迹、观测样本及估计轨迹的比较                                                (b)平均误差

                                                                        图4.3 卡尔曼滤波(噪声方差取1)

    从图4.3(a)可以看出,利用卡尔曼滤波器通过测量信号估计真实信号的波形,得到的估计轨迹接近于真实轨迹。当噪声方差为1时,估计轨迹与真实轨迹间的误差很小。

    从图4.3(b)可以看出,利用卡尔曼滤波器通过测量信号估计真实信号的波形,当噪声方差为1时,得到的估计轨迹与真实轨迹的平均误差较小,而且刚开始的时候平均误差相对较大,随着时间的推移,平均误差有逐渐减小的趋势。

    (2)卡尔曼滤波器的阶数取101阶,观测点数取100,噪声方差为4时,通过仿真得到真实轨迹,观测样本及估计轨迹的比较图像,和通过仿真得到的平均误差如图4.4所示。

       

                        (a)真实轨迹、观测样本及估计轨迹的比较                                                             (b)平均误差

                                                                                   图4.4 卡尔曼滤波(噪声方差取4)

    从图4.4(a)可以看出,当噪声方差为4时,估计轨迹与真实轨迹间的误差变大了。

    从图4.4(b)可以看出,当噪声方差为4时,估计轨迹与真实轨迹间的平均误差变大了。而且刚开始的时候平均误差相对较大,随着时间的推移,平均误差有逐渐见效的趋势。

    五、总结与展望

    通过此次实验,学会运用MATLAB工具实现构建维纳滤波器和卡尔曼滤波器,并对AR模型进行实验。从实验中,我了解到维纳滤波器是以LMS算法为核心的,而卡尔曼滤波器是从RLS算法演变而来的。从这两点则可以衍生出很多解题思路。最后通过MATLAB演示了信号在加入噪声后,通过维纳滤波器和卡尔曼滤波器进行估计,比较估计值和真实值,发现差别非常的小。希望后面若有遇到维纳滤波和卡尔曼滤波问时,能够回忆起课堂上和实验中学到的知识,充分的利用起来。

    六、附录

    6.1 维纳滤波器程序

    (1)噪声方差为1时:

    clc;
    
    clear all;
    
    maxlag=100;
    
    N=100;%观测点数取100
    
    x=zeros(N,1);
    
    y=zeros(N,1);
    
    %%列出状态方程%%
    
    x(1)=randn(1,1);
    
    x(2)=randn(1,1)+1.352*x(1);
    
    x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);
    
    x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);
    
    for n=5:N
    
        x(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1);%x为真实值
    
    end
    
    v=randn(N,1);
    
    y=x+v;%z_x为观测样本值=真值+噪声
    
    %%滤波%%
    
    x=x';
    
    y=y';
    
    xk_s(1)=y(1);%赋初值
    
    xk_s(2)=y(2);
    
    xk_s(3)=y(3);
    
    xk_s(4)=y(4);
    
    xk=[y(1);y(2);y(3);y(4)];
    
    %%维纳滤波器的生成
    
    [rx,lags]=xcorr(y,maxlag,'biased');%观测信号的自相关函数
    
    rx1=toeplitz(rx(101:end));%对称化自相关函数矩阵使之成为方阵,滤波器的阶数为101阶
    
    rx2=xcorr(x,y,maxlag,'biased');%观测信号与期望信号的互相关函数
    
    rx2=rx2(101:end);
    
    h=inv(rx1)*rx2';%维纳霍夫方程
    
    xk_s=filter(h,1,y);%加噪信号通过滤波器后的输出
    
    e_x=0;
    
    eq_x=0;
    
    e_x1=N:1;
    
    %计算滤波的均值,计算滤波误差的均值
    
    for i=1:N
    
        e_x(i)=x(i)-xk_s(i);%误差=真实值-滤波估计值
    
    end
    
    %%作图%%
    
    t=1:N;
    
    figure(1);
    
    plot(t,x,'r-',t,y,'g:',t,xk_s,'b-.');
    
    legend('真是轨迹','观测样本','估计轨迹');
    
    figure(2);
    
    plot(e_x);
    
    legend('平均误差');

    (2)噪声方差为4时:

    clc;
    
    clear all;
    
    maxlag=100;
    
    N=100;%观测点数取100
    
    x=zeros(N,1);
    
    y=zeros(N,1);
    
    %%列出状态方程%%
    
    x(1)=randn(1,1);
    
    x(2)=randn(1,1)+1.352*x(1);
    
    x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);
    
    x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);
    
    for n=5:N
    
        x(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1);%x为真实值
    
    end
    
    v=4*randn(N,1);
    
    y=x+v;%z_x为观测样本值=真值+噪声
    
    %%滤波%%
    
    x=x';
    
    y=y';
    
    xk_s(1)=y(1);%赋初值
    
    xk_s(2)=y(2);
    
    xk_s(3)=y(3);
    
    xk_s(4)=y(4);
    
    xk=[y(1);y(2);y(3);y(4)];
    
    %%维纳滤波器的生成
    
    [rx,lags]=xcorr(y,maxlag,'biased');%观测信号的自相关函数
    
    rx1=toeplitz(rx(101:end));%对称化自相关函数矩阵使之成为方阵,滤波器的阶数为101阶
    
    rx2=xcorr(x,y,maxlag,'biased');%观测信号与期望信号的互相关函数
    
    rx2=rx2(101:end);
    
    h=inv(rx1)*rx2';%维纳霍夫方程
    
    xk_s=filter(h,1,y);%加噪信号通过滤波器后的输出
    
    e_x=0;
    
    eq_x=0;
    
    e_x1=N:1;
    
    %计算滤波的均值,计算滤波误差的均值
    
    for i=1:N
    
        e_x(i)=x(i)-xk_s(i);%误差=真实值-滤波估计值
    
    end
    
    %%作图%%
    
    t=1:N;
    
    figure(1);
    
    plot(t,x,'r-',t,y,'g:',t,xk_s,'b-.');
    
    legend('真是轨迹','观测样本','估计轨迹');
    
    figure(2);
    
    plot(e_x);
    
    legend('平均误差');

    6.2 卡尔曼滤波器程序

    (1)噪声方差为1时:

    N=100;
    
    x=zeros(N,1);
    
    y=zeros(N,1);
    
    var=1;
    
    I=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];
    
    x(1)=randn(1,1);
    
    x(2)=randn(1,1)+1.352*x(1);
    
    x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);
    
    x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);
    
    for n=5:N
    
        x(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1);
    
    end
    
    v=randn(N,1);
    
    y=x+v;
    
    xk_s(1)=y(1);
    
    xk_s(2)=y(2);
    
    xk_s(3)=y(3);
    
    xk_s(4)=y(4);
    
    Ak=[1.352,-1.338,0.662,0.240;1,0,0,0;0,1,0,0;0,0,1,0];
    
    Ck=[1,0,0,0];
    
    Rk=[1];
    
    Pk=[1 0 0 0
    
        0 1 0 0
    
        0 0 1 0
    
        0 0 0 1];
    
    xk=[y(1);y(2);y(3);y(4)];
    
    Qk=[1];
    
    for r=5:N
    
        yk=y(r);
    
        Pk1=Ak*Pk*Ak'+Qk;
    
        Hk=Pk1*Ck'*inv(Ck*Pk1*Ck'+Rk);
    
        xk=Ak*xk+Hk*(yk-Ck*Ak*xk);
    
        Pk=(I-Hk*Ck)*Pk1;
    
        xk_s(r)=xk(1,1);
    
    end
    
    e_x=0;
    
    for i=1:N
    
        e_x(i)=x(i)-xk_s(i);
    
    end
    
    t=1:N;
    
    figure(1);
    
    plot(t,x,'r-',t,y,'g:',t,xk_s,'b-');
    
    legend('真实轨迹','观测样本','估计轨迹');
    
    figure(2);
    
    plot(e_x);
    
    legend('平均误差');

    (2)噪声方差为4时:

    N=100;
    
    x=zeros(N,1);
    
    y=zeros(N,1);
    
    var=2;
    
    I=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];
    
    x(1)=randn(1,1);
    
    x(2)=randn(1,1)+1.352*x(1);
    
    x(3)=randn(1,1)+1.352*x(2)-1.338*x(1);
    
    x(4)=randn(1,1)+1.352*x(3)-1.338*x(2)+0.602*x(1);
    
    for n=5:N
    
        x(n)=1.352*x(n-1)-1.338*x(n-2)+0.602*x(n-3)-0.24*x(n-4)+randn(1,1);
    
    end
    
    v=4*randn(N,1);
    
    y=x+v;
    
    xk_s(1)=y(1);
    
    xk_s(2)=y(2);
    
    xk_s(3)=y(3);
    
    xk_s(4)=y(4);
    
    Ak=[1.352,-1.338,0.662,0.240;1,0,0,0;0,1,0,0;0,0,1,0];
    
    Ck=[1,0,0,0];
    
    Rk=[1];
    
    Pk=[1 0 0 0
    
        0 1 0 0
    
        0 0 1 0
    
        0 0 0 1];
    
    xk=[y(1);y(2);y(3);y(4)];
    
    Qk=[1];
    
    for r=5:N
    
        yk=y(r);
    
        Pk1=Ak*Pk*Ak'+Qk;
    
        Hk=Pk1*Ck'*inv(Ck*Pk1*Ck'+Rk);
    
        xk=Ak*xk+Hk*(yk-Ck*Ak*xk);
    
        Pk=(I-Hk*Ck)*Pk1;
    
        xk_s(r)=xk(1,1);
    
    end
    
    e_x=0;
    
    for i=1:N
    
        e_x(i)=x(i)-xk_s(i);
    
    end
    
    t=1:N;
    
    figure(1);
    
    plot(t,x,'r-',t,y,'g:',t,xk_s,'b-');
    
    legend('真实轨迹','观测样本','估计轨迹');
    
    figure(2);
    
    plot(e_x);
    
    legend('平均误差');

     

    展开全文
  • 频域:最小二乘滤波、维纳滤波; 时域:卡尔曼滤波。 最小二乘滤波与维纳滤波均是对系统的频域特性进行估计,得到系统的传递函数。 卡尔曼滤波对系统的描述,是基于状态空间的描述,即估计系统的内部状态量,...

    卡尔曼滤波、最小二乘法、维纳滤波之我见

    算法要抓住三个方面:模型、策略、求解的方法。

    算法 模型 策略 求解的方法
    最小二乘法 传递函数 误差平方和最小 目标函数导数为0
    卡尔曼滤波 状态空间 均方差最小 先预测,后修正
    维纳滤波 传递函数 均方差最小 目标函数导数为0

    其实,你说最小二乘可以用在状态空间上吗?完全可以,只要知道状态及其导数就可以。但是跟卡尔曼滤波本质不同在于,需要均方差最小,所以可以定义一种概率分布,使二者相等。
    最小二乘法被玩出花的是序贯最小二乘,最容易让人和卡尔曼滤波分不清,因为大家都是迭代求解。但是其实仔细看,其实序贯最小二乘是假设:
    wk+1=[wk,Δ\Deltaw]
    然后寻找Δ\Deltaw使
    min||wk+1x-y||2=min||Δ\Deltaw x-Δ\Deltay||2
    可以看见,与卡尔曼滤波比较:1.这里没有预测;2.这里的策略不同。反而使通过简单变形,就可以变成在线的关于差值系统的最小二乘,因此尽管是不断迭代,还是算作最小二乘方法。

    同时,也有人用卡尔曼滤波解决最小二乘可以解决的问题,当然可以,但是需要先把模型换成状态空间模型,这个是刚需啊。。
    比如拟合曲线,就需要用离散的状态空间模型,只不过:
    xk+1=xk+v
    也就是状态空间的A矩阵变成了单位阵。

    详细一点:
    频域:最小二乘滤波、维纳滤波;
    时域:卡尔曼滤波。

    最小二乘滤波与维纳滤波均是对系统的频域特性进行估计,得到系统的传递函数。
    卡尔曼滤波对系统的描述,是基于状态空间的描述,即估计系统的内部状态量,进而对系统进行描述。

    从方法上讲,最小二乘寻找的是一种投影方法,可以把输出空间投影到输入空间张成的子空间上,并使误差最小(误差与子空间垂直时,最小),关注误差的二范数最小
    维纳滤波关注误差的均方误差,即二范数的均值最小。
    卡尔曼滤波亦为关注误差的均方误差,但是可以适应于非平稳过程。

    展开全文
  • 图象恢复——(逆滤波,维纳滤波

    万次阅读 多人点赞 2018-04-06 10:18:35
    目的:对获取图像在频域用高斯函数进行退化并叠加白噪声,对退化图像进行逆滤波和维纳滤波恢复,比较原始图像和恢复图像,对利用逆滤波和维纳滤波恢复方法恢复图像进行比较。 一、基本原理  图像复原是一种客观的...
  • 维纳滤波一般会在频域中进行分析,因为时域或空域中的卷积在频域中可表示为乘积,为了方便推导一般将信号的频谱展开为一维的列向量,而频响函数则表示为对角矩阵,那么系数的输出就可表示为频响矩阵与信号向量的相乘...
  • 产生模糊冲击函数的模板,对图像尺寸进行延拓,在频域对图像进行模糊,对模糊后的图像叠加高斯噪声,逆滤波,维纳滤波
  • 【图像处理】线性、位置不变退化图像的复原基础(维纳滤波,最小均方滤波,几何滤波)
  • 去相关与维纳滤波

    千次阅读 2011-10-04 16:01:36
    去相关与维纳滤波    说到滤波,我们最容易想到的是频率选择的滤波,比如低通滤波,高通滤波。然后就是FIR与IIR滤波器。维纳滤波器则从另外一个角度来深化了滤波的概念。引用维基百科关于维纳滤波的一段表述如下...
  • 信号检测与处理的一个十分重要的内容就是从噪声中提取信号,本文针对chirp信号加入不同强度的噪声的问题,采用了频域非因果的方法实现了维纳滤波。当伴有信噪比为20的高斯白噪声的chirp信号通过维纳滤波器,它可以对...
  • %维纳滤波 %数值的线性图像复原方法 %维纳滤波器寻找一个是统计误差函数 %e^2=E{(f-f^)^2} %最小估计值 f^ E是期望值操作符 f是未退化的图像 该表达式在频域 %matlab中 维纳滤波是使用函数 deconvwnr 实现的 %fr=...
  • 维纳滤波(1)

    2017-10-24 20:37:00
    当系统中的有效信号和噪声都是随机过程,信号和噪声的频谱还可能重叠(比如有效信号是高斯-马尔可夫过程,噪声是白噪声),根据频域参数设计滤波器的方法就不再适用。 维纳滤波器可以在一些场合解决上述为题,其...
  • 引用维基百科关于维纳滤波的一段表述如下: “仅仅在频域进行滤波的滤波器,仍然会有噪声通过滤波器。维纳设计方法需要额外的关于原始信号所包含频谱以及噪声的信息,维纳滤波器具有以下一些特点: 1、假设...
  • 图像复原4.1.图像退化/复原处理的模型4.2.噪声模型4.2.1.用imnoise函数为图像添加噪声4.2.2....维纳滤波4.8.有约束的最小二乘法滤波4.9.利用露西-理查德森算法的迭代非线性复原4.10.盲去卷积4.11.来自
  • 维纳滤波在图像复原中的应用

    万次阅读 2016-04-25 16:06:17
     g(x,y) = h(x,y)*f(x,y)+n(x,y) 频域:G(u,v) = H(u,v)F(u,v) +N(u,v)   其中f(x,y)为原始图像,h(x,y)为退化函数,n(x,y)为噪声函数,目标就是根据观测图像g(x,y)以及一些先验或者估计信息复原f(x,y)  图像...
  • 基于matlab GUI小波、中值、维纳频域上的滤波 二、源代码 function varargout = dsp1_2(varargin) % DSP1_2 MATLAB code for dsp1_2.fig % DSP1_2, by itself, creates a new DSP1_2 or raises the existing % ...
  • 虽然理论上完全去混响可以通过在某些条件下的逆滤波以及在已知房间脉冲响应(RIR)的条件下实现,但在实践中,RIR的盲辨识在时变和噪声环境中不够精确和鲁棒。因此,成功的去混响方法在时频域得到发展,实际中往往将...
  • 一、简介 在语音去噪中最常用的方法是谱减法,谱减法是一种发展较早且应用较为成熟的语音去噪算法,该算法利用加性噪声与语音不相关的特点,在假设噪声是...转换到频域后,这些峰值听起来就像帧与帧之间频率随机变化
  • 基于matlab GUI小波、中值、维纳频域上的滤波 二、源代码 function varargout = dsp1_2(varargin) % DSP1_2 MATLAB code for dsp1_2.fig % DSP1_2, by itself, creates a new DSP1_2 or raises the existing % ...
  • WienerFilter.m

    2020-03-25 17:06:18
    时域维纳滤波频域维纳滤波,其中频域有直接法和间接法实现对信号的滤波,根据参数flag来选择使用时域维纳滤波还是频域维纳滤波
  • Matlab 实现信号滤波

    千次阅读 2020-11-30 10:07:26
    文章目录项目介绍代码实现1、导入信号2、加入噪声3、绘制原始信号的时域、频域4、滤波4.1 移动平均滤波4.2 中值滤波4.3 维纳滤波4.4 自适应滤波4.5 巴特沃斯滤波4.5.1 低通滤波4.5.2 高通滤波4.5.3 带通滤波 ...
  • LMMSE是频域维纳滤波方法,其减小了MMSE的复杂度,但只适用于慢衰落信道,针对时变信道,本文提出卡尔曼滤波的信道估计方法,仿真结果表明,卡尔曼滤波的信道估计方法在时变信道中具有良好的性能。
  • 自适应中值滤波与维纳滤波都是图像处理...维纳滤波通过将时域图像数据转换至频域并除去降质函数来减少噪声影响,实现平滑效果。文中尝试将上述两种图像处理算法应用于电磁兼容扫描结果的处理,从而改善图像的显示效果。
  • 采用空域迭代反卷积和频域维纳滤波对多纵模激光器光谱进行数据仿真,并在不同光谱仪采样率条件下,比较了迭代反卷积和维纳滤波结果。仿真结果表明,迭代反卷积和维纳滤波可以有效消除光谱仪仪器响应函数引起的光谱...
  • 传统频域LMMSE维纳滤波插值算法性能优异,但其计算复杂度太高而难以实现,针对该问题,提出了一种基于LMMSE的分块滑窗信道估计算法,通过对信道自相关矩阵实行分块确定维纳滤波抽头数,降低了计算复杂度,通过滑窗...
  • 目录Kalman 滤波简介Kalman 滤波的应用领域参考文献 ...经典最优滤波理论包括 Wiener(维纳滤波理论和 Kalman(卡尔曼)滤波理论。前者采用频域方法,后者采用时域空间方法。 采用频域设计法是造成 Wiener 滤...
  • 采取首先对光子晶体光纤的电镜扫描图像使用阈值分割和维纳滤波等图像处理算法获得光纤截面的几何图像,然后使用图像插值及构造均值移动窗口滤波等方法实现了网格划分和介电系数平均等步骤,最后结合频域有限差分法对...
  • 摘 要: 分析了多级维纳滤波器的工作原理和改进方法,引入一种可靠的秩选方法,得到了一种稳妥的多级维纳滤波实现方式。与传统方式相比,引入本秩选方向后的算法很容易找到一个门限,使得输出SINR达到。对算法有限...

空空如也

空空如也

1 2 3 4
收藏数 65
精华内容 26
关键字:

维纳滤波频域