Eriksson_method.zip_fxLMS的建模_fxlms 实现_fxlms算法_次级建模_次级通道

% a single channel feed-forward active noise control system based on % On-line secondary path modellingthe FxLMS alogrithms % You can find many information about this simulation in "Use of random % noise for online transducer estimate in an adaptive attenuation system " % written by L,J Eriksson and M.A Allie in 1989. % Developed by BoZhong <bozhong316@163.com> % Data: 2012.11.1 % HIT % Harbin, China %-------------------------------------------------------------------------- % ----------------------Let's start the code ------------------------------ clc; clear all; close all; DataLong= 50000; % We do not know P(z) and S(z) in reality. So we have to make dummy paths Tap_Pw = 10; % the number of Primary path filter tap Tap_Sw = 5; % the number of Secondary path filter tap Pw=[0.05 -0.001 0.001 0.8 0.6 -0.2 -0.5 -0.1 0.4 -0.05]; % Primary path filter weight vector Sw=[0.05 -0.01 0.95 0.01 -0.9]; % Secondary path filter weight vector; %----------------online secondary path modelling par----------------------- Tap_ancc = 64; %(Tap_ancc >= Tap_spc in this program) the number of anc controler filter tap size_step_ancc = 0.002; % anc controler size step X_noise=randn(1,DataLong); % control noise Tap_spc = 64; % the number of secondary path modeling controler filter tap size_step_spc = 0.0002; % secondary path modeling controler size step V_iden = sqrt(0.01)*randn(1,DataLong); % Secondary path need online modelling. So, we can generate a white noise signal to online model sp,Generate values from a normal distribution with mean 0 and standard deviation 0.01; %----------------------------reality physics use-------------------------- ANCCx=zeros(1,Tap_ancc); % the state of ANCC(z) ANCCw=zeros(1,Tap_ancc); % the weight of ANCC(z) ActualSPx=zeros(Tap_Sw); % the dummy state for the secondary path ActualSP_Fx=zeros(1,Tap_Sw); % the state of the filtered x(k) % --------------------------Fxlms alogrithm use--------------------------- Shx=zeros(1,Tap_spc); % the state of Sh(z) Shw=zeros(1,Tap_spc); % the weight of Sh(z) s is weight of secondary path Sh_FilterVector=zeros(1,Tap_ancc); % the state of the filtered x(k) % ------------------secondary path identification use---------------------- SPCv=zeros(1,Tap_spc); % the state of SPC(z) SPCw=zeros(1,Tap_spc); % the weight of SPC(z) %bb=zeros(1,Tap_Sw); e_cont=zeros(1,DataLong); % data buffer for the anc controler control error e_sp = zeros(1,DataLong); % data buffer for secondary path modelling error Yp = filter(Pw, 1, X_noise); % and measure the arriving noise at the sensor position epoch=1; while epoch < DataLong %-----------------------calculate anc controler signal--------------- ANCCx = [X_noise(epoch) ANCCx(1:(Tap_ancc-1))]; % update the anc controller state ANCCy = sum(ANCCx.*ANCCw); % calculate the ANC controller output %bb = [V_iden(epoch) bb(1:(Tap_Sw-1))]; %bv = sum(bb.*Sw); MixValue = ANCCy - V_iden(epoch); % secondary loudspeaker output ActualSPx = [MixValue ActualSPx(1:(Tap_Sw-1))]; e_cont(epoch) = Yp(epoch)-sum(ActualSPx.*Sw);%+ bv; % mix signal propagate to secondary path %---------------calculate anc controller filter vector---------------- Shx=[X_noise(epoch) Shx(1:(Tap_spc-1))]; % update the state of Sh(z) Sh_FilterVector=[sum(Shx.*Shw) Sh_FilterVector(1:(Tap_ancc-1))]; % calculate the filtered x(k) %----------------calculate secondary path model error signal--------- SPCv = [V_iden(epoch) SPCv(1:(Tap_spc-1))]; % update the secondary path identifiy control state SPCy = sum(SPCv.*SPCw); % sencondary path filter output e_sp(epoch)=e_cont(epoch)-SPCy; % calculate secondary path modeling error %--------------------------weight updata------------------------------ ANCCw = ANCCw + size_step_ancc*e_cont(epoch)*Sh_FilterVector; % adjust the anc controller weight SPCw = SPCw + size_step_spc*e_sp(epoch)*SPCv; % adjust the secondary path filter weight Shw = SPCw; epoch=epoch+1; end % Lets check the result stem(Sw) hold on stem(Shw, 'r*') ylabel('Amplitude'); xlabel('Numbering of filter tap'); legend('Coefficients of S(z)', 'Coefficients of Sh(z)'); figure(2) freqz(Sw,1); title('Sw(z)的频率响应'); figure(3); freqz(Shw,1); title('Shw(z)的频率响应'); figure(4); subplot(2,1,1); plot([1:DataLong], e_cont); ylabel('Amplitude'); xlabel('Discrete time k'); legend('Noise residue'); subplot(2,1,2); plot([1:DataLong], Yp); hold on ; plot([1:DataLong], Yp-e_cont, 'r'); ylabel('Amplitude'); xlabel('Discrete time k'); legend('Noise signal', 'Control signal');