皮皮网

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