%% Active Noise Control Using a Filtered-X LMS FIR Adaptive Filter.% This demonstration illustrates the application of adaptive filters to the% attenuation of acoustic noise via active noise control.%%   Author(s): Scott C. Douglas%   Copyright 1999-2005 The MathWorks, Inc.%% Active Noise Control% In active noise control, one attempts to reduce the volume of an unwanted% noise propagating through the air using an electro-acoustic system using% measurement sensors such as microphones and output actuators such as% loudspeakers.  The noise signal usually comes from some device, such as a% rotating machine, so that it is possible to measure the noise near its% source. The goal of the active noise control system is to produce an% "anti-noise" that attenuates the unwanted noise in a  desired quiet% region using an adaptive filter.  This problem differs from traditional% adaptive noise cancellation in that:%        - The desired response signal cannot be directly measured; %          only the attenuated signal is available.%        - The active noise control system must take into account the%          secondary loudspeaker-to-microphone error path in its %          adaptation.%% For more implementation details on active noise control tasks, see% S.M. Kuo and D.R. Morgan, "Active Noise Control Systems:  Algorithms% and DSP Implementations",  Wiley-Interscience, New York, 1996.%% The Secondary Propagation Path% The secondary propagation path is the path the anti-noise takes from% the output loudspeaker to the error microphone within the quiet zone.% The following commands generate a loudspeaker-to-error microphone% impulse response that is bandlimited to the range 160 - 2000 Hz and% with a filter length of 0.1 seconds. For this active noise control task, % we shall use a sampling frequency of 8000 Hz.Fs = 8000;N = 800;delayS = 7;Fd = fdesign.bandpass('N,Fst1,Fst2,Ast',8,0.04,0.5,20);Hd = design(Fd,'cheby2','FilterStructure','df2tsos');H = filter(Hd,[zeros(1,delayS) log(0.99*rand(1,N-delayS)+0.01).* ...    sign(randn(1,N-delayS)).*exp(-0.01*(1:N-delayS))]);H = H/norm(H);t = 1/Fs:1/Fs:N/Fs;plot(t,H,'b');xlabel('Time [sec]');ylabel('Coefficient value');title('True Secondary Path Impulse Response');%% Estimating the Secondary Propagation Path% The first task in active noise control is to estimate the impulse% response of the secondary propagation path.  This step is usually % performed prior to noise control using a synthetic random signal played % through the output loudspeaker while the unwanted noise is not present. % The following commands generate 3.75 seconds of this random noise as well % as the measured signal at the error microphone.  ntrS = 30000;s = randn(1,ntrS);dS = filter(H,1,s) + 0.01*randn(1,ntrS);%% Designing the Secondary Propagation Path Estimate% Typically, the length of the secondary path filter estimate is not as% long as the actual secondary path and need not be for adequate control% in most cases.  We shall use a secondary path filter length of 250% taps, corresponding to an impulse response length of 31 msec. % While any adaptive FIR filtering algorithm could be used for this % purpose, the normalized LMS algorithm is often used due to its % simplicity and robustness. Plots of the output and error signals show% that the algorithm converges after about 10000 iterations.M = 250;muS = 0.1; offsetS = 0.1;h = adaptfilt.nlms(M,muS,1,offsetS);[yS,eS] = filter(h,s,dS);n = 1:ntrS;plot(n,dS,n,yS,n,eS);xlabel('Number of iterations');ylabel('Signal value');title('Secondary Identification Using the NLMS Adaptive Filter');legend('Desired Signal','Output Signal','Error Signal');%% Accuracy of the Secondary Path Estimate% How accurate is the secondary path impulse response estimate?  This% plot shows the coefficients of both the true and estimated path.% Only the tail of the true impulse response is not estimated % accurately.  This residual error does not significantly harm the% performance of the active noise control system during its operation% in the chosen task. Hhat = h.Coefficients;plot(t,H,t(1:M),Hhat,t,[H(1:M)-Hhat(1:M) H(M+1:N)]);xlabel('Time [sec]');ylabel('Coefficient value');title('Secondary Path Impulse Response Estimation');legend('True','Estimated','Error');%% The Primary Propagation Path% The propagation path of the noise to be cancelled can also be% characterized by a linear filter.  The following commands% generate an input-to-error microphone impulse response that is% bandlimited to the range 200 - 800 Hz and has a filter length of % 0.1 seconds.delayW = 15;Fd2 = fdesign.bandpass('N,Fst1,Fst2,Ast',10,0.05,0.2,20);Hd2 = design(Fd2,'cheby2','FilterStructure','df2tsos');G = filter(Hd2,[zeros(1,delayW) log(0.99*rand(1,N-delayW)+0.01).*...    sign(randn(1,N-delayW)).*exp(-0.01*(1:N-delayW))]);G = G/norm(G);plot(t,G,'b');xlabel('Time [sec]');ylabel('Coefficient value');title('Primary Path Impulse Response');%% The Noise to be Cancelled% Typical active noise control applications involve the sounds of% rotating machinery due to their annoying characteristics.  Here,% we have synthetically generated 7.5 seconds of a noise that % might come from a typical electric motor.  Listening to its sound % at the error microphone before cancellation, it has the characteristic % industrial "whine" of such motors.  The spectrum of the sound% is also plotted.  ntrW = 60000;F0 = 60;n = 1:ntrW;A = [0.01 0.01 0.02 0.2 0.3 0.4 0.3 0.2 0.1 0.07 0.02 0.01];x = zeros(1,ntrW);for k=1:length(A);    x = x + A(k)*sin(2*pi*(F0*k/Fs*n+rand(1)));endd = filter(G,1,x) + 0.1*randn(1,ntrW);Hp = spectrum.welch; Hp.SegmentLength = 4444;Hp.FFTLength='UserDefined';Pd = psd(Hp,d(ntrW-20000:ntrW),'NFFT',8192,'Fs',Fs);plot(Pd)axis([0 2 -70 0]);title('Power Spectral Density of the Noise to be Cancelled');p8K = audioplayer(d/max(abs(d)),Fs);playblocking(p8K);%% Active Noise Control using the filtered-X LMS Algorithm% The most popular adaptive algorithm for active noise control is% the filtered-X LMS algorithm.  This algorithm uses the secondary% path estimate to calculate an output signal whose contribution% at the error sensor destructively interferes with the undesired% noise.  The reference signal is a noisy version of the undesired% sound measured near its source.  We shall use a controller filter% length of about 44 msec and a step size of 0.0001 for these% signal statistics.  The resulting algorithm converges after about% 5 seconds of adaptation. Listening to the error signal, the % annoying "whine" is reduced considerably.  xhat = x + 0.1*randn(1,ntrW);L = 350;muW = 0.0001;h = adaptfilt.filtxlms(L,muW,1,Hhat);[y,e] = filter(h,xhat,d);plot(n,d,'b',n,y,'g',n,e,'r');xlabel('Number of iterations');ylabel('Signal value');tstr = ['Active Noise Control Using', ...    ' the Filtered-X LMS Adaptive Controller'];title(tstr);legend('Original Noise','Anti-Noise','Residual Noise');p8K = audioplayer(e/max(abs(e)),Fs);playblocking(p8K);%% Residual Error Signal Spectrum% Comparing the spectrum of the residual error signal with that of the% original noise signal, we see that most of the periodic components % have been attenuated considerably.  The steady-state cancellation % performance may not be uniform across all frequencies, however.  % Such is often the case for real-world systems applied to active noise % control tasks.Pe = psd(Hp,e(ntrW-20000:ntrW),'NFFT',8192,'Fs',Fs);plot(Pd.Frequencies,10*log10(Pd.Data(:,1)),'b',...    Pe.Frequencies,10*log10(Pe.Data(:,1)),'r');axis([0 2000 -70 0]);grid onxlabel('Frequency (Hz)');ylabel('Power/frequency (dB/Hz)');tstr = ['Power Spectral Density of the', ...    ' Original and Attenuated Noise'];title(tstr);legend('Original','Attenuated');displayEndOfDemoMessage(mfilename)