%-----------------------------------------------------------------------------
% Author - Phillip De Leon
% Date - 10/11/96
% Required Files - SP_Toolbox, tf_plot.m, lms.m, nlms.m, speech.mat
% Reference - Haykin p. 387
%-----------------------------------------------------------------------------
% ALE.M - this code simulates an adaptive line enhancer as described in the
%   above reference.  The idea is to adaptively cancel periodic interference
%   in a broadband signal.
%-----------------------------------------------------------------------------
clear
pack
clg

%
% User params
%
speech = 1;  % 1) speech broadband signal 0) white noise broadband signal
f = 1000;  % frequency of interference (1000)
M = 8;   % length of adaptive filter (8)
adaptive_algorithm = 'nlms';   % either 'lms' or 'nlms'
delta = 1000;  % decorrelation delay (in samples) (1000 for speech, 10 for white noise)
SNR = 100.0;   % SNR at the input to the ALE (sinusoid power/broadband power)
spectral_window_length = 2048;  % length of FFT for spectral analysis

if (strcmp(adaptive_algorithm,'lms'))
   mu = 1/(120000*M);  % (1/(12000*M),1/(120000*M) etc...) PLAY WITH THESE!!!
elseif (strcmp(adaptive_algorithm,'nlms'))
   mu_tilde = 0.0005;  % (1,0.5,0.05,0.0005) PLAY WITH THESE!!!
else
   error('Error: Adaptive algorithm must be either ''lms'' or ''nlms.''')
end;
w_hat = zeros(M,1); 

%
% Load broadband signal
%
if (speech)
   load speech;  % load speech file and sampling rate 
else
   fs = 8000;  % sample rate (8000)
   x = randn(5*fs,1);
end;
N = length(x);  % length of broadband signal

%
% Build periodic interference
%
p = sqrt(SNR)*sinusoid(fs,f,N/fs);
p = p(1:N);  % make sure sinusoid is same length as broadband signal

%
% Add broadband signal + periodic interference, plot signal and spectrum, and play it
%
u = x + p;
U = length(u);
time_freq_plot(u(U-spectral_window_length+1:U),fs);
subplot(211);title('Broadband signal plus periodic interference');ylabel('u(k)');
subplot(212);title(['Signal spectrum (',num2str(spectral_window_length),' samples)']);ylabel('U(w)');
disp('Hit any key to continue...');
pause;
disp('Playing broadband signal plus sinusoidal interference...')
sound(u,fs);
disp('Now cancelling sinusoidal interference...')

%
% Build delayed version of u for input to ALE
%
u_delay = u(delta:N);
u = u(1:N-delta+1);  % make signals the same length

%
% Cancel interference
%
if strcmp(adaptive_algorithm,'lms')
   [e,w_hat] = lms(w_hat,u_delay,u,mu);
else
   [e,w_hat] = nlms(w_hat,u_delay,u,mu_tilde);
end;

%
% Plot error signal and spectrum, play e(n) = u(n) - y(n) = u(n) - ALE output
%
disp('Playing broadband signal plus sinusoidal interference after ALE processing...')
sound(e,fs);
E = length(e);
time_freq_plot(e(E-spectral_window_length+1:E),fs);
subplot(211);title('Broadband signal plus periodic interference (after ALE)');ylabel('e(k)');
subplot(212);title(['Error spectrum (',num2str(spectral_window_length),' samples)']);ylabel('E(f)');
