svd.m

谐波恢复的ARMA建模算法(基于SVD-TLS)（张贤达《现代信号处理》）。当谐波信号为复谐波的时

n = 256;         %采样个数
m = 1:n;
w = randn(1,n);
w = w/std(w);        %保证为零均值，方差为1
w = w - mean(w);
x(m) = 20^0.5*sin(2*pi*0.2*m)+2^0.5*sin(2*pi*0.213*m);   %每个余弦波的信噪比(SNR)定义为该余弦波的功率与噪声平均功率即方差之比
                                                         %而每个余弦拨的功率等于该余弦波幅值有效值的平方
%figure(1);plot(x);
%figure(2);plot(w)
y(n) = x(n) + w(n);
%figure(3);plot(y);
[c,lags] = xcorr(y);   % c里从 R(-n+1)，…,R(0),…，R(n-1) 一共(2n-1)个值
b = zeros(1,n-1);
b = c(1,(n+1):(2*n-1));     % 取R(1)到R(n-1)一共n-1个值 
M = 8;Pe = 32;
R = zeros(M,Pe+1);          % M = 6，Pe = 13,M×(Pe+1)的矩阵
for f = 1:Pe+1;                %分别给R的各列赋值，Page110
    R(:,f) = b(1,Pe+1+1-f:Pe+M+1-f);
end

                                       % R(:,1) = c(1,14:19);      %分别给R的各列赋值
                                       % R(:,2) = c(1,13:18);
                                       % R(:,3) = c(1,12:17);
                                       % R(:,4) = c(1,11:16);
                                       % R(:,5) = c(1,10:15);
                                       % R(:,6) = c(1,9:14);
                                       % R(:,7) = c(1,8:13);
                                       % R(:,8) = c(1,7:12);
                                       % R(:,9) = c(1,6:11);
                                       % R(:,10) = c(1,5:10);
                                       % R(:,11) = c(1,4:9);
                                       % R(:,12) = c(1,3:8);
                                       % R(:,13) = c(1,2:7);
                                       % R(:,14) = c(1,1:6);
                                       
[U,S,V] = svd(R);            %求R的奇异值,同时保存求奇异值所需的两个酉矩阵
A = diag(S);L = size(A);     %保存所有奇异值,并保存奇异值的个数
AF = A'*A;                   %得到所有奇异值的平方和，即S的F范数
for i = 1:L;
B(1:i,:) = A(1:i,:)         %取A的前i个元素 
BF = B'*B;                  
  if (BF/AF > 0.995)
      p = i                 %增广矩阵的有效秩
   break
  end
end
for j = 1:p;                   %式(3.4.56)
    for k = 1:(Pe+1-p);          %n = 13,n+1 = 14
        sum1 = zeros(p+1);      %定义一个(p+1)*(p+1)的矩阵
        sum2 = zeros(p+1);
        sum2 = A(j)^2*V(k:k+p,j)*V(k:k+p,j)';        
        sum1 = sum1 + sum2;     %式(3.4.56)中第一个连加号内的和
    end
   sum = zeros(p+1);
   sum = sum + sum1;            %式(3.4.56)中第二个连加号内的和，即得Sp
end
Spn = inv(sum);                  %Sp的逆
X = Spn(:,1)/Spn(1,1)           
a = [1 X(2:p+1,1)' X(p:-1:2,1)' 1]    %未知参数的最小二乘解,αi = α2p-i 式(3.7.4),Page116
Z = roots(a)
F1 = abs(atan(imag(Z(1))/real(Z(1))))/2/pi
F2 = abs(atan(imag(Z(2))/real(Z(2))))/2/pi
F3 = abs(atan(imag(Z(3))/real(Z(3))))/2/pi
F4 = abs(atan(imag(Z(4))/real(Z(4))))/2/pi
F5 = abs(atan(imag(Z(5))/real(Z(5))))/2/pi
F6 = abs(atan(imag(Z(6))/real(Z(6))))/2/pi
F7 = abs(atan(imag(Z(7))/real(Z(7))))/2/pi
F8 = abs(atan(imag(Z(8))/real(Z(8))))/2/pi
F = [F1 F2 F3 F4 F4 F5 F6 F7 F8]
clear