Airplane Noise Removal Demonstration

% Demonstration of LMS algorithm for noise cancellation.
% Rich Kozick, Spring 1997
% Rob Frohne's modifications for Aircraft Noise Cancellation.
% modifications for Macintosh 2000. (commented out here)
% Rob Frohne's modifications for Linux. 
% Aircraft noise cancellation simulation.
% Desired signal
clear all
Totaltime=3;
system("aoss espeak 'Hit a key and speak the signal.'");
Fs=8000; %sample rate.
st=record(Totaltime,Fs);
s = st';
Ls = length(s);
% Interference + random noise
system("aoss espeak 'Hit a key and make the noise!'");
%pause();
nt = record(Totaltime, Fs);
n = nt';
%Sign = 0.01;
%Dn=20;  % Delay of the noise that appears in y.
%n = n(1:Ls) + Sign*randn(Ls,1);
%an = [0 .01 -.5 1 -.5 .1 .01 0];
%bn = 4*[0 0 0 0 0 0 0 0 0 0 0 0 0 0 .5 1 .5];
bn=rand(1,20);
an = [1];
%sys = tf(an,bn,1/Fs);
%bode(sys);
%figure(1);
nf = filter(bn,an,n);
%nf = [nf(Dn:(Ls)); zeros(Dn-1,1)] + Sign*randn(Ls,1);
y = s + nf;
%Sigx = 0.01;  
%ax = [1];       % bx and ax are filtering on n to produce x
%bx = [0 0 0 .1 1 .1];
%Dx = 1;        % Delay of n that appears in x
%x = filter(bx, ax, n);%x = [x(Dx:(Ls)); zeros(Dx-1,1)] + Sigx*randn(Ls,1);
%x = x + Sigx*randn(1,Ls);

system("aoss espeak 'Here is the noisy signal.'");
soundsc(y,Fs);

N = 20;         % Length of adaptive filter

% LMS algorithm for adaptive noise cancellation

h = zeros(N,1);
%mu = 1/(10*N*var(n));
mu = .05;
for k=N:Ls
   xk = n(k:-1:(k-N+1));
   nhat(k) = h'*xk';
   e(k) = - y(k) + nhat(k);
   h = h - mu*e(k)*xk'/(xk*xk'); % For the previous line.
end

% The signal estimate is in the vector e
%speak('Here is the cleaned signal.');
system("aoss espeak 'Here is the cleaned signal.'");
soundsc(e(1,6000:end),Fs);
%playaudio(e);  (For Linux without octave-forge.)

The record.m function in Ubuntu Maverick is broken. You can use this one instead.

## Copyright (C) 1995, 1996, 1997, 1998, 1999, 2000, 2003, 2004, 2005,
##               2006, 2007 John W. Eaton
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.   

## -*- texinfo -*-
## @deftypefn {Function File} {} record (@var{sec}, @var{sampling_rate})
## Records @var{sec} seconds of audio input into the vector @var{x}.  The
## default value for @var{sampling_rate} is 8000 samples per second, or
## 8kHz.  The program waits until the user types @key{RET} and then
## immediately starts to record.
## @seealso{lin2mu, mu2lin, loadaudio, saveaudio, playaudio, setaudio}
## @end deftypefn

## Author: AW <Andreas.Weingessel@ci.tuwien.ac.at>
## Created: 19 September 1994
## Adapted-By: jwe
## And adapted again 11/25/2010 by Rob Frohne

function X = record (sec, sampling_rate)

  if (nargin == 1)
    sampling_rate = 8000;
  elseif (nargin != 2)
    print_usage ();
  endif

  file = tmpnam ();
  file= [file,".wav"];

  input ("Please hit ENTER and speak afterwards!\n", 1);

  cmd = sprintf ("rec -c1 -r%d %s trim 0 %d",
                   sampling_rate, file, sec)

  system (cmd);

  X = wavread(file);

endfunction