基于松弛思想的RELAX算法——matlab代码(附详细注释)

% relax.m
% Author:迷茫打工人
% Date：2021/03/16

clc;
clear;
close all;
%% 设定初始参数，产生原始信号
Num=512;        % 512个采样点
Num_fft=1024;   % 补0至1024个采样点
Fs=2*pi;        % 采样频率为2pi
n=0:Num-1;  
xn=2*cos(0.2*pi*n)+5*cos(0.4*pi*n)+randn(size(n));    % 原始时域信号
% randn(size(n))返回一个和n有同样维数大小的符合正态分布的随机数组

%% 原始时域信号波形
figure(1);      % 第一幅图
subplot(2,1,1); % 多图形绘制，将图形排成2行1列
plot(n,xn);     % 绘制信号的X坐标和Y坐标
title('xn=2*cos(0.2*pi*n)+5*cos(0.4*pi*n)+randn(size(n))'); 
axis([0 Num,-10 10]);   % 限制X轴范围0至Num，Y轴范围-10至10
grid on;    % 加入网格线

%% RELAX算法
KK=4; 
y(1,:)=xn;      % 将时域信号xn赋值给y矩阵的第一行
yyk=fftshift(fft(y(1,:),Num_fft));     % 将原始信号补0后进行FFT变换
                                       % fftshift将零频点移到频谱的中间  
                                       % 需要手动将频率序列向左移fs/2                                 
%%%%%%%%%%%%%%% 频率计算公式为：fk=argmax arg||DFT(yk)||^2 %%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% 复包络计算公式为：ak=w'(fk)*yk/N %%%%%%%%%%%%%%%
[ywk(1),m]=max(abs(yyk));    % 将频谱中最大幅值赋给ywk(1),对应的位置索引赋给m
w(1)=-pi+(m-1)*2*pi/Num_fft;
b(1)=y(1,:)*(exp(-j*w(1)*(0:Num-1)'))/Num;  % 将该正弦分量的复包络赋给b(1)
b1=b;
w1=w;  
g(1)=0; 

%% 当信号中的正弦分量的数量>=2
for K=2:KK  
    cast=1;     % 迭代阈值变量初始化
    c=K-1;      % 第一次循环，此时c=1
    while (cast>=0.0001)   % 当变化量小于0.0001时，收敛结束
        k=mod(c,K)+1;      % 取余运算，c被K除过后的余数+1，后面用于选取分量进行迭代
        % 将(K-1)个正弦分量累和到sum()中
        for ii=1:Num 
           sum(ii)=0;    % 清零
           for i=1:K-1
              sum(ii)=sum(ii)+b1(i)*exp(j*w1(i)*(ii-1));
           end 
        end 
        y(k,:)=xn-sum; %将原始信号xn去除(K-1)个正弦分量后，留下一分量赋给y矩阵的第k行
%% 开始计算该分量的频率及复包络
        yyk=fftshift(fft(y(k,:),Num_fft));  % 剩余分量补0进行FFT变换
        yyk2=(abs(yyk))/Num; 
        [ywk(K),m]=max(yyk2);   % 将频谱的最大幅值存放ywk()，对应的位置索引赋给m                       
        w(k)=-pi+(m-1)*2*pi/Num_fft;    % 该剩余分量的的频率值
        b(k)=y(k,:)*(exp(-j*w(k)*(0:Num-1)'))/Num; 
        f_b=b;
        f_w=w;
        % 选取分量组，用一组分量计算另一个分量
        if k==K
           w1=f_w(2:K); 
           b1=f_b(2:K); 
        end 
        if k==K-1 
           w1=f_w(1:K-1); 
           b1=f_b(1:K-1); 
        end 
        if (k~=K-1)&&(k~=K)   
           w1=[f_w(1:k) f_w(k+2:K)]; 
           b1=[f_b(1:k) f_b(k+2:K)]; 
        end 
        % 将K个正弦分量累和到be()中
        for ii=1:Num 
           be(ii)=0;    % 清零
           for i=1:K 
              be(ii)=b(i)*exp(j*w(i)*ii)+be(ii);
           end 
        end  
        g(c+1)=(xn-be)*(xn-be)';    % (xn-be)表示原始信号除去所有计算出的分量，剩下噪声及其他数据
        cast=abs(g(c+1)-g(c));      % 计算本次迭代变化量
        c=c+1;
    end 
    w1=w;   % 迭代收敛，w1存储各正弦分量的频率信息
    b1=b;   % 迭代收敛，b1存储各正弦分量的幅值+相位信息
end 
RelaxPSD=10*log(2*abs(b));      % 幅值转换为dB值
%% 得到的功率谱图
subplot(2,1,2);   % 多图形绘制，该图在自上而下的第二个位置
w1=w1/pi;    % w1存储各分量频率值，为方便表示将频率值除以pi
stem(w1,RelaxPSD); % 在w1的指定点处将数据序列RelaxPSD从x轴到数据值按茎状形式画出，以圆圈终止
xlabel('Normalized Frequency(*pi rad/sample)');
ylabel('power spectral(dB)');
title('RELAX power spectral');
grid on;