【正文】 R濾波器合適的階數(shù) % 判據(jù)是當階數(shù)增加而均方誤差沒有明顯下降時，則認為階數(shù)足夠 while abs(e0e1)1e6 % e1 和 e0 不夠接近則循環(huán)，確定閾值 N = N+1。 %使 N遞增，使其逐步靠近需要滿足的條件 e0 = e1。 Rxs = Rss(M:(M+N1))。 Rxx = Rww(M:(M+N1))+Rss(M:(M+N1))。 %輸入信號總的自相關函數(shù) R_xx = zeros(N)。 %產生一個 N*N 的全零矩陣作為仿真信號的自相關函數(shù) for j = 1:N for n = 1:N R_xx(j,n) = Rxx(abs(jn)+1)。 end end h = inv (R_xx)*Rxs39。 %求解互相關函數(shù)轉置的逆矩陣 % xsxxRRH 1?? e1 = Rss(M)h39。*Rxs39。 % end N % 顯示 N的最終值 e = e1。 測試信號霍夫方程源程序： % 主程序 clear。 clc。 M = input(39。信號的長度 M = 39。)。 n = 1:M。 %確定數(shù)據(jù)的輸入是 1 到 M s = exp(*n).*sin(pi*n/50)。 % 仿真信號，可以自己生成，任意形式 %q=load(39。D:\39。)。 %s = q(1:M)39。 % load eedata。 % 實際腦電信號 % s = eegdata (1:M)39。 w = *randn(1,M)。 % 白噪聲，系數(shù)代表噪聲相對強度 x = s+w。 % 仿真信號 Rss = xcorr(s,s)。 %估計信號自相關函數(shù) Rww = xcorr(w,w)。 % 估計噪聲自相關函數(shù) [h,e] = WH(Rss,Rww,M)。%自定義的功能函數(shù)，用來求解維納霍夫 ss = filter(h,1,x)。 %用維納濾波器濾波 ,輸入 x 為濾波前序列 % ss 為濾波結果序列， h/1 提供濾波器系數(shù)， h為分子， 1為分母 figure。 subplot(2,2,1)。plot(n,s)。title(39。信號 39。)。 subplot(2,2,2)。plot(n,w)。title(39。噪聲 39。)。 subplot(2,2,3)。plot(n,s+w)。title(39。觀測值 39。)。 subplot(2,2,4)。plot(n,ss)。title(39。信號估計 39。)。 figure。 % 控制窗口數(shù)量，在下一個窗口里面 plot(n,sss)。title(39。估計誤差 39。)。 error=mean((sss).^2) %求最小均方誤差 腦電圖霍夫方程源程序： % 主程序 clear。 clc。 M = input(39。信號的長度 M = 39。)。 n = 1:M。 %s = exp(*n).*sin(pi*n/50)。 % 仿真信號，可以自己生成，任意形式 % load ecgdata。 % 實際腦電信號 %s =ecgdata(1:M)39。 q=load(39。D:\39。)。 s = q(1:M)39。 w = *randn(1,M)。 % 白噪聲，系數(shù)代表噪聲相對強度 x = s+w。 % 仿真信號 Rss = xcorr(s,s)。 %估計信號自相關函數(shù) Rww = xcorr(w,w)。 % 估計噪聲自相關函數(shù) [h,e] = WH(Rss,Rww,M)。 ss = filter(h,1,x)。 %用維納濾波器濾波 figure。 subplot(2,2,1)。plot(4*n,s)。title(39。信號 39。)。 xlabel(39。時間 /ms 39。) 。 ylabel(39。電壓 /uV 39。) %此處每一個格代表一個數(shù)據(jù)點，下同 subplot(2,2,2)。plot(4*n,w)。title(39。噪聲 39。)。 xlabel(39。時間 /ms 39。) 。 ylabel(39。電壓 / uV 39。) subplot(2,2,3)。plot(4*n,x)。title(39。觀測值 39。)。 xlabel(39。時間 /ms 39。) 。 ylabel(39。電壓 / uV 39。) subplot(2,2,4)。plot(4*n,ss)。title(39。信號估計 39。)。 xlabel(39。時間 /ms 39。) 。 ylabel(39。電壓 / uV 39。) figure。 plot(4*n,sss)。title(39。估計誤差 39。)。 xlabel(39。時間 /ms 39。) 。 ylabel(39。電壓 / uV 39。) error=mean((sss).^2) 心電圖霍夫方程源程序： % 主程序 clear。 clc。 M = input(39。信號的長度 M = 39。)。 n = 1:M。 %s = exp(*n).*sin(pi*n/50)。 % 仿真信號，可以自己生成，任意形式 q=load(39。D:\39。)。%腦電數(shù)據(jù)為存在 D盤的 文件 s = q(1:M)39。 % load eedata。 % 實際腦電信號 % s = eegdata (1:M)39。 w = *randn(1,M)。 % 白噪聲，系數(shù)代表噪聲相對強度 x = s+w。 % 仿真信號 Rss = xcorr(s,s)。 %估計信號自相關函數(shù) Rww = xcorr(w,w)。 % 估計噪聲自相關函數(shù) [h,e] = WH(Rss,Rww,M)。 ss = filter(h,1,x)。 %用濾波器濾波 figure。 subplot(2,2,1)。plot(8*n,s)。title(39。信號 39。)。 xlabel(39。時間 /ms 39。) 。 ylabel(39。電壓 /mv 39。) %此處每一個格代表一個數(shù)據(jù)點，下同 subplot(2,2,2)。plot(8*n,w)。title(39。噪聲 39。)。 xlabel(39。時間 /ms 39。) 。 ylabel(39。電壓 /mv 39。) subplot(2,2,3)。plot(8*n,x)。title(39。觀測值 39。)。 xlabel(39。時間 /ms 39。) 。 ylabel(39。電壓 /mv 39。) subplot(2,2,4)。plot(8*n,ss)。title(39。信號估計 39。)。 xlabel(39。時間 /ms 39。) 。 ylabel(39。電壓 /mv 39。) figure。 plot(8*n,sss)。title(39。估計誤差 39。)。 xlabel(39。時間 /ms 39。) 。 ylabel(39。電壓 /mv 39。) error=mean((sss).^2