LMS2

LMS fn used in the noise cancellation

Author: tahany682

LMS2.m

@@ -0,0 +1,36 @@
+function [av_e, av_w, w]=LMS2(L,mu,n)%%u is the input signal & n is the number of iteriation
+
+w=zeros(1,L); 
+e=zeros(1,L);
+for mu=mu
+    for i2=1:100
+                 f0=1;
+                 f1=300;
+                 t1=2;
+                 t=0:2./9999:2;
+                 S= 0.9929*chirp(t,f0,t1,f1);
+                 L=length(S);
+                 power_of_signal=0.4965;                                
+                power_of_noise=0.01*power_of_signal;
+                v=sqrt(power_of_noise)*randn(1,1000+L);
+                v=v(1001:1000+L);
+                d=S+v;
+                h = fir1(11, 0.25);
+                % Apply the filter to the input sequence
+                V2 = filter(h, 1, v);
+                u=V2;
+
+              for i=1:n-1
+                        e(i)=d(i)-w(i).*u(i);
+                        w(i+1)=w(i)+mu.*u(i).*e(i);
+              end
+              
+              e_total(:,i2)=e.^2;
+              W_total(:,i2)=w;
+    end
+    av_w = mean(transpose(W_total));    
+    av_e =mean(transpose(e_total));
+    
+end
+ end
+   

SD.m

@@ -0,0 +1,25 @@
+function [Js Ws]=SD(u,d,m,mu)
+
+sgma_d=var(d);
+ 
+[Ps Rs]=correlation(u,d,m);
+
+Ws=zeros(m,1);
+ 
+L_max=max(eig(Rs));
+mu_max=2./L_max;
+
+for mu=mu
+Ws=zeros(m,1);
+for i=1:10000
+    
+ Ws=Ws+mu.*(Ps-(Rs*Ws));
+ Js(i)=sgma_d-transpose(Ws)*Ps-transpose(Ps)*Ws+transpose(Ws)*Rs*Ws;
+ 
+end
+%j=min(Js);
+figure
+plot(Js);
+title('Steepest descent ')
+hold on
+end

correlation.m

@@ -0,0 +1,21 @@
+
+function [P R]=correlation(u,d,m)
+% if length(u)>=length(d)   
+% L=length(d);
+% else
+L=length(u);
+% end
+%CALCULATE R MATRIX
+r=zeros(1,m);
+
+for n=1:m 
+t=sum(u(1:L+1-n).*u(n:L))./length(u(1:L+1-n).*u(n:L));
+r(n)=t;
+end
+
+R=toeplitz(r);
+P=zeros(m,1);
+for n=1:m     
+q=sum(u(1:L+1-n).*d(n:L))./length(u(1:L+1-n).*d(n:L));
+P(n)=q;
+end

noise_cancellation.m

@@ -1,36 +1,62 @@
-% [signal,Fs] = audioread('eric.wav');
-% L=length(signal);
-% m=100;
-% Ts=1/Fs;
-%  windowLength = round(0.01*Fs);
-%  overlapLength = round(0.005*Fs);
-%  S = stft(signal,"Window",hann(windowLength,"periodic"),"OverlapLength",overlapLength,"FrequencyRange","onesided");
-%  S = abs(S);
 
+%-----------------------------------------------------------------------------------------
  f0=1;
  f1=300;
  t1=2;
- t=0:2./9999:2
- S=2*chirp(t,f0,t1,f1);
+ t=0:2./9999:2;
+ S= 0.9929*chirp(t,f0,t1,f1);
  L=length(S);
+ figure
+ plot(S);
+ 
 % power_of_signal=var(signal);
 power_of_signal=var(S);                                
-power_of_noise=0.5*power_of_signal;
-v=sqrt(power_of_noise)*randn(1,1000+L);
-v=v(1001:1000+L);
-b=[1];
+power_of_noise=0.01*power_of_signal;
+v=sqrt(power_of_noise)*randn(1,5000+L);
+v=v(5001:5000+L);
+%% Generate input sequence u(n)
+% Design the half-band LPF
+order = 11;
+f_cutoff = 0.25; % normalized cutoff frequency (0.5 = Nyquist frequency)
+h = fir1(order, f_cutoff);
+% Apply the filter to the input sequence
+V2 = filter(h, 1, v);
+% n=length(v);
+% Plot the results
+% subplot(2,1,1);
+% plot( v);
+% title('Input sequence u(n)');
+% xlabel('n');
+% ylabel('v(n)');
+% 
+% subplot(2,1,2);
+% plot(V2);
+% title('Output sequence v(n)');
+% xlabel('n');
+% ylabel('v(n)');
 
-a=[1 0.5];
-V2=filter(b,a,v);
+%%
+n=4000;
 d=S+v;
 F=100;
-m=11;
-% for i=0:Fs/F           %from 0 to 480 "480 window"
-
-
-[j w]= winner_filter(V2,d,m);
-% 
-% y(i:100)=s(i:F*i)-transpose(j);
-% end
+m=100;
+mu=0.1;
+[j w]= winner_filter(V2,d,m,n);
+figure
+ plot(j);
+ title('winner filter');
+ MMSE_w=min(j)
+%%
+[Js Ws]=SD(V2,d,m,mu);
 %figure
-%stem(w);
+% plot(Js);
+% title('Steepest descent filter');
+MMSE_SD=min(Js)
+p=var(d)
+
+%%
+ 
+ [av_e, av_w, w]=LMS2(L,mu,10000);
+ figure
+ plot(av_e);
+ title("LMS Learning Curves")

winner_filter.m

@@ -0,0 +1,14 @@
+function [j w]= winner_filter(u,d,m,n)
+%calculating the variance of the desired signal
+sgma_d=var(d);
+
+%CALCULATE R MATRIX
+%THE CROSS CORRELATION VECTOR P
+
+[P R]=correlation(u,d,m);
+j=zeros(1,n);
+ for i=1:n
+ w=inv(R)*P;
+ j(i)=sgma_d-transpose(P)*inv(R)*P;
+ end 
+