气象要素奇异谱分析

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

function [signalFiltered]=SSA(signal,windowLen)

%========================================================================

% 参看http://www.ilovematlab.cn/thread-301868-1-1.html柳絮艳的分享代码ssa改写

% signal 原始信号

% windowLen 窗口长度

% signalFiltered   重构时间序列

%=========================================================================

% Step1 : 建立轨迹矩阵

N=length(signal);

if windowLen>N/2;

    windowLen=N-windowLen;

end

K=N-windowLen+1;

X=zeros(windowLen,K);

for i=1:K

    X(1:windowLen,i)=signal(i:windowLen+i-1);

end

% Step 2: 奇异值分解

S=X*X';

[U,autoval]=eig(S);%eig返回矩阵的特征值和特征向量，U是特征向量，autoval是特征值

[d,i]=sort(diag(autoval),'descend');

sev=sum(d); %特征值的和,可根据占比选择有效信号

U=U(:,i);

V=(X')*U;

% Step 3:分组

I=[1:windowLen/2];%I的选择可根据信号特征选择

Vt=V';

rca=U(:,I)*Vt(I,:);

% Step 4: 对交平均化重构信号

y=zeros(N,1);

Lp=min(windowLen,K);

Kp=max(windowLen,K);

%重构 1~Lp-1

for k=0:Lp-2

    for m=1:k+1;

        y(k+1)=y(k+1)+(1/(k+1))*rca(m,k-m+2);

    end

end

%重构 Lp~Kp

for k=Lp-1:Kp-1

    for m=1:Lp;

        y(k+1)=y(k+1)+(1/(Lp))*rca(m,k-m+2);

    end

end

%重构 Kp+1~N

for k=Kp:N-1

    for m=k-Kp+2:N-Kp+1;

        y(k+1)=y(k+1)+(1/(N-k))*rca(m,k-m+2);

    end

end

signalFiltered = y;

end