close all
%channel system order
sysorder = 5
% Number of system points
N=2000
inp = randn(N,1)
n = randn(N,1)
[b,a] = butter(2,0.25)
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.0976
0.2873
0.3360
0.2210
0.0964]
y = lsim(Gz,inp)
%add some noise
n = n * std(y)/(10*std(n))
d = y + n
totallength=size(d,1)
%Take 60 points for training
N=60
%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 <20
mu=0.32
else
mu=0.15
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.05 0.35])
% RLS 算法
randn('seed', 0)
rand('seed', 0)
NoOfData = 8000 % Set no of data points used for training
Order = 32 % Set the adaptive filter order
Lambda = 0.98 % Set the forgetting factor
Delta = 0.001 % 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(20*log10(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')
我恰好也做此题啊!你运气真好!下面的是C语言写的,不知是否可行?
那矩阵的名字可能不符合你的要求,你自己改一下吧,呵呵~~~~~
#include<stdio.h>
#define MAXSIZE 100 /* 非零元个数的最大值 */
typedef struct triple
{
int i,j/* 行下标,列下标 */
int e/* 非零元素值 */
}triple
typedef struct tsmatrix
{
triple data[MAXSIZE+1]/* 非零元三元组表,data[0]未用 */
int mu,nu,tu/* 矩阵的行数、列数和非零元个数 */
/* 各列第一个非零元的位置表rpos[0]未用 */
}rlsmatrix
createsmatrix(rlsmatrix *M)
{ /* 创建稀疏矩阵M */
int e,i,m,n
M->data[0].i=0/* 为以下比较顺序做准备 */
printf("请输入矩阵的行数，列数，和非零元素的个数:")
scanf("%d",&M->mu)scanf("%d",&M->nu)scanf("%d",&M->tu)
for(i=1i<=M->tui++)
{
printf("请按行序顺序输入第%d个非零元素所在的行(1～%d),列(1～%d),元素值：",i,M->mu,M->nu)
scanf("%d",&m)scanf("%d",&n)scanf("%d",&e)
if(m<1||m>M->mu||n<1||n>M->nu) /*行或列超出范围 */
{printf("行或列超出范围")getch()exit()}
if(m<M->data[i-1].i||m==M->data[i-1].i&&n<=M->data[i-1].j) /*行或列的顺序有错*/
{printf("行或列的顺序有错")getch()exit()}
M->data[i].i=m
M->data[i].j=n
M->data[i].e=e
}
}
/* 求矩阵的快速转置 */
void transposesmatrix(rlsmatrix M,rlsmatrix *T)
{ /* cpos存放每列的第一个非零元素的地址，temp中间变量 */
int i,m,*cpos,*temp,k=0
T->mu=M.nu
T->nu=M.mu
T->tu=M.tu
cpos=(int *)malloc(M.mu*sizeof(int))
if(cpos==NULL)exit()
temp=(int *)malloc(M.mu*sizeof(int))
if(temp==NULL)exit()
/* 对cpos对初始化，初值为0 */
*(cpos+1)=0
for(i=1i<=M.nui++)
{
for(m=1m<=M.tum++)
{
if(M.data[m].j==i)
k++
}
temp[i]=k
if(i==1&&k!=0)
*(cpos+i)=1/* 为cpos赋值 */
if(i>1)
*(cpos+i)=*(temp+i-1)+1
}
free(temp)
for(i=1i<=M.tui++)/* 进行转置 */
{T->data[*(cpos+M.data[i].j)].i=M.data[i].j
T->data[*(cpos+M.data[i].j)].j=M.data[i].i
T->data[*(cpos+M.data[i].j)].e=M.data[i].e
(*(cpos+M.data[i].j))++}
free(cpos)
}
multsmatrix(rlsmatrix M,rlsmatrix N,rlsmatrix *T)
{
int i,j,Qn=0
int *Qe
if(M.nu!=N.mu)
{printf("两矩阵无法相乘")getch()exit()}
T->mu=M.mu
T->nu=N.nu
Qe=(int *)malloc(M.mu*N.nu*sizeof(int))/* Qe为矩阵Q的临时数组*/
for(i=1i<=M.mu*N.nui++)
*(Qe+i)=0/* 矩阵Q的第i行j列的元素值存于*(Qe+(M.data[i].i-1)*N.nu+N.data[j].j)中，初值为0 */
for(i=1i<=M.tui++) /* 结果累加到Qe */
for(j=1j<=N.tuj++)
if(M.data[i].j==N.data[j].i)
*(Qe+(M.data[i].i-1)*N.nu+N.data[j].j)+=M.data[i].e*N.data[j].e
for(i=1i<=M.mui++) /*Qe矩阵中，因为M的每一行和N的每一列相乘都是T的一个元素，不管它是零或非零 */
for(j=1j<=N.nuj++) /*当M的第一行和N的第一列相乘则得T的第一个元素；当M的第一行和N的第二列相乘则得T的第二个元素；…… */
if(*(Qe+(i-1)*N.nu+j)!=0) /*……当M的第i行和N的第j列相乘则得T的第p个元素；根据归纳法得p=(i-1)*N的列数+j */
{
Qn++//非零元个数加一
T->data[Qn].e=*(Qe+(i-1)*N.nu+j)
T->data[Qn].i=i
T->data[Qn].j=j
}
free(Qe)
T->tu=Qn
return
}
void printmatrix(rlsmatrix M)//输出
{
int i,m=1,n,k
printf("矩阵的简化模式是:\n")
for(i=1i<=M.tui++)
printf("%d,%d,%d\n",M.data[i].i,M.data[i].j,M.data[i].e)
printf("矩阵的行数是:%d\n",M.mu)
printf("矩阵的列数是:%d\n",M.nu)
printf("矩阵中非零元素个数是:%d\n",M.tu)
printf("矩阵的完整模式的:\n")
for(n=1n<=M.mun++)
{
printf("|")
for(k=1k<=M.nuk++)
{
if(M.data[m].i==n&&M.data[m].j==k)
{printf("%3d",M.data[m].e)
m++}
else
printf(" 0")
}
printf(" |\n")
}
return
}
void main()
{
rlsmatrix M,N,T,K
printf("请为矩阵M赋值!\n")
createsmatrix(&M)
printmatrix(M)
printf("请为矩阵N赋值!\n")
createsmatrix(&N)
printmatrix(N)getch()
printf("矩阵M*N得矩阵T!\n")
multsmatrix(M,N,&T)
printmatrix(T)
getch()
printf("对矩阵T进行转置!")getch()
transposesmatrix(T,&K)
printf("转置后的矩阵是\n")
printmatrix(K)
getch()
}
