adaptfxlms.m

Wide-band noise cancellation

close all;
clc;
clear;
%% 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:seyed ahmad javadi num:86448103


%% Active Noise Control

% we shall use a sampling frequency of 8000 Hz.

Fs = 8000;
N = 800;

delayS = 7;

%
fs1=160;
fp1=250;
fp2=1800;
fs2=2000;
App=.1;
Ass=50;
fcuts = [fs1 fp1 fp2 fs2];
mags = [0 1 0];
Gpb=10^(-App/20);
Gsb=10^(-Ass/20);
devs = [Gsb Gpb Gsb];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,Fs);
Hd = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
H = filter(Hd,1,[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;
figure(1);
plot(t,H,'b');
xlabel('Time [sec]');
ylabel('Coefficient value');
title('True Secondary Path Impulse Response');

%% Estimating the Secondary Propagation Path


ntrS = 30000;
s = randn(1,ntrS);
dS = filter(H,1,s) + 0.01*randn(1,ntrS);

%% Designing the Secondary Propagation Path Estimate

M = 250;
muS = 0.1; offsetS = 0.1;
h = adaptfilt.nlms(M,muS,1,offsetS);
[yS,eS] = filter(h,s,dS);

n = 1:ntrS;
figure(2);
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

Hhat = h.Coefficients;
figure(3);
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

delayW = 15;

fs1=200;
fp1=250;
fp2=750;
fs2=800;
App=.1;
Ass=50;
fcuts = [fs1 fp1 fp2 fs2];
mags = [0 1 0];
Gpb=10^(-App/20);
Gsb=10^(-Ass/20);
devs = [Gsb Gpb Gsb];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,Fs);
Hd2 = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale')
G = filter(Hd2,1,[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);
figure(4);
plot(t,G,'b');
xlabel('Time [sec]');
ylabel('Coefficient value');
title('Primary Path Impulse Response');


%% The Noise to be Cancelled.  

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)));
end
d = 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);
figure(5);
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
 

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);
figure(6);
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

Pe = psd(Hp,e(ntrW-20000:ntrW),'NFFT',8192,'Fs',Fs);
figure(7);
plot(Pd.Frequencies,10*log10(Pd.Data(:,1)),'b',...
    Pe.Frequencies,10*log10(Pe.Data(:,1)),'r');
axis([0 2000 -70 0]);
grid on
xlabel('Frequency (Hz)');
ylabel('Power/frequency (dB/Hz)');
tstr = ['Power Spectral Density of the', ...
    ' Original and Attenuated Noise'];
title(tstr);
legend('Original','Attenuated');