anc.zip_FXLMS 噪声_anc matlab_fxlms noise_fxlms算法_有源噪声控制

% ANC based on FxLMS % u(k) +--------+ d(k) + e(k) % U ------+--->| Pw |------D----------------------------------------------------> sum --+--->e % | +--------+ ^+ | % | | | % | \ r(k)| | % | +---------+ y(k) + +---------+ | | % +-Wx->| Ww |Wy-------------->sum-----------------------Sx>| Sw |Sy-+ | % | +---------+ ^+ +---------+ | % | \ +-------+ | +-------+ - + | % | \ | WGN |V---+----Vx->| Cw_0 |Vy------->sum<---------------+ % | \ +-------+ | +-------+ | | % | \ | | | % | \ | | | % | \ | +---------+ | | % | \ +----Vx>| LMS |<f--------+ | % | \-----------\ +---------+ | % | \ | % | \ | % | +---------+ x(k) +-----------+ | % +-Cx->| Cw_1 |Cy--Xx--->| LMS |<-----------------------------------------+ % +---------+ +-----------+ % use FIR to model P(z),W(z),S(z),C(z) clear all clc fs = 10000;%采样频率 L = 3000;% 长度 m = 0:L-1; t = m/fs; %时间 N = 7; % order of primary and secondary fir filter M = 17; % order of adaptive fir filter W(z),M 应取 N 的2-3之间 %noise sourse % 输入的估计 Z = 5; % 叶片数 n = 60; % 风扇转速(r/s) fk = n*Z; % 旋转噪声信号基频(Hz)(500Hz以下为低频) R = 0.1; % 叶片长度(m) Vm = 2*pi*fk*R; % 风扇圆周线速度最大值 h = 0.2; % 常值，0.15-0.22 d = 0.01; % 叶片在气流入射方向的厚度(m) fm = h*Vm/d; % 絮流噪声频率峰值 % 风扇噪声由旋转噪声和絮流噪声组成，旋转噪声用正弦信号代表 % 絮流噪声频率峰值为fm，用白噪声代表 u1 = sin(2*pi*fk*t)+0.5*sin(4*pi*fk*t)+0.2*sin(6*pi*fk*t); % 旋转噪声 u2 = wgn(1,L,5); % 加入白噪声 % 用低通滤波器将高于fm的噪声滤除 fk1 = fk*2/fs; fm1 = fm*2/fs; [order1,wn] = buttord(fm1,1.2*fm1,3,30); [b1,a1] = butter(order1,wn,'low'); U1 = filter(b1,a1,u2); U = u1+U1; % 在线辨识的辨识信号输入,白噪声 V = randn(1,L); %设置p(z)和s(z)的权值，以下用fir滤波器进行拟合 Pw1 = [0 1 0 -0.07 0.006]; Pw2 = 1; Sw = [0 1 4.4 3 -3.6 0 0 0]; % Sw = rand(1,N+1); % D = filter(Pw1, Pw2, U); % r = filter(Sw, 1, V); % % NFFT = 2^nextpow2(L); % f_1 = fs/2*linspace(0,1,NFFT/2+1); % y_1 = fft(U,NFFT)/L; % figure; % plot(f_1,2*abs(y_1(1:NFFT/2+1))) % xlabel('Frequency (Hz)') % ylabel('|y_1(f)|') % y_1 = fft(V,NFFT)/L; % figure; % plot(f_1,2*abs(y_1(1:NFFT/2+1))) % xlabel('Frequency (Hz)') % ylabel('|y_1(f)|') % y_1 = fft(D,NFFT)/L; % figure; % plot(f_1,2*abs(y_1(1:NFFT/2+1))) % xlabel('Frequency (Hz)') % ylabel('|y_1(f)|') % % 初始化 Cx = zeros(1,N+1); % the state of C(z) Vx = zeros(1,N+1); Cw_0 = zeros(1,N+1); % the weight of C(z) % Cw_1 = ones(1,N+1); Wx = zeros(1,M); % the state of C(z) Ww = zeros(1,M); % the weight of C(z) Sx = zeros(size(Cw_0)); % the dummy state for the secondary path e = zeros(1,L); % data buffer for the control error f = zeros(1,L); % WGN noise error Xx = zeros(1,M); % the state of the filtered u(k) for k = 1:L, Wx = [U(k) Wx(1:M-1)]; Wy = sum(Wx.*Ww); Vx = [V(k) Vx(1:N)]; Vy = sum(Vx.*Cw_0); Sx = [Wy Sx(1:length(Sx)-1)]+Vx; e(k) = D(k)-sum(Sx.*Sw); f(k) = e(k)-Vy; deta_iden = Vx*Vx'; mu_iden = 0.001/deta_iden; Cw_0 = Cw_0-mu_iden*(f(k)+D(k)+Wy)*Vx; % Cw_0 = Cw_0-mu_iden*f(k)*Vx; Cx = [U(k) Cx(1:N)]; % update the state Cy = sum(Cx.*Cw_0); % calculate output of C(z) Xx = [Cy Xx(1:M-1)]; % deta_main_1 = Wx*Wx'; % deta_main_2 = Xx*Xx'; % % if deta_main_1 <= deta_main_2 % mu_main = 0.005/deta_main_2; % else % mu_main = 0.005/deta_main_1; % end mu_main = 0.000001; Ww = Ww-mu_main*e(k)*Xx; % if k<N+1 % Cw_1 = [Cw_0(1) Cw_1(1:N)]; % else % Cw_1 = Cw_0; % end end % 检验结果 figure subplot(2,2,1) plot([1:L], f) ylabel('Amplitude'); xlabel('Discrete time k'); legend('Identification error'); subplot(2,2,2) stem(Sw, 'g*') hold on stem(Cw_0, 'r<') % hold on % stem(Cw_1, 'b>') ylabel('Amplitude'); xlabel('Numbering of filter tap'); legend('Coefficients of S(z)','Coefficients of C0(z)'); subplot(2,2,3) plot([1:L], e) ylabel('Amplitude'); xlabel('Discrete time k'); legend('Noise residue') subplot(2,2,4) plot([1:L], D) hold on plot([1:L], D-e, 'r:') ylabel('Amplitude'); xlabel('Discrete time k'); legend('Noise signal', 'Control signal')