信号谱分析matlab实验

clear all ;
close all ;
clc;
cd F:\MATLAB\spectralanalysis ;%file path of experiment_data.mat; load experiment_data ;
figure;
plot(data);
N1=64;%data length;
N2=256;
data64=[data(65:65+N1-1) data(129:129+N1-1) data(193:193+N1-1)];%the three statistical samples of data using starting points of 64,128,and 192 for N=64;
data256=[data(257:257+N2-1) data(513:513+N2-1)
data(769:769+N2-1)];%the three statistical samples of data using starting points of 256,512,and 768 for N=256;
Ts=1;%sampling interval is 1s;
%1.Periodogram Methods
%9-6a
K=8;
figure(1);
for time1=1:3
plot(linspace(0,1/Ts,K*N1),10*log10(abs(fft(data64(:,time1),K*N1)).^2/length(data64(:,time1))));
hold on ;
end
hold off ;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-6a periodogram,N=64,rect window,no smooth');
%9-6(b)
我们对这个模型采样，获得1024个采样点，结合极限谱密度的图，我们可以知道我们的数据中包含三个sine 波形以及附加的高斯有色噪声（colored Gaussian noise ）
wn1=2*hanning(N1);
figure(2);
for time1=1:3
datatapering64(:,time1)=data64(:,time1).*wn1;
plot(linspace(0,1/Ts,K*N1),10*log10(abs(fft(datatapering64(:,time1),K *N1)).^2/length(datatapering64(:,time1))));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-6b periodogram,N=64,hanning window,no smooth');
%--------------------------------------------------------------------%9-6c
figure(3);
for time1=1:3
Sf=(abs(fft(data64(:,time1),K*N1)).^2/length(data64(:,time1)));
for n=1:(length(Sf)-16)
Sfmat(n,:)=Sf(n:n+15,:);
end
Sfsmooth=Sfmat*ones(16,1)./16;
plot(linspace(0,1/Ts,n),10*log10(Sfsmooth));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-6c periodogram,N=64,rect window,smooth M=2');
%9-6d
wn1=2*hanning(N1);
figure(4);
for time1=1:3
datatapering64(:,time1)=data64(:,time1).*wn1;
Sf=(abs(fft(datatapering64(:,time1),K*N1)).^2/length(datatapering64(: ,time1)));
for n=1:(length(Sf)-16)
Sfmat(n,:)=Sf(n:n+15,:);
end
Sfsmooth=Sfmat*ones(16,1)./16;
plot(linspace(0,1/Ts,n),10*log10(Sfsmooth));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-6d periodogram,N=64,hanning window,smooth M=2');
从图9-6a中的周期图中可以看出，频谱假峰很多，可靠性差、分辨力低，难以从中得到正确的频谱信息。

其中的补零处理，只是为了减小做DFT时带来的栅栏效应，对频谱的分辨力并没有帮助。

加升余弦窗处理后的频谱（图9-6b）情况好了些，这是由于矩形窗（no data tapering）的旁瓣电平起伏大，卷积运算会产生较大的失真，而升余弦窗的频谱旁瓣显著减小，由此截断得到的周期图的泄露明显减少。

对比9-6c、9-6d的频率平滑处理，从频谱结果可以看到相较于没有平滑，平滑后的图形起伏明显平坦多了，假峰明显减少，可靠性得到了提升。

