clc
clear
close?all
N?=?128*64;
n?=?（0:N-1）;
Nfft?=?128;
Fs?=?1;?%?采樣頻率假定為1
t?=?n/Fs;
f1?=?0.6381/2/pi;
fai1?=?pi/6;
s0?=?10*cos（2*pi*f1*t?+?fai1）;
%?%?得到SNR為-20時候對應的雙譜切片
SNR?=?-20;
noi?=?awgn（s0SNR（1）0）?-?s0;
s0?=?s0?+?noi;

%?使用直接雙譜估計方法（基于fft），bspec:雙譜估計，一個Nfft*Nfft的數組;
[bspecwaxis]?=?bispecd（s0?Nfft?1?64?0）;
close?all

%雙譜信號（包含相干信號的幅度和相位信息，幅度受信號的功率譜影響）
SS?=?Nfft/2+1;
waxis1?=?waxis（SS:end）*Fs;
bspec1?=?bspec（SS:endSS:end）;
[X1Y1]?=?meshgrid（waxis1waxis1）;
mesh（X1Y1abs（bspec1））
axis?tight
xlabel（‘f1/Hz‘）
ylabel（‘f2/Hz‘）
zlabel（‘幅度‘）
title（‘雙譜估計三維圖‘）;

figure（）
diagBspec?=?diag（fliplr（bspec））;?%?輸出主對角線上的元素
plot（waxis1abs（diagBspec（SS:end）））;
grid?on;
xlabel（‘f/Hz‘）
ylabel（‘幅度‘）
title（‘雙譜估計切片圖‘）
SNR=（-20:10）;
Pf=zeros（1length（SNR））;
Pd=zeros（1length（SNR））;
threshhold?=?0.004;?%?為了取得較大的檢測概率，選取較小的閥