main2.m

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

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

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

reala = [1, -1.6408, 2.2044, -1.4808, 0.8145];
realb = [1, 1.5857, 0.9604];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%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);
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(12)
for t=1:times
    x = xx(reala, realb, dataDrop, dataLength);
    [coa, coe] = AR(x, 12);
    coa = [1;coa];
    cob = sqrt(coe);
    [w(:,t),P(:,t)] = spectrum(coa, cob, plotNum);
end
%20次实验结果
figure(5)
title('AR(12)算法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('AR(12)算法对数后平均(兰)、平均后对数(绿)和真实谱(红)');
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(16)
for t=1:times
    x = xx(reala, realb, dataDrop, dataLength);
    [coa, coe] = AR(x, 16);
    coa = [1;coa];
    cob = sqrt(coe);
    [w(:,t),P(:,t)] = spectrum(coa, cob, plotNum);
end
%20次实验结果
figure(7)
title('AR(16)算法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('AR(16)算法对数后平均(兰)、平均后对数(绿)和真实谱(红)');
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(4,2)
for t=1:times
    x = xx(reala, realb, dataDrop, dataLength);
    [coa, cob, coe] = ARMA(x, 4, 2);
    cob = sqrt(coe)*cob;
    [w(:,t),P(:,t)] = spectrum(coa, cob, plotNum);
end
%20次实验结果
figure(9)
title('ARMA(4,2)算法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('ARMA(4,2)算法对数后平均(兰)、平均后对数(绿)和真实谱(红)');
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,4)
for t=1:times
    x = xx(reala, realb, dataDrop, dataLength);
    [coa, cob, coe] = ARMA(x, 8, 4);
    cob = sqrt(coe)*cob;
    [w(:,t),P(:,t)] = spectrum(coa, cob, plotNum);
end
%20次实验结果
figure(11)
title('ARMA(8,4)算法20次实验结果-功率谱');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
for t=1:times
    plot(w(:,t), 10*log(P(:,t)));
end
%平均结果和真实谱
figure(12)
title('ARMA(8,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(12,6)
for t=1:times
    x = xx(reala, realb, dataDrop, dataLength);
    [coa, cob, coe] = ARMA(x, 12, 6);
    cob = sqrt(coe)*cob;
    [w(:,t),P(:,t)] = spectrum(coa, cob, plotNum);
end
%20次实验结果
figure(13)
title('ARMA(12,6)算法20次实验结果-功率谱');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
for t=1:times
    plot(w(:,t), 10*log(P(:,t)));
end
%平均结果和真实谱
figure(14)
title('ARMA(12,6)算法对数后平均(兰)、平均后对数(绿)和真实谱(红)');
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(15)
title('周期图算法20次实验结果-功率谱');
xlabel('角频率w (rad)');
ylabel('功率谱密度P-对数 (dB)');
grid;
hold on;
for t=1:times
    plot(w(:,t), 10*log(P(:,t)));
end
%平均结果和真实谱
figure(16)
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');