(1) 编程求解，并与用数学软件求解的结果对比。
(2) 考察迭代法的收敛性，若均收敛，对比两种方法的收敛速度。
解：源程序：
①雅克比迭代法：建立函数文件jacobi.m
function [n,x]=jacobi(A,b,X,nm,w)
%用雅克比迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1
m=length(A)
D=diag(diag(A)) %令A=D-L-U,计算矩阵D
L=tril(-A)+D %令A=D-L-U,计算矩阵L
U=triu(-A)+D %令A=D-L-U,计算矩阵U
M=inv(D)*(L+U) %计算迭代矩阵
g=inv(D)*b %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g%用迭代格式进行迭代
if norm(x-X,2)<w
disp('迭代次数为')n
disp('方程组的解为')x
return
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=xn=n+1
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!')
disp('最大迭代次数后的结果为')x
②高斯赛德尔迭代法：建立函数文件gaussseidel.m
function [n,x]=gaussseidel(A,b,X,nm,w)
%用高斯-赛德尔迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1
m=length(A)
I=eye(m) %生成m*m阶的单位矩阵
D=diag(diag(A))%令A=D-L-U,计算矩阵D
L=tril(-A)+D %令A=D-L-U,计算矩阵L
U=triu(-A)+D %令A=D-L-U,计算矩阵U
M=inv(D-L)*U %计算迭代矩阵
g=inv(I-inv(D)*L)*(inv(D)*b) %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g %用迭代格式进行迭代
if norm(x-X,2)<w
disp('迭代次数为')n
disp('方程组的解为')x
return
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=xn=n+1
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!')
disp('最大迭代次数后的结果为')x
运行过程及结果：
①雅克比迭代法的运行过程及结果为：
>>A=[10,3,12,-11,31,3,12]
b=[2-54]
X=[000]
nm=50
w=10^-6
>>jacobi(A,b,X,nm,w)
迭代次数为
n =
14
方程组的解为
x =
0.0254
0.5144
0.2026
②高斯赛德尔迭代法的运行过程及结果为：
>>A=[10,3,12,-11,31,3,12]
b=[2-54]
X=[000]
nm=50
w=10^-6
>>gaussseidel(A,b,X,nm,w)
迭代次数为
n =
6
方程组的解为
x =
0.0254
0.5144
0.2026
③用matlab中的A\b命令的运行过程及结果如下：
>>A=[10,3,12,-11,31,3,12]
b=[2-54]
>>A\b
ans =
0.0254
0.5144
0.2026
(1)由运行结果可知，编程得到的方程组的解为[0.02540.51440.2026]。而用matlab中的A\b命令求出的方程组的解为[0.02540.51440.2026]，完全一致。
(2)由运行过程可知，两种方法均收敛，雅克比迭代次数为14，高斯赛德尔迭代次数为6，说明后者的迭代速度比前者快。
2.用超松弛迭代法求解方程组，分别取松弛因子 ，取迭代初值[000]，迭代30次或满足 时停止计算。
(1) 编程求解，并与用数学软件求解的结果对比。
(2) 考察迭代是否收敛，若收敛，松弛因子取何值时收敛最快，与有关定理的结论对照，看结果是否一致。
解：源程序：
①超松弛迭代法：建立函数文件sor22.m
function [n,x]=sor22(A,b,X,nm,w,ww)
%用超松弛迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1
m=length(A)
D=diag(diag(A)) %令A=D-L-U,计算矩阵D
L=tril(-A)+D %令A=D-L-U,计算矩阵L
U=triu(-A)+D %令A=D-L-U,计算矩阵U
M=inv(D-ww*L)*((1-ww)*D+ww*U) %计算迭代矩阵
g=ww*inv(D-ww*L)*b %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g%用迭代格式进行迭代
if norm(x-X,'inf')<w
disp('迭代次数为')n
disp('方程组的解为')x
return
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=xn=n+1
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!')
disp('最大迭代次数后的结果为')x
②向后超松弛迭代法：建立函数文件sorxh22.m
function [n,x]=sorxh22(A,b,X,nm,w,ww)
%用向后超松弛迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1
m=length(A)
D=diag(diag(A)) %令A=D-L-U,计算矩阵D
L=tril(-A)+D %令A=D-L-U,计算矩阵L
U=triu(-A)+D %令A=D-L-U,计算矩阵U
M=inv(D-ww*U)*((1-ww)*D+ww*L) %计算矩阵迭代矩阵
g=ww*inv(D-ww*U)*b %计算迭代格式中的常数项
%下面是迭代过程
while n<=nm
x=M*X+g%用迭代格式进行迭代
if norm(x-X,'inf')<w
disp('迭代次数为')n
disp('方程组的解为')x
return
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=xn=n+1
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!')
disp('最大迭代次数后的结果为')x
③对称超松弛迭代法：建立函数文件ssor22.m
function [n,x]=ssor22(A,b,X,nm,w,ww)
%用对称超松弛迭代法求解方程组Ax=b
%输入：A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：x为求得的方程组的解构成的列向量，n为迭代次数
n=1
m=length(A)
I=eye(m) %生成m*m阶的单位矩阵
D=diag(diag(A))%令A=D-L-U,计算矩阵D
L=tril(-A)+D %令A=D-L-U,计算矩阵L
U=triu(-A)+D %令A=D-L-U,计算矩阵U
%下面计算迭代矩阵和迭代格式中的常数项
M=inv(D-ww*U)*((1-ww)*D+ww*L)*inv(D-ww*L)*((1-ww)*D+ww*U)g=ww*inv(D-ww*U)*(I+((1-ww)*D+ww*L)*inv(D-ww*L))*b
%下面是迭代过程
while n<=nm
x=M*X+g%用迭代格式进行迭代
if norm(x-X,'inf')<w
disp('迭代次数为')n
disp('方程组的解为')x
return
%上面：达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=xn=n+1
end
%下面：如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp('在最大迭代次数内不收敛!')
disp('最大迭代次数后的结果为')x
%lms算法源程序clear all
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')
1、循环码编码与解码Matlab源程序（实验以（7，4）循环码进行分析）m
=
3
n
=
2^m-1
%定义码长
k
=
n-m
%信息位长
msg
=
randint(k*4,1,2)
%随机提取信号，引起一致地分布的任意整数矩阵
subplot(2,2,1)
stem(msg)
title('编码器输入信号')
p=cyclpoly(n,k)
%循环码生成多项式，n=7，k=4
code
=
encode(msg,n,k,'cyclic',p)
%编码函数，对信号进行差错编码
subplot(2,2,2)
stem(code)
title('编码器输出信号')
recode=decode(code,n,k,'cyclic',p)
%对信号进行译码，对接收到的码字进行译码，恢复
出原始的信息，译码参数和方式必须和编码时采用的严格相同
subplot(2,2,3)
stem(recode)
title('译码器输出信号')
t=-1:0.01:1
x=recode
%将recode赋值给x，并进行长度与fft设定
N=length(x)
fx=fft(x)
df=100/N
n=0:N/2
f=n*df
subplot(2,2,4)
plot(f,abs(fx(n+1))*2/N)
grid
title('频谱图')
2、误码率与信噪比之间的关系程序（以（3，2）循环码进行测试）
m
=
2
n
=
2^m-1
%定义码长
k
=
n-m
%信息位长
Fs=40
%系统采样频率
Fd=1
%码速率
N=Fs/Fd
M=2
%进制数
for
SNRpBit=1:100%信噪比
SNR=SNRpBit/log2(M)
%制造100个信息组，每组k位
msg
=
randint(100,k,[0,1])
code
=
encode(msg,n,k,'cyclic/binary')
%加入噪声
%在已调信号中加入高斯白噪声
noisycode=awgn(code,SNR-10*log10(0.5)-10*log10(N),'measured',[],'dB')
%将浮点数转化为二进制，波形整形过程
for
i=1:100
for
a=1:k+1
if
noisycode(i,a)<0.5
noisycode(i,a)
=
0
else
noisycode(i,a)
=
1
end
end
end
%译码
newmsg
=
decode(noisycode,n,k,'cyclic')
%计算误码率
[number,ratio]=biterr(newmsg,msg)
result(SNRpBit)=ratio
disp(['The
bit
error
rate
is',num2str(ratio)])
end
%不同信噪比下循环码经过加性高斯白噪声信道的误码率
figure(1)
stem(result)
title('循环码在不同信噪比下的误码率')
legend('误码率','*')
xlabel('信噪比')
ylabel('在加性高斯白噪声下的误码率')
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）