%--------------------------------------------------------------------%9-7a
figure(5)
for time2=1:3
plot(linspace(0,1/Ts,K*N2),10*log10(abs(fft(data256(:,time2),K*N2)).^ 2/length(data256(:,time2))));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-7a periodogram,N=256,rect window,no smooth');
%9-7c
figure(7);
for time2=1:3
Sf=abs(fft(data256(:,time2),K*N2)).^2/length(data256(:,time2));
for n=1:(length(Sf)-16)
Sfmat(n,:)=Sf(n:n+15,:);
end
Sfsmooth=Sfmat*ones(16,1)./16;
plot(linspace(0,1/Ts,n),10*log10(Sfsmooth));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-7c periodogram,N=256,rect window,smooth M=2');
%--------------------------------------------------------------------% %9-7b
wn2=2*hanning(N2);
figure(6);
for time2=1:3
datatapering256(:,time2)=data256(:,time2).*wn2;
plot(linspace(0,1/Ts,K*N2),10*log10(abs(fft(datatapering256(:,time2), K*N2)).^2/length(datatapering256(:,time2))));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-7b periodogram,N=256,hanning window,no smooth');
%9-7d
wn2=2*hanning(N2);
figure(8);
for time2=1:3
datatapering256(:,time2)=data256(:,time2).*wn2;
Sf=abs(fft(datatapering256(:,time2),K*N2)).^2/length(datatapering256(
:,time2));
for n=1:(length(Sf)-16)
Sfmat(n,:)=Sf(n:n+15,:);
end
Sfsmooth=Sfmat*ones(16,1)./16;
plot(linspace(0,1/Ts,n),10*log10(Sfsmooth));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-7d periodogram,N=256,hanning window,smooth M=2')
从上面四个不同处理的周期图中也可以看到，作升余弦窗处理后的周期图的泄露明显减少。

频率平滑处理提高了可靠性。

在处理数据点数N=64和N=256纵向比较来看，利用的数据越多，谱分析的结果可靠性越高，这符合常理。

%2.sine-wave removal
clear;clc;clf;
%to load data
load('experiment_data');
% data_to_pro=[data(65:128) data(129:192) data(193:256)];t=64;%chosen
data length!!!!!
data_to_pro=[data(257:512) data(513:768)
data(769:1024)];t=256; %chosen data length!!!!!
N=8;%padding number
for index=1:3;%to control which segment to use
freq=(0:1/(N*t-1):1)';
pxx=10*log10((abs(fft(data_to_pro(:,index),t*N))).^2./t);
figure(2)
subplot(221)
plot(freq,pxx);hold on;
axis([0,0.5,-20,20]);grid;
title('PSD of the original data');
n=(1:t)';
%-----------------the first estimate--------------
[a,i]=max(pxx);%estimated freq using formula 234
pulse=zeros(length(pxx),1);pulse(i)=a;istore=ones(length(pxx),1);isto re(i)=0;%to record pulse in order to reinsert afterwards
a=10^(a/10);%estimated a is not in dB
fcap=i/(t*N)*2*pi;
sincap=sin(fcap.*n-(t-1)/2*fcap);%estimated sin wave in formula 233 coscap=cos(fcap.*n-(t-1)/2*fcap);%estimated cos wave in formula 233 pusi=-atan((sincap'*data_to_pro(:,index))/(coscap'*data_to_pro(:,inde x)))/pi*180;%estimate phase
acap=2*sqrt(a/t);%formula 235
estimate=acap.*cos(fcap.*n+pusi/180*pi-(t-1)/2*fcap);% this is estimated cosine wave use a, f, pusi
%--------------the first removal--------------------
removal1_add=data_to_pro(:,index)+estimate;%since the phase estimate will probably introduce a 'pi' error,whether add or substract the estimate sine wave is unclear
removal1_sub=data_to_pro(:,index)-estimate;
if(removal1_add'*removal1_add>removal1_sub'*removal1_sub)%correct removal will decrease the power of vector 'data' %
removal1=removal1_sub;%therefore after compare the two possible result we can decide which one to chose
else
removal1=removal1_add;
end
figure(2)
subplot(222)
pxx=10*log10((abs(fft(removal1,t*N))).^2./t);
plot(freq,pxx.*istore+pulse);axis([0,0.5,-20,20]);hold on;grid on;
title('PSD after the first removal')
%-----------the second estimate----------------
[a,i]=max(pxx);
pulse(i)=a;istore(i)=0;
a=10^(a/10);
fcap=i/(t*N)*2*pi;
sincap=sin(fcap.*n-(t-1)/2*fcap);
coscap=cos(fcap.*n-(t-1)/2*fcap);
pusi=-atan((sincap'*removal1)/(coscap'*removal1))/pi*180;
acap=2*sqrt(a/t);
estimate=-acap.*cos(fcap.*n+pusi/180*pi-(t-1)/2*fcap);
%------------the second removal-------------------------
removal2_add=removal1+estimate;
removal2_sub=removal1-estimate;
if(removal2_add'*removal2_add>removal2_sub'*removal2_sub)
removal2=removal2_sub;
else
removal2=removal2_add;
end
figure(2)
subplot(223)
pxx=10*log10((abs(fft(removal2,t*N))).^2./t);
plot(freq,pxx.*istore+pulse);axis([0,0.5,-20,20]);hold on;grid on; title('PSD after second removal')
%-----------------the third estimate--------around 0.1fs----------[a,i]=max(pxx(1:length(pxx)/8));
pulse(i)=a;istore(i)=0;
a=10^(a/10);
fcap=i/(t*N)*2*pi;
sincap=sin(fcap.*n-(t-1)/2*fcap);
coscap=cos(fcap.*n-(t-1)/2*fcap);
pusi=-atan((sincap'*removal2)/(coscap'*removal2))/pi*180;
acap=2*sqrt(a/t);
estimate=-acap.*cos(fcap.*n+pusi/180*pi-(t-1)/2*fcap);
removal3_add=removal2+estimate;
removal3_sub=removal2-estimate;
if(removal3_add'*removal3_add>removal3_sub'*removal3_sub)
removal3=removal3_sub;
else
removal3=removal3_add;
end
figure(2)
subplot(224)
pxx=10*log10((abs(fft(removal3,t*N))).^2./t);
plot(freq,pxx.*istore+pulse);axis([0,0.5,-20,20]);hold on;grid on; title('9-9a PSD after third removal')
end
可以看到用这个方法可以在很大程度上减少泄露。

