main1.m

arma模型的功率谱估计代码。还包括系统的极零点分析等完整软件包。

%main for 9-(1)
close all;
clear;

dataDrop = 200;     %丢弃点数,200
dataLength = 256;   %数据长度,256
plotNum = 128;      %作图取点数
times = 20;         %独立实验数

reala = [1, -1.3817, 1.5632, -0.8843, 0.4096];
realb = [1, 0.3544, 0.3508, 0.1736, 0.2401];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%AR(4)
for t=1:times
    x = xx(reala, realb, dataDrop, dataLength);
    [coa, coe] = AR(x, 4);
    coa = [1;coa];
    cob = sqrt(coe);
    [w(:,t),P(:,t)] = spectrum(coa, cob, plotNum);
end
%20次实验结果
figure(1)
title('AR(4)算法20次实验结果-功率谱');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
for t=1:times
    plot(w(:,t), 10*log(P(:,t)));
end
%平均结果和真实谱
figure(2)
title('AR(4)算法对数后平均(兰)、平均后对数(绿)和真实谱(红)');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
meanw = mean(w, 2);
meanP = mean(P, 2);
logP = 10*log(P);
meanlogP = mean(logP, 2);
[realw, realP] = spectrum(reala, realb, plotNum);
plot(meanw, meanlogP, 'b');
plot(meanw, 10*log(meanP), 'g');
plot(realw, 10*log(realP), 'r');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%AR(8)
for t=1:times
    x = xx(reala, realb, dataDrop, dataLength);
    [coa, coe] = AR(x, 8);
    coa = [1;coa];
    cob = sqrt(coe);
    [w(:,t),P(:,t)] = spectrum(coa, cob, plotNum);
end
%20次实验结果
figure(3)
title('AR(8)算法20次实验结果-功率谱');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
for t=1:times
    plot(w(:,t), 10*log(P(:,t)));
end
%平均结果和真实谱
figure(4)
title('AR(8)算法对数后平均(兰)、平均后对数(绿)和真实谱(红)');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
meanw = mean(w, 2);
meanP = mean(P, 2);
[realw, realP] = spectrum(reala, realb, plotNum);
plot(meanw, meanlogP, 'b');
plot(meanw, 10*log(meanP), 'g');
plot(realw, 10*log(realP), 'r');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ARMA(4,4)
for t=1:times
    x = xx(reala, realb, dataDrop, dataLength);
    [coa, cob, coe] = ARMA(x, 4, 4);
    cob = sqrt(coe)*cob;
    [w(:,t),P(:,t)] = spectrum(coa, cob, plotNum);
end
%20次实验结果
figure(5)
title('ARMA(4,4)算法20次实验结果-功率谱');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
for t=1:times
    plot(w(:,t), 10*log(P(:,t)));
end
%平均结果和真实谱
figure(6)
title('ARMA(4,4)算法对数后平均(兰)、平均后对数(绿)和真实谱(红)');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
meanw = mean(w, 2);
meanP = mean(P, 2);
logP = 10*log(P);
meanlogP = mean(logP, 2);
[realw, realP] = spectrum(reala, realb, plotNum);
plot(meanw, meanlogP, 'b');
plot(meanw, 10*log(meanP), 'g');
plot(realw, 10*log(realP), 'r');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ARMA(8,8)
for t=1:times
    x = xx(reala, realb, dataDrop, dataLength);
    [coa, cob, coe] = ARMA(x, 8, 8);
    cob = sqrt(coe)*cob;
    [w(:,t),P(:,t)] = spectrum(coa, cob, plotNum);
end
%20次实验结果
figure(7)
title('ARMA(8,8)算法20次实验结果-功率谱');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
for t=1:times
    plot(w(:,t), 10*log(P(:,t)));
end
%平均结果和真实谱
figure(8)
title('ARMA(8,8)算法对数后平均(兰)、平均后对数(绿)和真实谱(红)');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
meanw = mean(w, 2);
meanP = mean(P, 2);
logP = 10*log(P);
meanlogP = mean(logP, 2);
[realw, realP] = spectrum(reala, realb, plotNum);
plot(meanw, meanlogP, 'b');
plot(meanw, 10*log(meanP), 'g');
plot(realw, 10*log(realP), 'r');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Periodic
for t=1:times
    x = xx(reala, realb, dataDrop, dataLength);
    [w(:,t),P(:,t)] = Periodic(x, plotNum);
end
%20次实验结果
figure(9)
title('周期图算法20次实验结果-功率谱');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
for t=1:times
    plot(w(:,t), 10*log(P(:,t)));
end
%平均结果和真实谱
figure(10)
title('周期图算法对数后平均(兰)、平均后对数(绿)和真实谱(红)');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
meanw = mean(w, 2);
meanP = mean(P, 2);
logP = 10*log(P);
meanlogP = mean(logP, 2);
[realw, realP] = spectrum(reala, realb, plotNum);
plot(meanw, meanlogP, 'b');
plot(meanw, 10*log(meanP), 'g');
plot(realw, 10*log(realP), 'r');