MATLAB谐波信号双谱估计

clear

close all

N = 128*;

n = [0:N-1];

Nfft = 128;

Fs = 1;

t = n/Fs;

f1 = 0.6381/2/pi;

f2 = 0.8345/2/pi;

f3 = f1+f2;

f4 = 0.4909/2/pi;

f5 = 1.7671/2/pi;

f6 = f4+f5;

fai1 = pi/6;

fai2 = pi/3;

fai3 = fai1+fai2;

fai4 = 5*pi/4;

fai5 = 6*pi/4;

fai6 = pi;

s1 = cos(2*pi*f1*t + fai1);

s2 = cos(2*pi*f2*t + fai2);

s3 = cos(2*pi*f3*t + fai3);

s4 = cos(2*pi*f4*t + fai4);

s5 = cos(2*pi*f5*t + fai5);

s6 = cos(2*pi*f6*t + fai6);

s = s1+s2+s3+s4+s5+s6;

% s = awgn(s,0,'measured');

noi = zeros(1,N);

%高斯白噪声

noi = randn(1,N);

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];

nd = length(Nd)-1;

nc = length(Nc)-1;

xik = zeros(nc,1);

ek = zeros(nd,1);

xi = randn(Ln,1);

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;

noi = 0.25*e;

end

%信噪比

S_N= 10*log10(sum(abs(s).^2)/sum(abs(noi).^2));

[bspec,waxis] = bispecd(s,Nfft,1,,32);

close all

%双谱信号（包含相干信号的幅度和相位信息，幅度受信号的功率谱影响）

SS = Nfft/2+1;

waxis1 = waxis(SS:end)*Fs;

figure

subplot(121)

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)

mesh(X1,Y1,abs(bspec1))

axis tight

xlabel('f1/Hz')

ylabel('f2/Hz')

zlabel('幅度')

title('双谱估计三维图');

figure

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()

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()

subplot(121)

contour(waxis1,waxis1,abs(Bspec_d(SS:end,SS:end)),4);

grid on

title('双相干谱估计二维等高线图')

xlabel('f1/Hz'), ylabel('f2/Hz')

subplot(122)

mesh(X1,Y1,abs(Bspec_d(SS:end,SS:end)));

axis tight

xlabel('f1/Hz')

ylabel('f2/Hz')

zlabel('幅度')

title('双相干谱估计三维图');

figure()

subplot(211)

diagBspec = diag(fliplr(bspec));

plot(waxis1,abs(diagBspec(SS:end)));

grid on;

xlabel('f/Hz')

ylabel('幅度')

title('双谱估计切片图')

subplot(212)

diagBspec_d = diag(fliplr(Bspec_d));

plot(waxis1,abs(diagBspec_d(SS:end)));

grid

on;

xlabel('f/Hz')

ylabel('幅度')

title('双相干谱估计切片图')