1.用MATLAB产生回声的死机源码易语言源代码
2.在matlab中用lu分解求逆矩阵的源代码
用MATLAB产生回声的源代码
clear all
close all
%channel system order
sysorder = 5 ;
% Number of system points
N=;
inp = randn(N,1);
n = randn(N,1);
[b,a] = butter(2,0.);
Gz = tf(b,a,-1);
%This function is submitted to make inverse Z-transform (Matlab central file exchange)
%The first sysorder weight value
%h=ldiv(b,a,sysorder)';
% if you use ldiv this will give h :filter weights to be
h= [0.;
0.;
0.;
0.;
0.;];
y = lsim(Gz,inp);
%add some noise
n = n * std(y)/(*std(n));
d = y + n;
totallength=size(d,1);
%Take points for training
N= ;
%begin of algorithm
w = zeros ( sysorder , 1 ) ;
for n = sysorder : N
u = inp(n:-1:n-sysorder+1) ;
y(n)= w' * u;
e(n) = d(n) - y(n) ;
% Start with big mu for speeding the convergence then slow down to reach the correct weights
if n <
mu=0.;
else
mu=0.;
end
w = w + mu * u * e(n) ;
end
%check of results
for n = N+1 : totallength
u = inp(n:-1:n-sysorder+1) ;
y(n) = w' * u ;
e(n) = d(n) - y(n) ;
end
hold on
plot(d)
plot(y,'r');
title('System output') ;
xlabel('Samples')
ylabel('True and estimated output')
figure
semilogy((abs(e))) ;
title('Error curve') ;
xlabel('Samples')
ylabel('Error value')
figure
plot(h, 'k+')
hold on
plot(w, 'r*')
legend('Actual weights','Estimated weights')
title('Comparison of the actual weights and the estimated weights') ;
axis([0 6 0. 0.])
% RLS 算法
randn('seed', 0) ;
rand('seed', 0) ;
NoOfData = ; % Set no of data points used for training
Order = ; % Set the adaptive filter order
Lambda = 0. ; % Set the forgetting factor
Delta = 0. ; % R initialized to Delta*I
x = randn(NoOfData, 1) ;% Input assumed to be white
h = rand(Order, 1) ; % System picked randomly
d = filter(h, 1, x) ; % Generate output (desired signal)
% Initialize RLS
P = Delta * eye ( Order, Order ) ;
w = zeros ( Order, 1 ) ;
% RLS Adaptation
for n = Order : NoOfData ;
u = x(n:-1:n-Order+1) ;
pi_ = u' * P ;
k = Lambda + pi_ * u ;
K = pi_'/k;
e(n) = d(n) - w' * u ;
w = w + K * e(n) ;
PPrime = K * pi_ ;
P = ( P - PPrime ) / Lambda ;
w_err(n) = norm(h - w) ;
end ;
% Plot results
figure ;
plot(*log(abs(e))) ;
title('Learning Curve') ;
xlabel('Iteration Number') ;
ylabel('Output Estimation Error in dB') ;
figure ;
semilogy(w_err) ;
title('Weight Estimation Error') ;
xlabel('Iteration Number') ;
ylabel('Weight Error in dB') ;
在matlab中用lu分解求逆矩阵的源代码
function X=Ni(A)
%Input - A is an N x N matrix
%Output - I is an N x N inverse matrix of A
%and I(j,:)containing the solution to AX(:,j) =E(:,j).
%Initialize X, Y,the temporary storage matrix C, and the row
% permutation information matrix R
[N,N]=size(A);
B=eye(N); %B is an N x N identity matrix
X=zeros(N,N);
Y=zeros(N,N);
C=zeros(1,N);
R=1:N;
%the next steps is to find the factorization(factorize for only once)
for p=1:N-1
%Find the pivot row for column p
[max1, j]=max(abs(A(p:N,p)));
%Interchange row p and j
C=A(p,:);
A(p,:)=A(j+p-1,:);
A(j+p-1,:)=C;
d=R(p);
R(p)=R(j+p-1);
R(j+p-1)=d;
if A(p,p)==0
'A is singular. No unique solution'
break
end
%Calculate multiplier and place in subdiagonal portion of A
for k=p+1:N
mult=A(k,p)/A(p,p);
A(k,p) = mult;
A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);
end
end
for j=1:N
%when j is fixed then the method is similar to the Program 3.3
%Solve for Y(:,j)
Y(1,j) = B(R(1),j);
for k=2:N
Y(k,j)= B(R(k),j)-A(k,1:k-1)*Y(1:k-1,j);
end
%Solve for X(:,j)
X(N,j)=Y(N,j)/A(N,N);
for k=N-1:-1:1
X(k,j)=(Y(k,j)-A(k,k+1:N)*X(k+1:N,j))/A(k,k);
end
end