但是观察N=64的图，可以发现在0.2/T
附近的谱线的位置还是很多变，虽然相对之前的方法，我们可以很明显地看到一些分离开来的谱线，但是还是不理想。

N=256的图则效果好很多，0.2/T附近的谱线基本都是一致的，没有很大的变化，从这里也看出了样本的数目对结果还是有很大影响的。

% 3.minimum-leakage spectrum estimates
M=16;
f=linspace(0,0.5,1024);
s=exp(1j*2*pi*(0:(M-1))'*f);
for time1=1:3
index=N1:-1:M;
for time2=1:length(index)
X(time2,:)=data64(index(time2):-1:(index(time2)-M+1),time1).'; end
[temp1,temp2]=size(X);
R=1/temp1*(X'*X);
Sf(:,time1)=abs(M*1./diag((s.'*(R\conj(s)))));
end
figure;
for time1=1:3
plot(linspace(0,0.5,1024),10*log10(Sf(:,time1)));
hold on;
end
ylim([-30 10]);
xlim([0 0.5]);
title('9.10a minimum-leakage spectrum estimates with N=64 M=16');
hold off
比较图9-10a和9-10b，可以看出相同阶数时增大数据量能提高可靠性。

图9-10b~d可以看出滤波器阶数的增加会增大分辨力，但会使可靠性降低，高斯过程的起伏增大，0.1处的谱线峰值降低，这样我们选取滤波器阶数时需要一个折中，使频谱具有相对较好的分辨力和可靠性。

对比谱分析的结果，可以看到minimum-leakage spectrum estimates的优势还是比较明
显的，要好于参量方法中的AR建模方法。

%4.AR Method
function Sf=yw_AR_zty(data,M,f)
%write by for Y-W AR Method
%example£ºdata=data64(:,1);M=16;f=linspace(0,1,1000);plot(f,Sf);
f=f(:);
datalength=length(data);
%计算自相关
data=data(:);
temp=xcorr(data,'biased');
Rx=temp(datalength:(2*datalength-1),1)';
%初始化
b2=zeros(1,M);
a=zeros(M,M);
a(1,1)=-Rx(2)/Rx(1); b2(1)=(1-a(1,1)^2)*Rx(1);
%递推
for N=2:M
a(N,N)=-1*(Rx(N+1)+a(N-1,:)*[fliplr(Rx(2:N))
zeros(1,M-N+1)]')/b2(N-1);
b2(N)=(1-(a(N,N))^2)*b2(N-1);
for p=1:N-1
a(N,p)=a(N-1,p)+a(N,N)*a(N-1,N-p);
end
end
b2(M)=Rx(1)+Rx(2:(M+1))*a(M,:)';
%谱估计
Sf=b2(M)./abs(1+a(M,:)*exp(-1j*2*pi*(1:M).'*f.')).^2;
---------------------------------------------------------------------%9-11a
figure(2);
for time1=1:3
M=16;
f=linspace(0,1,1000);
Sf=yw_AR_zty(data64(:,time1),M,f);
plot(f,10*log10(Sf));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-11a Y-W AR spectrum estimates,N=64,M=16')
%9-11b
figure(3);
for time2=1:3
M=16;
f=linspace(0,1,1000);
Sf=yw_AR_zty(data256(:,time2),M,f);
plot(f,10*log10(Sf));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-11b Y-W AR spectrum estimates,N=256,M=16') %9-11c
figure(4);
for time2=1:3
M=32;
f=linspace(0,1,1000);
Sf=yw_AR_zty(data256(:,time2),M,f);
plot(f,10*log10(Sf));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-11b Y-W AR spectrum estimates,N=256,M=16')
用yw方法进行谱估计，结果显然不如最小泄漏方法。

对于N=64的数据，谱图（9-11a）可以看到在0.1处有尖峰，符合理想结果，但是对于0.2处的频谱分辨力低下。

当N=256,M=16时，频谱的分辨力还是没有得到大的提升，见图9-11b。

当N=256,M=32时（9-11c），频谱分辨力有了提升，但是也可以看到可靠性下降了，谱估计的能力不比之前的周期图和ML方法好。

比较低阶数时的状况（图9-15a、图9-15b），可以看出，频谱的分辨力很低，当M=7时，不能分辨出f=0.1处的频谱尖峰，从这里可以看到对于AR谱估计方法，模型阶数的确定对于谱分析结果很重要。

%5.Burg AR Method
%Burg function
function [Sf_burg a_now,sigma2]=burg_zty(data,model_order,f)
%written by to compute AR spectrum using burg method;
%model_order means model order;f:denote searching frequency range;
% for example,f=linspace(0,0.5,1024),model_order=16 or 32;
%Sf_burg=burg_zty(data64(:,time1),16,f);
%plot(f,10*log10(Sf_burg));
data=data(:);
f=f(:);
N=length(data);
a_pre=zeros(model_order,1);
a_now=a_pre;
ef_pre=data;%forward prediction error;
eb_pre=data;%backward prediction error;
ef_now=zeros(size(ef_pre));
eb_now=zeros(size(eb_pre));
for M=1:model_order
a_now(M)=-2*(ef_pre(M+1:N).'*eb_pre(M:N-1))/...
(ef_pre(M+1:N).'*ef_pre(M+1:N)+eb_pre(M:N-1).'*eb_pre(M:N-1));%(82)
if M>=2
a_now(1:M-1)=a_pre(1:M-1)+a_now(M)*a_pre(M-1:-1:1);%(79)
end
ef_now(M+1:N)=ef_pre(M+1:N)+a_now(M)*eb_pre(M:N-1);%(80a)
eb_now(M+1:N)=eb_pre(M:N-1)+conj(a_now(M))*ef_pre(M+1:N);%(80b)
ef_pre=ef_now;
eb_pre=eb_now;
a_pre=a_now;
end
sigma2=1/(N-model_order)*(ef_now(model_order+1:N).'*ef_now(model_orde r+1:N)+eb_now(model_order+1:N).'*eb_now(model_order+1:N))*0.5;%(78)
Sf_burg=sigma2./abs(1+a_now'*exp(-1j*2*pi*(1:model_order).'*f.')).^2; %(61)
---------------------------------------------------------------------%9-12a
figure(2);
for time1=1:3
model_order=16;
f=linspace(0,0.5,1024);
sf_burg=burg_zty(data64(:,time1),model_order,f);
plot(f,10*log10(sf_burg));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
% ylim([-20 20]);
ylim([-20 max(10*log10(sf_burg))+5]);
title('9-17a Burg AR spectrum estimates,N=64,M=16');
Burg 方法在估计谱线的相对强弱高度上，好于周期图法。

在N=64下（图9-12a），Burg方法对谱线有更好地分辨能力。

在N=256下，周期图法对谱线有更好地分辨能力。

当M=44时(图9-16a)，此时Burg方法的精确性就不如周期图法了。

比较来看对AR方法，阶数M=32是比较合适的估计阶数，有着相对好的分辨力和可靠性。

%6.FB AR Method
%FB AR Method function
%Example: plot figure 9-13 and 9-14
%9-13-a
%fb_ls(data,64,16,[64 128 192],3,'figure 9-13-a');
%9-13-b
%fb_ls(data,256,16,[256 512 768],3,'figure 9-13-b');
%9-13-c
%fb_ls(data,256,32,[256 512 768],3,'figure 9-13-b');
%9-14 <--This one is different with 9-14 in book
%fb_ls(data,64,24,64,1,'figure 9-14');
function fb_ls(data,N,M,START,times,title_name)
K=N-M;
xb=zeros(M+1,N-M);
rx=zeros(M+1,M+1,times);
for time = 1:times
xf=toeplitz(data(START(time)+M+1:-1:START(time)+1),data(START(time)+M +1:START(time)+N));
for time1=1:M+1
xb(time1,:)=data(START(time)+time1:START(time)+N-M+time1-1);
end
rx(:,:,time)=1/K*(xf*xf'+xb*xb');
end
f=0:0.001:0.5;
exponent=zeros(M,length(f));
for p=1:M
exponent(p,:)=exp(-1i*2*pi*f*p);
end
figure
a=zeros(M,1,times);
b=zeros(times);
P=zeros(1,length(f),times);
for time=1:times
a(:,:,time)=rx(2:M+1,2:M+1,time)\(-rx(2:M+1,1,time));
b(time)=rx(1,:,time)*[1 (a(:,:,time))']';
P(:,:,time)=b(time)./((abs((1+(a(:,:,time)'*exponent)))).^2);
plot(f,10*log10(abs(P(:,:,time))));
hold on;
end
hold off
xlabel('frequency');
ylabel('spectral density estimate(dB)'); xlim([0 (1/2)]);
ylim([-20 30]);
title(title_name);
end
FB LS模型能较好地分辨出0.1和0.2处的三条谱线，且谱线的位置也比较准确，但是三个谱线的相对高度有误差，高斯噪声的估计存在较多的假峰，可见这种方法的可靠性并不高。

在N=256时性能不如改进的周期图法。

对于FB LS模型，它也是一种AR建模方法，因此和之前的几种AR方法一样，模型阶数估计对谱估计结果至关重要
%7.ODNE AR Method
%ODNE AR function(biased)
function Sf=odne_biased_zty(data,M,Q,f)
%written by to compute AR spectrum by solving Q-set overdetermined normal equation;
%M:modelorder;f:denote searching frequency range;
%for example: Q=48,M=16 or 32, f=linspace(0,1,1000);
%Sf=odne_biased_zty(data256(:,time1),32,Q,f);
%plot(f,10*log10(Sf));
data=data(:);
f=f(:);
N=length(data);
r=xcorr(data,'biased');%(67)
r=r(N:N+Q);%~Rx(0:Q);
Rx=toeplitz(r(1:Q),r(1:M));
rx=r(2:end);
a=-(Rx'*Rx)\Rx'*rx;%(90)
b2=r(1)+r(2:M+1)'*a;%(69)
Sf=b2./(abs(1+a'*exp(-1j*2*pi*(1:M).'*f.'))).^2;
--------------------------------------------------------------------- %ODNE AR function(unbiased)
function Sf=odne_unbiased_zty(data,M,Q,f)
data=data(:);
f=f(:);
N=length(data);
r=xcorr(data,'unbiased');%(67)
r=r(N:N+Q);%~Rx(0:Q);
Rx=toeplitz(r(1:Q),r(1:M));
rx=r(2:end);
a=-(Rx'*Rx)\Rx'*rx;%(90)
b2=r(1)+r(2:M+1)'*a;%(69)
Sf=b2./(abs(1+a'*exp(-1j*2*pi*(1:M).'*f.'))).^2;
---------------------------------------------------------------------%9-17a
figure(2);
for time1=1:3
M=16;
Q=48;
f=linspace(0,1,1000);
Sf=odne_biased_zty(data64(:,time1),M,Q,f);
plot(f,10*log10(Sf));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-17a ODNE AR spectrum estimates,N=64,Q=48,M=16');
%9-18b
figure(3);
for time1=1:3
M=16;
Q=48;
f=linspace(0,1,1000);
Sf=odne_unbiased_zty(data256(:,time1),M,Q,f);
plot(f,10*log10(Sf));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-30 max(10*log10(Sf))+10]);
title('9-17a ODNE AR spectrum estimates,N=64,Q=48,M=16');
随着N所取的值变大（如图9-17a、9-17c），从N=64变化到N=256，可以明显看见正弦函数的估计谱线从平滑不可见到逐渐变得清晰可辨，分辨力提升。

由于是用估计的协方差矩阵来代替，当M越大可取的数据越少，方差越高，但是偏置会减小；而N越大，所取的数据越多，方差越小。

因此为了获得更小的方差，常常要求M<<N。

由于从有偏估计到无偏估计，其均值的估计产生了偏差，但是可以换来方差的减小。

图9-18a~c为ODNE无偏估计，可以看到图中两个尖峰的中心位置是发生了挪移的，但是其波动要比无偏估计的要好，这是由于从有偏估计到无偏估计，其均值的估计产生了偏差，但是可以换来方差的减小。

同时可以看到，有偏估计由于主瓣中心离开了正确位置导致了其宽度的变大。

可以看到无偏自相关矩阵估计替代有偏估计，在尖峰估计中，其分辨力更好。

%8.SVD AR Method
function Sf=svd_AR_zty(data,Q,M0,f)
%write by for computing SVD AR Method
%example M0=16;Q=48;data=data64(:,1);f=0:0.001:1;plot(f,Sf);
%(1)自相关矩阵计算
N=length(data);
data=data(:);
for pq=1:N
Rx(1,pq)=(data(1:(N-pq+1),1)'*data(pq:N,1))/N;
end
matrixRx=toeplitz(Rx(1,1:Q));%相关矩阵
[W,A,V]=svd(matrixRx);%SVD分解
%(2)参量估计
rx=-(Rx(1,2:Q+1))';
lamda=diag(A);
lamda=[1./lamda(1:M0);zeros(Q-M0,1)];
a=W*diag(lamda)*V'*rx;%(100);
b2=Rx(1)+Rx(2:Q+1)*a;%(69);
%(3)谱估计
f=f(:);
Sf=b2./abs(1+a'*exp(-1j*2*pi*(1:Q).'*f.')).^2;
---------------------------------------------------------------------%9-19a
figure(2)
for time1=1:3
Q=48;
M0=16;
f=linspace(0,1,1000);
Sf=svd_AR_zty(data64(:,time1),Q,M0,f);
plot(f,10*log10(Sf));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-19a SVD AR spectrum estimates,N=64,M0=16')
%9-19b
figure(3)
for time1=1:3
Q=48;
M0=32;
f=0:0.001:pi;
Sf=svd_AR_zty(data64(:,time1),Q,M0,f);
plot(f,10*log10(Sf));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-19b SVD AR spectrum estimates,N=64,M0=32')
%9-19c
figure(4)
for time2=1:3
Q=48;
M0=16;
f=0:0.001:pi;
Sf=svd_AR_zty(data256(:,time2),Q,M0,f);
plot(f,10*log10(Sf));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-19a SVD AR spectrum estimates,N=256,M0=16') %9-19d
figure(5)
for time2=1:3
Q=48;
M0=32;
f=0:0.001:pi;
Sf=svd_AR_zty(data256(:,time2),Q,M0,f);
plot(f,10*log10(Sf));
hold on;
end
hold off;
xlabel('frequency');
ylabel('spectral density estimate(dB)');
xlim([0 (1/Ts/2)]);
ylim([-20 20]);
title('9-19a SVD AR spectrum estimates,N=256,M0=32')
从图9-19a~d来看，谱分析结果不如之前其他的参量估计和非参量估计：频谱的分辨力不高，高斯过程段假峰很多。

SVD方法主要适用范围：相关矩阵的特征值有陡降的变化，有利于得到正确的阶数估计，它对频谱尖峰的检测尤为适用。

而当前分析的实验数据的相关矩阵的特征值并没有这种陡降的特征，因此模型阶数估计并不准确，相关矩阵的对微扰的敏感性还是很大，如此得到的谱估计结果的稳定性很差，分辨力也不高。

%9.Hybrid Method
clear all;
close all;
clc;
load experiment_data;
N1=64;%data length;
N2=256;
data64=[data(65:65+N1-1) data(129:129+N1-1) data(193:193+N1-1)];%the three statistical samples of data using starting points of 64,128,and 192 for N=64;
data256=[data(257:257+N2-1) data(513:513+N2-1)
data(769:769+N2-1)];%the three statistical samples of data using starting points of 256,512,and 768 for N=256;
figure;
subplot(3,1,1);
plot(data);
title('orignal data,1024 samples');
subplot(3,1,2);
for time1=1:3
plot(data64(:,time1)+time1*6);
hold on;
end
title('truncated data, 64 samples');
xlim([1 1024]);
subplot(3,1,3);
for time1=1:3
plot(data256(:,time1)+time1*6);
hold on;
end
title('truncated data, 256 samples');
xlim([1 1024]);
Ts=1;%sampling interval is 1s;
%%
%hybrid method with data length of 64;fig 9.20;
%for example: line spectrum number:3; then enter:217, 206, 103; repeat 3
%times.
close all;
K=16;%zero padding factor;
N=64;%data point;
f=linspace(0,1,K*N+1);
model_order=16;
h1=figure;%used to plot temp figure;
h2=figure;%used to plot final spectrum;
for time1=1:3
%step1.FB method is used to estimate frequency;
sf=fb_yzb(data64(:,time1),model_order,f);
figure(h1);
plot(1:(length(sf)+1)/2,10*log10(sf(1:(length(sf)+1)/2)));
freq_num = str2num(input('input the number of line spectrum:',
's'));%number of line spectrum;
f_index=zeros(freq_num,1);%value to save line spectrum frequency index; for time2=1:freq_num
f_index(time2)=str2num(input(['input the ',num2str(time2),' frequency index of line spectrum:'], 's'));%number of line spectrum;
end
%step2:subtract line spectrum from data;
data_s=data64(:,time1);%data after subtraction;
for time2=1:freq_num
sf=(abs(fft(data64(:,time1),K*N1)).^2/length(data64(:,time1)));
% Sf_64=2*(1/(N*Ts)*per_gram(data_s,N*K,Ts)).^0.5;%(235) Sf_64=2*(1/(N*Ts)*sf).^0.5;
a_est(time1,time2)=(Sf_64(f_index(time2)));%estimation of a for f;
y=-sum(data_s.*sin(2*pi*f(f_index(time2))*[1:length(data_s)].'*Ts));% (233)
x=sum(data_s.*cos(2*pi*f(f_index(time2))*[1:length(data_s)].'*Ts));
phi_est(time1,time2)=atan(y/x);%(233)
if x<0
phi_est(time1,time2)=phi_est(time1,time2)+pi;
end
data_s=data_s-a_est(time1,time2)*cos(2*pi*f(f_index(time2))*[1:length (data_s)].'*Ts+phi_est(time1,time2));
end
%step3:M=4,K=8 smooth periodagram;
M=4;
%frequency smoothing;
temp2=abs(fft(data_s,N1*K)).^2/N1;
c=[temp2((K*M/2):(length(temp2)/2+1));zeros(K*M/2-1,1)];
r=[flipud(temp2(1:(K*M/2)));zeros(K*M/2,1)];
temp=toeplitz(c,r);
Sf_recover=temp*ones(K*M,1)*1/(K*M);
%step 4: plot the total spectrum;
for time2=1:freq_num
Sf_recover(f_index(time2))=Sf_recover(f_index(time2))+(a_est(time1,ti me2)/2)^2*N;
end
figure(h2);
plot(((1:((N*K/2+1)))-1)*1/((N*K)*Ts),10*log10(Sf_recover(1:N*K/2+1)) );
hold on;
end
figure(h2);
xlabel('frequency');
ylabel('spectral density estimate(dB)');
ylim([-20 20]);
title('fig 9.20:hybrid spectrum estimate');
saveas(gcf,'9.20 hybrid spectrum estimate.emf');
方法，使得到的谱分析结果有较好的分辨力，又有较好的可靠性。