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');