|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
在沈善德书《电力系统辨识》中p39有伪随机序列的功率密度谱曲线,我是竭尽我所能了,就是得不到他的那种形状啊!在此请教此方面的高手指点。
伪随机序列我是由程序genepncode生成的,程序如下:
function [Y,m]=genpncode(sm,a,ts,num)
% 产生长度等于 pow2(sm)-1 m序列,序列的幅值为a,Y为输出的M序列,m为M序列的周期等于 pow2(sm)-1
l=pow2(sm)-1;
num=l*num;
n=[1 0 1 0 1 0 1 0 1 0 1 0];
if sm<3
errordlg('错误:sm不应小于3');
return;
end
if sm>10
errordlg('错误:sm不应大于10');
return;
end
for i=1:num
temp=0;
% 选择多项式
if sm==3
temp=n(1)+n(3);
end
if sm==4
temp=n(1)+n(4);
end
if sm==5
temp=n(2)+n(5);
end
if sm==6
temp=n(1)+n(6);
end
if sm==7
temp=n(3)+n(7);
end
if sm==8
temp=n(2)+n(3)+n(4)+n(8);
end
if sm==9
temp=n(9)+n(4);
end
if sm==10
temp=n(10)+n(3);
end
% mod(2)
temp=mod(temp,2);
for j=sm:-1:2
n(j)=n(j-1);
end
n(1)=temp;
y(i)=n(sm);
if y(i)==0
y(i)=-a;
else
y(i)=a;
end
end
Y=y;
m=l;
请高人对我的分析谱的程序指点下:
clc;
clear all;
sm=7; %寄存器个数
a=1; %幅值
ts=1/12.7; %间隔时间
%n=2; %产生n个周波的伪随机信号
[u,m]=genpncode(sm,a,ts,1); %得到M序列为u,序列的周期为m,M序列的幅值为a,寄存器的个数为sm
%l=length(u)
%自己加的几个语句,作用是将每个采样点重复K次再频谱分析
x=u;
w1=u;
k=10;
n=length(u);
for i=1:n
for j=0:(k-1)
x(1+k*(i-1)+j)=w1(i);
end
end
%间接法,现求取自相关函数,再对自相关函数fft,求取功率谱===========================
fs=127;
N=length(u);
Nfft=256;
maxlags=100;
[cxn,lags]=xcorr(u,maxlags,'unbiased'); %计算序列的自相关函数
figure(1)
plot(lags/fs,cxn)
CXk=fft(cxn,Nfft);
pxx4=abs(CXk);
f=(0:(Nfft/2-1))*fs/Nfft;
plot_Pxx=10*log10(pxx4(1:(Nfft/2)));
figure(2)
plot(f,plot_Pxx);
谢谢了先!!! |
|