new-lms-arithmetic-code .m

In 1960, R.E. Kalman published his famous paper describing a recursive solution to the discrete-dat

clear all;close all;
RF = load('RFrost.txt');

refGain = 1;
worder = 8;


%N = 2048;
%t=1:N;
%signal = sin(2*pi.*t.*t/N/N*8); %.*fliplr(cos(2*pi.*t.*t/N/N*15))
%signal = 0*t;
%signal = randn(1,N);
%signal = RF(3000:7096)';
%signal = 0*t;
signal = RF';
%clear RF;

N = length(signal);

wreal = randn(1,worder);
additivenoise = randn(1,N);
%load additivenoise;

ref = conv(additivenoise,wreal);
primary = signal + ref(1:length(signal));
fref = additivenoise*refGain; %real reference mic


w(1,:) = ones(1,worder);
mu = .1;
%Zero pad so we can start filter at 0 and not throw of the index
frefpad = [zeros(1,worder -1) fref];


start = flops;
for n = 1:N;
   %offset n so we can reference the correct value in zero-padded fref
   m = n + worder -1;
   frefblock = frefpad(m-worder+1:1:m)';
   refP(n) = w(n,:)*(frefblock);   %ADDING EXTRA RANDOM NOISE TO REFERENCE MIC
   output(n) = primary(n) - refP(n);
   w(n+1,:) = w(n,:) + mu.*frefblock'.*output(n); %we are using the output as our error signal
end;

work = flops-start;
w(length(w),:);
WErrorNorm = norm(wreal - w(length(w),:))


%*Plot of w vs. time
my = [max(max(w)), min(min(w))];
mx = [worder-1, worder+1];
%plot(mx,my,'g')
%figure;
for ii = 1:worder;
   %rv = wreal(ii)*ones(length(w));
   %rv = wreal(:,ii);
  % plot(rv,'blue');
   %hold on
   %plot(w(:,ii),'r');
end;


figure;
subplot(2,2,1);
plot(primary,'r');
axis([0 length(primary) min(primary) max(primary)]);
title('primary microphone signal');
subplot(2,2,2);
plot(output,'y');
axis([0 length(primary) min(signal)-.5 max(signal)+.5]);
title('filtered output');
subplot(2,2,3);
plot((ref(1:length(refP))-refP).^2,'b');
axis([0 length(primary) min(primary) max(primary)]);
title('error');
subplot(2,2,4);
plot(signal,'g');
axis([0 length(signal) min(signal) max(signal)]);
title('signal');

%We start calculating the noise at 4 since the early values of the output
%can be VERY large, and bias our SNR measurement.  


sv = 2*worder;
sw = length(signal);

SNRpre = norm(signal(sv:sw))/norm(ref(sv:sw))
%SNRpre = 10*log10(SNRpre)

SNRpost = norm(signal(sv:sw))/norm(signal(sv:sw) - output(sv:sw))
%SNRpost = 10*log10(SNRpost)


wavwrite(signal,10000,8,'signal.wav');
wavwrite(primary,10000,8,'primary.wav');
wavwrite(fref,10000,8,'reference.wav');
wavwrite(output,10000,8,'output.wav');