ANCbyFeedbackLMS.rar_ancbyfeedbacklms_fxlms noise_反馈_反馈ANC算法_有源噪

% ANC 反馈控制，参考信号 Xref = e -mu_ref * yw*sh(z) % mu_ref 构建参考信号的修正系数，0则完全用误差信号 %------------------------------------------------------------------------- %反馈控制与前馈控制主要的不同点在于参考信号的获取，前馈控制需要通过声传感器来获取， %而反馈控制则通过线性预测器来产生参考信号。反馈控制适合于窄带噪声的消除。 %问题 %（1） x(k)信号的数量级对系统是否收敛有很大的影响，过大则发散，过小则不能消声 %（2）降噪后的分贝数反而更大，why？ % 系统的框架： % % +-----------+ + % x(k) ------->| P(z) |--yp(k)----------------> sum --+---> e(k) % +-----------+ ^- | % | | % +-------------------------------+ | | % | | | | % | \ | ys(k) | % | +-----------+ | +-----------+ | | % | +--->| C(z) |--yw(k)-+->| S(z) |---+ | % | | +-----------+ +-----------+ | % | | \ | % | | \-----------------\ | % | | \ | % | | +-----------+ +-----------+ | % | +--->| Sh(z) |--xs(k)-+->| LMS |<------+ % | | +-----------+ +-----------+ | % | xh(z) | % | | + | % | +----------------------- sum <-------------------+ % | ^+ % | +-----------+ | % +--------->| Sh(z) |--------+ % +-----------+ % % 与前馈控制类似，由FIR滤波模拟初级通道P(z)，次级通道S(z)，次级通道估计Sh(z)， % 和滤波器C(z)。 % 次级噪声x(k)通过P(z)，施加到声传感器，传感器得到的声信号yp（k）。 % 次级通道的控制信号yw(k)通过S(z)，施加到声传感器，传感器得到的声信号ys（k）。 % 误差信号 e(k) = yp(k) + ys(k) 或 e(k) = yp(k) - ys(k) (实际系统应该是相加的) % e(k) = yp(k) + ys(k) % 使用 “+” 和 “-” 时 注意 参考信号的构建 Xref = e - yw * sh % 滤波器系数更新 w(n+1) = w(n) + u*x(k)*sh(z)*e(k) % % e(k) = yp(k) - ys(k) % 使用 “+” 和 “-” 时 注意 参考信号的构建 Xref = e + yw * sh % 滤波器系数更新 w(n+1) = w(n) + u*x(k)*sh(z)*e(k) % % 采用LMS 算法来更新滤波器的权系数。 % 由于参考信号的处理涉及到次级通道，所以系统需要辨识次级通道的参数sh（z） % 参考文献 "Active Noise Control Systems - % Algorithms and DSP Implementations," % written by S. M. Kuo and % D. R. Morgan in 1996. %-------------------------------------------------------------------------- clear all close all clc % 相关参数的定义 NN=2500; % 数据长度 C = 0; nFilter = 16; % 滤波器的阶数 nSecPath = 16; % 估计的次级通道的阶数， load dataFile; % 变量为 PrimaryData， tt, samplingFrequency % 变更数据，仿真用，加一随机小信号 temp =1.0e-3*sin(2*3.14156*100*tt+3.14156/180*10); PrimaryData = PrimaryData +0.0*temp; temp1=rand(1); temp2=rand(1); temp3=temp1*1e-5*sin(2*3.14156*200*tt+3.14156/180*temp2); PrimaryData=PrimaryData+temp; % 设定初级和次级通道的传递函数 primaryPathWeight=[0.01 0.25 0.5 1 0.5 0.25 0.01]; secondaryPathWeight=primaryPathWeight*0.25; %--------------（1）次级通道参数的辨识（离线辨识，施加一信号，依据其响应确定参数） % 随机产生输入信号 x_iden=randn(1,NN); % 作用与次级通道，得到其响应信号（实际的系统有背景噪声，如何考虑？） y_iden=filter(secondaryPathWeight, 1, x_iden); % 辨识过程 Shx=zeros(1,nSecPath); % the state of Sh(z) Shw=zeros(1,nSecPath); % 估计次级通道的系数 Sh(z) e_iden=zeros(1,NN); % 估计的误差信号 % 应用LMS算法辨识次级通道的参数 mu=0.1; % 参数辨识过程的修正系数（学习率） for k=1:NN, % discrete time k Shx=[x_iden(k) Shx(1:nSecPath-1)]; % update the state Shy=sum(Shx.*Shw); % calculate output of Sh(z) e_iden(k)=y_iden(k)-Shy; % calculate error Shw=Shw+mu*e_iden(k)*Shx; % adjust the weight end % 显示辨识的结果 subplot(2,1,1) plot([1:NN], e_iden) ylabel('Amplitude'); xlabel('Discrete time k'); legend('Identification error'); subplot(2,1,2) stem(secondaryPathWeight) hold on stem(Shw, 'r*') ylabel('Amplitude'); xlabel('Numbering of filter tap'); legend('Coefficients of S(z)', 'Coefficients of Sh(z)') %------------(2) ANC by feedback LMS mu=0.05; % ANC过程的滤波器权系数更新修正系数（学习率） % 初级信号 X=0.1*randn(1,NN)+0.9*cos(2*pi*(1/5)*[0:NN-1]); X = 1.0e2*PrimaryData'; % 通过初级通道到传感器产生的初级噪声信号 primaryNoise=filter(primaryPathWeight, 1, X); % 采用反馈控制，缺少参考信号。需要构建参考信号。 % 给定参考信号初值，长度与初级信号相同 Xh=zeros(size(X)); % 初始化系统 %----------------------------------------------- inputOfFilter=zeros(1,nFilter); % 滤波器的输入，实际应该为参考信号，为计算表达式方便而构造的 filterWeight=zeros(1,nFilter); % 滤波器的权系数； % 滤波器的阶数为16时，可以用以下作为滤波器系数的初值 %filterWeight =[0.4217 0.3852 0.3434 0.3042 0.2808 0.2657 0.2541 0.2416 0.2271 0.2103 0.1933 0.1770 0.1625 0.1491 0.1367 0.1242]; Cyx=zeros(1,nFilter); % the state for the estimate of control signal inputSecPath=zeros(size(secondaryPathWeight)); % 次级通道的输入 inputLMS=zeros(1,nFilter); % LMS的输入（状态） errorNoise=zeros(1,NN); % 声传感器的信号 secondaryNoise=zeros(1,NN); % 次级声源到声传感器的信号 % 控制器的第一个信号 k=1; inputOfFilter=[Xh(k) inputOfFilter(1:nFilter-1)]; % 控制器的输入 outputOfFilter=sum(inputOfFilter.*filterWeight); % 计算控制器的输出 inputSecPath=[outputOfFilter inputSecPath(1:length(inputSecPath)-1)]; % 构建更新次级通道的输入 secondaryNoise(k) = sum(inputSecPath.*secondaryPathWeight); % 次级声源到声传感器的信号 errorNoise(k)=primaryNoise(k)+secondaryNoise(k); % 声传感器的信号 % and start the FbLMS algorithm mu=0.1; % learning rate mu_ref =1.0; %构建参考信号的修正系数，0则完全用误差信号， for k=2:NN % discrete time k Cyx=[outputOfFilter Cyx(1:nSecPath-1)]; % update the state for control signal Xh(k)=errorNoise(k-1)-mu_ref*sum(Cyx.*Shw); % 构建参考信号 inputOfFilter=[Xh(k) inputOfFilter(1:nFilter-1)]; % update the controller state outputOfFilter=sum(inputOfFilter.*filterWeight); % calculate the controller output inputSecPath=[outputOfFilter inputSecPath(1:length(inputSecPath)-1)]; % propagate to secondary path secondaryNoise(k) = sum(inputSecPath.*secondaryPathWeight); errorNoise(k) = primaryNoise(k) + secondaryNoise(k); % measure the residue Shx=[Xh(k) Shx(1:nSecPath-1)]; % update the state of Sh(z) inputLMS=[sum(Shx.*Shw) inputLMS(1:nFilter-1)]; % calculate the filtered x(k) filterWeight=filterWeight - mu*errorNoise(k)*inputLMS; % adjust the controller weight end % 数据的频谱分析 [FFTprimaryNoise,dBAprimaryNoise] = getFFTData(primaryNoise,samplingFrequency,C); [FFTerrorNoise,dBAerrorNoise] = getFFTData(errorNoise,samplingFrequency,C); f = samplingFrequency/NN*[0:(NN-1)]; % 将频率0到samplingFrequency按信号长度分段 fmax = 1000; Nmax = ceil(fmax*NN/samplingFrequency); figure plot(f(1:Nmax),FFTprimaryNoise(1:Nmax),f(1:Nmax),FFTerrorNoise(1:Nmax),'r'); title('声信号的频谱分析'); legend('Primary noise FFT','Noise residue FFT') dBAprimaryNoise dBAerrorNoise