第一周信号处理实验课
上课学习
学校东区计算机学院708机房没有安装好MATLAB，而且有些机位不能联网，上课简单介绍了一下MATLAB和课后作业。
下次上课带自己的电脑去。
课后作业
一些MATLAB的简单操作，课程报告第二周后再提交。
回去在自己电脑上安装好MATLAB，学校网速太慢了，而且还不稳定，百度云限速。
MATLAB安装办法自己网上找。
MATLAB学习
百度云下了我两天两夜，终于下载安装好了。
明天开始学。
学了半天，看了两节课，还是不会做作业……
第五周信号处理实验课
Command conv
First rule: The two signals (input and impulse response) should be defined in the same time interval.
Second rule: When a signal consists of multiple parts, the time intervals in which each part is defined must not overlap.
Third rule: The output of the conv command has to be multiplied with the time step used in the definition of the input and impulse response signals, in order to correctly compute the output of the system.
Fourth rule: The output of the system is plotted in the double time interval of the one in which the input and impulse response signals are defined.
%1st way
t1=0:.01:1-0.01;
t2=1:.01:2;
t3=2.01:.01:10;
t=[t1 t2 t3];
x1=zeros(size(t1));
x2=ones(size(t2));
x3=zeros(size(t3));
x=[x1 x2 x3];
h=x;
y=conv(x,h)*.01;
plot(0:.01:20,y);
axis([0 6 -.1 1.1])
legend('y(t)')
%the two signals
figure
plot(t,x,0-t,h,':')
legend('x(t)','h(t-\tau)=h(0-\tau)')
ylim([-.1 1.1])
figure
plot(t,x,2-t,h,':')
ylim([-.1 1.1])
legend('x(t)','h(t-\tau)=h(2-\tau)')
%2nd way
figure
t=1:.01:2;
x=ones(size(t));
h=ones(size(t));
y=conv(x,h)*.01;
plot(2:.01:4,y)
legend('y(t)')
Other command
% syms
% Create symbolic variables and functions
% 符号变量
hold on是当前轴及图像保持而不被刷新，准备接受此后将绘制的图形，多图共存。即，启动图形保持功能，当前坐标轴和图形都将保持，从此绘制的图形都将添加在这个图形的基础上，并自动调整坐标轴的范围。
hold off使当前轴及图像不再具备被刷新的性质，新图出现时，取消原图。即，关闭图形保持功能。
area(x, y)：该函数以参数x和y绘制面积图。如果x和y为向量，则相当于函数plot(x, y)，并将0到y之间进行了填充。如果参数y为矩阵，则将y的每一列绘制面积图并进行累计求和。
area(y)：如果参数y为向量，则绘制面积图；如果y为矩阵则绘制每一列的面积图之和。
simplify的调用格式为：simplify(S)，对表达式S进行化简
矩阵代数上机实验
第一周学习matlab
MATLAB 使用撇号运算符 (') 执行复共轭转置，使用点撇号运算符 (.') 执行无共轭的转置。对于包含所有实数元素的矩阵，这两个运算符返回相同结果。
非共轭复数转置（其中每个元素的复数部分保留其符号）表示为 z.'：
eye(n) 返回 n×n 单位方阵
inv 函数和表达式 A^-1 均可对矩阵求逆
有些矩阵接近奇异矩阵，虽然存在逆矩阵，但计算容易出现数值误差。cond 函数计算逆运算的条件数，它指示矩阵求逆结果的精度。条件数的范围是从 1（数值稳定的矩阵）到 Inf（奇异矩阵）。
很少需要为某个矩阵构造显式逆矩阵。求解线性方程组 Ax = b 时，常常会误用 inv。从执行时间和数值精度方面而言，求解此方程的最佳方法是使用矩阵反斜杠运算符，即 x = A\b。有关详细信息，请参阅 mldivide。
Kronecker 张量积
两个矩阵的 Kronecker 乘积 kron(X,Y) 为 X 的元素与 Y 的元素的所有可能乘积构成的较大矩阵。如果 X 为 m×n 且 Y 为 p×q，则 kron(X,Y) 为 mp×nq。元素以特定方式排列，呈现 X 的每个元素分别与整个矩阵 Y 相乘的结果。
向量范数和矩阵范数
向量 x 的 p-范数，
使用 norm(x,p) 进行计算。此运算是为 p > 1 的任意值定义的，但最常见的 p 值为 1、2 和 ∞。默认值为 p = 2，这与欧几里德长度或向量幅值对应：
v = [2 0 -1];
[norm(v,1) norm(v) norm(v,inf)]
ans =
3.0000 2.2361 2.0000矩阵 A 的 p-范数，
可以针对 p = 1、2 和 ∞ 通过 norm(A,p) 进行计算。同样，默认值也为 p = 2：
A = pascal(3);
[norm(A,1) norm(A) norm(A,inf)]
ans =
10.0000 7.8730 10.0000如果要计算矩阵的每行或每列的范数，可以使用 vecnorm：
vecnorm(A)
ans =
1.7321 3.7417 6.7823第二周InverseMatrix
function B= InverseMatrix (A)
[m,n] = size(A);
if m~=n
disp('A is not a square matrix');
return
end
if rank(A) < n
disp('A is not a full rank matrix');
return
end
C = [A ;eye(n)];
%初等列变换代码
for i=1:n
%第一步对角线元素变1
%第二步变A为单位矩阵
d_1=C(i,i);
C(:,i)=C(:,i)/d_1;
for j=1:n
if i~=j
d_2=C(i,j);
C(:,j)=C(:,j)-d_2.*C(:,i);
end
end
end
B = C(n+1:end,:);
end忘了除数不为0
第三周InverseMatrixCramer
n=5;
A=rand(n,n);
B=InverseMatrixCramer(A)
inv(A)
function B = InverseMatrixCramer (A)
[m,n] = size(A);
if m~=n
disp('A is not a square matrix');
return
end
if rank(A) < n
disp('A is not a full rank matrix');
return
end
B = zeros(size(A));
%Cramer's rule
C=eye(n);
for i=1:n
c=C(:,i);
B_det=det(A);
for j=1:n
D=A;
D(:,j)=c;
B_det_j=det(D);
b=B_det_j/B_det;
B(i,j)=b;
end
end
B=B';
end第四周SloveEquationKron
%% 可观察四个变量不同取值的变化
m=2; p=4; q=4; n=2;
A = rand(m,p); X= zeros(p,q);
B = rand(q,n); D = rand(m,n);
X = SloveEquationKron(A,B,D);
%%% 验证程序是否正确
norm(kron(B',A)\D(:)-X(:),'fro') %正确输出0
norm(A*X*B-D,'fro')
function X = SloveEquationKron(A,B,D)
[m,p] = size(A);
[q,n] = size(B);
[m,n] = size(D);
%K=kron(B',A)
B_=B';
[ma,na] = size(B_);
[mb,nb] = size(A);
[ia,ib] = meshgrid(1:ma,1:mb);
[ja,jb] = meshgrid(1:na,1:nb);
K = B_(ia,ja).*A(ib,jb);
%Solve equation via kron
C=K;
XX = C\D(:);
X = reshape(XX,[p,q]);
end
第五周SVD_my
function [U,D,V] = SVD_my(A)
%[m,n] = size(A);
%U=rand(m,m);
%V=rand(n,n);
%D=A;
T1 = null(A');%计算A共轭转置零空间的标准正交基
B=A'*A;
[V,D]=eig(B);%求特征值
[m0,n0]=size(B);
for i=1:m0
for j=1:n0
if i+j==m0+1
E(i,j)=1;
else
E(i,j)=0;
end
end
end
V=V*E;
D=E*D*E;%
V1=V;
T=diag(D);
T0=T;
D0=D;
[m1,n1]=size(T);
for i=m1:-1:1
if abs(T0(i)-0)< 10^(-10)
T(i)=[];
V1(:,i)=[];
D0(i,:)=[];
end
end
S=sqrt(D0);
D1=diag(sqrt(T));
U=A*V1*inv(D1);
[m,n] = size(A);
[md,nd]=size(D);
z=zeros(m-md,n);
D=[S;z];
U=[U T1];
end第七周满秩分解求广义逆
close all
clear all
%m=10; n=8;
%A = rand(m,n);
A=[1 2 0;0 1 1;1 3 1];
[A_MP,L,R] = fullrank_my(A)
%%% 验证程序是否正确
norm(A_MP-pinv(A),'fro')%%正确为接近于0
%%
function [A_MP,L,R] = fullrank_my(A)
A1 = rref(A); %将矩阵A化成行最简形式保存在A1中
[m,n]=size(A); %获取矩阵A的大小：m行n列
B0= [];%生成一个空向量
C0= [];%生成一个空向量
for i=1:m %依次扫描矩阵m行
flag=1;
for j=1:n %依次扫描矩阵n列
if A1(i,j)==1 %若A1(i, j)等于1
for k=1:i-1 %固定j列，扫描此列的第1行到i-1行元素
if A1(k,j)~=0 %判断是否全为0
flag=0; %若不全为0，则将flag置为0（说明此列不是单位矩阵的列）
break;
end
end
for k=i+1:m %固定j列，扫描此列的第i+1行到m行（即最后一行）元素
if A1(k,j)~=0 %判断是否全为0
flag=0; %若不全为0，则将flag置为0（说明此列不是单位矩阵的列）
break;
end
end
if flag==1 %若flag为1（不为0），则说明此列是【矩阵的行最简形式矩阵】的单位矩阵的列
B0=[B0,A(:,j)]; %将矩阵A的j列加到B0列向量之后
C0=[C0;A1(i,:)]; %将矩阵A1的i行加到C0行向量之后，
end
end
end
end
[m1,n1]=size(B0); %获取矩阵B0的大小：m1行n1列
[m2,n2]=size(C0); %获取矩阵C0的大小：m2行n2列
B=B0(:,1:n1); %将矩阵B0的第1列到最后一列赋值给矩阵B
C=C0(1:m2,:); %将矩阵C0的第1行到最后一行赋值给矩阵C
L=B;
R=C;
%A_MP=R'*(R*R')^(-1)*(L'*L)^(-1)*L';
A_MP=R'*inv(R*R')*inv(L'*L)*L';
end
