clc
clear
close all
N = 128*64; %定义总采用点数
n = [0:N-1]; %定义总采用点数n,是一个[1×8192]的向量
Nfft = 128; %定义fft的点数
Fs = 1; %定义信号采用速率
t = n/Fs; %定义信号持续时间t,是一个[1×8192]的向量
f1 = 0.6381/2/pi; %定义有用信号的第1个频率f1
f2 = 0.8345/2/pi; %定义有用信号的第2个频率f2
f3 = f1+f2; %定义有用信号的第3个频率f3
f4 = 0.4909/2/pi; %定义有用信号的第4个频率f4
f5 = 1.7671/2/pi; %定义有用信号的第5个频率f5
f6 = f4+f5; %定义有用信号的第6个频率f6
fai1 = pi/6; %定义有用信号的第1个初相fai1
fai2 = pi/3; %定义有用信号的第2个初相fai3
fai3 = fai1+fai2; %定义有用信号的第3个初相fai3
fai4 = 5*pi/4; %定义有用信号的第4个初相fai4
fai5 = 6*pi/4; %定义有用信号的第5个初相fai5
fai6 = pi; %定义有用信号的第6个初相fai6
s1 = cos(2*pi*f1*t + fai1); %定义有用信号的第1个分量信号
s2 = cos(2*pi*f2*t + fai2); %定义有用信号的第2个分量信号
s3 = cos(2*pi*f3*t + fai3); %定义有用信号的第3个分量信号
s4 = cos(2*pi*f4*t + fai4); %定义有用信号的第4个分量信号
s5 = cos(2*pi*f5*t + fai5); %定义有用信号的第5个分量信号
s6 = cos(2*pi*f6*t + fai6); %定义有用信号的第6个分量信号
s = s1+s2+s3+s4+s5+s6; %定义有用信号的表达式,即输入的有用信号
% s = awgn(s,0,'measured');
noi = zeros(1,N); %定义信号noi,是一个[1×8192]的0向量
%高斯白噪声 %rand函数产生由在(0, 1)之间正态分布的随机数组成的数组
noi = randn(1,N); %定义噪声信号noi,是一个[1×8192]的随机向量
s = s + noi; %定义输入的混合信号表达式,即有用信号与高斯噪声信号的叠加
%有色噪声
nys = 0;
if (nys == 1)
Ln = length(s);
Nd = [1 -0.5 0.7 0.1];
Nc = [1 0.5 0.2]; %[Bspec,waxis] = bispecd (y,nfft, wind, nsamp, overlap)
nd = length(Nd)-1; %y=s,表示输入信号
nc = length(Nc)-1; %nfft=128,表示fft的点数
xik = zeros(nc,1); %wind=0, 用Parzen窗,wind=1,用hexagonal窗(六角),wind=5,用最优窗。
ek = zeros(nd,1); %nsamp=64,每段样本的数量
xi = randn(Ln,1); %overlap=32,每段样本的重叠数量
for kn=1:Ln
e(kn)=-Nd(2:nd+1)*ek+Nc*[xi(kn);xik];
for nn=nd:-1:2
ek(nn)=ek(nn-1);
end
ek(1)=e(kn);
for nn=nc:-1:2
xik(nn)=xik(nn-1);
end
xik(1)=xi(kn);
end %产生有色噪声
s = s+0.25*e; %有用信号与高斯噪声叠加后的信号s再加上非高斯有色噪声
noi = 0.25*e; %定义噪声信号noi为混合噪声信号
end %有色噪声
%信噪比
S_N= 10*log10(sum(abs(s).^2)/sum(abs(noi).^2)); %定义输入信噪比
[bspec,waxis] = bispeci (s,16,64,32,0,Nfft,0); %调用间接法求双阶谱估计子函数
close all
%双谱信号(包含相干信号的幅度和相位信息,幅度受信号的功率谱影响)
SS = Nfft/2+1;
waxis1 = waxis(SS:end)*Fs;
figure
subplot(121) %11双谱估计二维等高线图
bspec1 = bspec(SS:end,SS:end);
contour(waxis1,waxis1,abs(bspec1),4);
grid on
title('双谱估计二维等高线图')
xlabel('f1/Hz'), ylabel('f2/Hz')
[X1,Y1] = meshgrid(waxis1,waxis1);
subplot(122) %12双谱估计三维图
mesh(X1,Y1,abs(bspec1))
axis tight
xlabel('f1/Hz')
ylabel('f2/Hz')
zlabel('幅度')
title('双谱估计三维图');
figure %2原始时域信号图
plot(t,s);
axis([0,200,-8,8])
xlabel('t')
ylabel('幅度')
title('时域信号')
P = 10*log10(abs(fft(s-mean(s),Nfft).^2)/(Nfft+1));
f = (0:length(P)-1)*Fs/(length(P));
figure %3功率谱估计图
plot(f(1:Nfft/2),P(1:Nfft/2));
title('功率谱估计');
xlabel('f/Hz')
ylabel('P/dB')
%相干双谱信号(频率f1,f2非线性相位耦合产生的能量在f1+f2处总能量中所占的比例)
mask = hankel([1:Nfft],[Nfft,1:Nfft-1] ); % the hankel mask (faster)
P_d = abs(fft(s,Nfft).^2)/(Nfft+1);
P_d = P_d';
bspec_d = zeros(Nfft,Nfft);
Bspec_d = zeros(Nfft,Nfft);
bspec_d = (P_d * P_d') .*reshape(P_d(mask), Nfft, Nfft);
Bspec_d = abs(bspec)./sqrt(bspec_d);
figure %4双相干谱估计二维等高线图
subplot(121)
contour(waxis1,waxis1,abs(Bspec_d(SS:end,SS:end)),4);
grid on
title('双相干谱估计二维等高线图')
xlabel('f1/Hz'), ylabel('f2/Hz')
subplot(122) %4双相干谱估计三维图
mesh(X1,Y1,abs(Bspec_d(SS:end,SS:end)));
axis tight
xlabel('f1/Hz')
ylabel('f2/Hz')
zlabel('幅度')
title('双相干谱估计三维图');
figure
subplot(211) %51双相干谱估计三维图
diagBspec = diag(fliplr(bspec));
plot(waxis1,abs(diagBspec(SS:end)));
grid on;
xlabel('f/Hz')
ylabel('幅度')
title('双谱估计切片图')
subplot(212) %51双相干谱估计三维图
diagBspec_d = diag(fliplr(Bspec_d));
plot(waxis1,abs(diagBspec_d(SS:end)));
grid on;
xlabel('f/Hz')
ylabel('幅度')
title('双相干谱估计切片图')
评论3