tfpm_homepage.m

用于模拟时变非平稳的ARMA过程

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The script for the TFARMA Toolbox Tutorial
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;tfpm;
N    = 256;

Aml = [
        0             0.2774 + 0.1809i  -0.0145 + 0.0404i   0.0541 + 0.0448i
        0             0.1762 + 0.1503i   0.0565 - 0.0074i   0.0807 - 0.0163i
   1.0000            -0.2140 - 0.1005i   0.2221 + 0.1832i  -0.1783 - 0.0609i
        0             0.0744 + 0.1909i   0.0529 + 0.1038i   0.1087 - 0.0698i
        0             0.1601 + 0.1545i  -0.0889 - 0.0562i  -0.0116 + 0.0304i
];

Bml = [

  -0.0455 + 0.0222i
   0.0008 - 0.0743i
   0.1617          
   0.0008 + 0.0743i
  -0.0455 - 0.0222i
]

x= tfarma_gen(randn(N, 1), Aml, Bml);
figure(1);clf
subplot(2, 1, 1);plot(real(x), 'Color', 'k', 'Linewidth', 2);axis([1 N -1.3 1.3])
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [ -1 1])
subplot(2, 1, 2);plot(imag(x), 'Color', 'k', 'Linewidth', 2);axis([1 N -1.3 1.3])
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [ -1 1])

[AmlEst, BmlEst]= tfar(x); 
[MEst, LEst]= param_dim(AmlEst)

PxEst= tfarma_evsp(AmlEst, BmlEst, N); 
figure(2);clf;imagesc(rot90(log(PxEst), 1));colormap(flipud(gray));
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])

Psi= tf_multiwin(N, 16, 16);
WxEst= real(ml_to_nk( ambi_est(x, x, -1).*conj(Psi) ));
figure(3);clf;imagesc(rot90(log(WxEst)));colormap(flipud(gray));
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])

Px= tfarma_evsp(Aml, Bml, N);
figure(4);clf;imagesc(rot90(log(Px), 1));colormap(flipud(gray));
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])








%VTFAR

N= 256;
AML(:,:,1,1)= [  0     0
                 0     0];
AML(:,:,2,1)= [  1     0
                 0     1];
AML(:,:,3,1)= [  0     0
                 0     0];
AML(:,:,1,2)= [ -0.0301 + 0.0801i   0.0422 + 0.1571i
                -0.0488 - 0.0455i   0.1999 + 0.0740i];
AML(:,:,2,2)= [  0.2794 - 0.0079i   0.1020 - 0.0698i
                -0.1316 - 0.1935i   0.3989 + 0.2268i];
AML(:,:,3,2)= [ -0.0735 - 0.2695i  -0.2218 + 0.1556i
                 0.2284 - 0.0626i  -0.0291 + 0.1997i];
AML(:,:,1,3)= [ -0.1848 + 0.1263i   0.0352 + 0.2071i
                -0.0592 + 0.0581i   0.0468 + 0.1225i];
AML(:,:,2,3)= [  0.1046 + 0.0451i  -0.0107 + 0.0301i
                 0.0745 + 0.0095i   0.0314 - 0.0086i];
AML(:,:,3,3)= [ -0.0042 + 0.0378i  -0.0433 + 0.1037i
                -0.0274 - 0.2137i  -0.0920 + 0.1374i];
B0L= eye(2);
X= vtfarma_gen(randn(N, 2), AML, B0L);

figure(5);clf
subplot(2, 2, 1);plot(real(X(:, 1)), 'Color', 'k', 'Linewidth', 1);axis([1 N -5 5])
subplot(2, 2, 2);plot(real(X(:, 2)), 'Color', 'k', 'Linewidth', 1);axis([1 N -5 5])
subplot(2, 2, 3);plot(imag(X(:, 1)), 'Color', 'k', 'Linewidth', 1);axis([1 N -5 5])
subplot(2, 2, 4);plot(imag(X(:, 2)), 'Color', 'k', 'Linewidth', 1);axis([1 N -5 5])



[AMLEst, B0LEst]= vtfar(X)
[MEst, LEst]= param_dim(AMLEst)






%load motorlow;x= X(:, 11)+j*.00001*randn(size(X(:, 2)));      
x= sim_nwssus_3;

%x= (randn(100, 1) + j*randn(100, 1))/sqrt(2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% the variable x needs to contain your complex-valued 
% signal of any length
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Into a column vector
x= x(:);

% Mean subtraction
x= x-mean(x);

% Real block length
N0= length(x);

% Histogram plot
figure(1);clf;
subplot(1, 2, 1);hist(real(x));colormap(gray);
subplot(1, 2, 2);hist(imag(x));colormap(gray);

% Embedding in a power-of-two-frame with SNR= 10dB
n= ceil(log2(N0));
N= 2^n;
varre= var(real(x));
varim= var(imag(x));

v1= sqrt(varre/10)*randn(ceil ((N-N0)/2), 1) + j*sqrt(varim/10)*randn(ceil ((N-N0)/2), 1);
v2= sqrt(varre/10)*randn(floor((N-N0)/2), 1) + j*sqrt(varim/10)*randn(floor((N-N0)/2), 1);

y= [v1; x; v2];
maxx= .1*ceil(max(abs(y))*11)

% Embedded signal plot
figure(2);clf;
subplot(2, 1, 1);plot(real(y), 'Color', 'k', 'Linewidth', 2);axis([1 N -maxx maxx]);
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
subplot(2, 1, 2);plot(imag(y), 'Color', 'k', 'Linewidth', 2);axis([1 N -maxx maxx]);
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])

% The FFT length for the spectrogram
NFFT= N/4;

% The window
w= symmpad(hanning(2*ceil(NFFT/6)-1), NFFT);

% Spectrogram
S= specgram([ zeros(N/8, 1); y; zeros(N/8-1, 1)], NFFT, [], w, NFFT-1);
S= [S(N/8+1:end, :); S(1:N/8, :)];
figure(3);imagesc(log(flipud(abs(S).^2)));colormap(flipud(gray));
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [1 N/16 N/8 3*N/16 N/4])
set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])

% Rihaczek distribution
Psi= tf_multiwin(N, 16, 16);
WxEst= real(ml_to_nk( ambi_est(y, y, -1).*conj(Psi) ));
figure(4);clf;imagesc(rot90(log(WxEst)));colormap(flipud(gray));
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])
drawnow

% TFAR Fit
[AmlEst, BmlEst]= tfar(y); 
[MEst, LEst]= param_dim(AmlEst)
PxEst= tfarma_evsp(AmlEst, BmlEst, N); 
figure(5);clf;imagesc(log(rot90(PxEst, 1)));colormap(flipud(gray));
set(gca, 'XTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'XTickLabel', [0 N/4-1 N/2-1 3*N/4-1 N-1])
set(gca, 'YTick', [1 N/4 N/2 3*N/4 N])
set(gca, 'YTickLabel', [N/2-1 N/4-1 0 -N/4 -N/2])