✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab仿真内容点击👇
⛄ 内容介绍
信息科学的一项重要任务是获取客观世界的真实信息,然而对于任何一个系统,必然存在噪声,而当所要检测的信号比较微弱且淹没在强噪声背景中时,用传统的检测方法检测信号是很困难的,因此如何把淹没于强噪声中的有用信号提取出来的问题越来越引起人们的关注.混沌理论表明一类混沌系统在一定条件下对小信号具有参数敏感性同时对噪声具有免疫,因此使得它在信号检测中非常具有潜力. 本文详细研究了利用Duffing混沌振子来进行微弱信号检测的原理和过程,运用大量实验数据说明了参数对Duffing振子用于微弱信号检测的影响,并选择合适的参数对微弱信号进行了检测.实验结果表明,本文的方法对于微弱信号,特别是淹没在强噪声背景下的信号,检测的效果明显优于传统的时域检测方法.本文的最后依据Duffing振子微弱信号检测原理,利用虚拟仪器实现了该检测系统.
⛄ 完整代码
% Duffing混沌模型状态方程如下：
% x'' + k * x' - x(t) + x(t)^3 = a * cos(ωt)
% - x(t) + x(t)^3是非线性恢复力
% a * cos(ωt)是周期策动力
% k代表阻尼比
% 降阶：y = x'
% 把原二阶微分方程化为一阶微分方程组：
% y' = -k * y + x - x^3 + a * cos(ωt)
% f_1(x,y) = y
% f_2(x,y) = -k * y + x - x^3 + a * cos(ωt)
% 在混沌信号基础上加入具有信噪比的待测信号
% 其幅度为a_sig，角频率为omega_sig，信噪比为SNR
% k = 0.5 omega = 1 时
% 无噪声的情况下 周期策动系数临界值a = 0.80207
% 噪声方差为0.08的情况下 周期策动系数临界值a = 0.8002 输入信号幅度敏感值a_sig = 0.0145
% 这时输入信号信噪比为SNR = -25.8035 dB 检测门限为 -18.3863 dB
% 下面是比较好的数据
% 噪声方差为0.01的情况下 周期策动系数临界值a = 0.798 输入信号幅度敏感值a_sig = 0.003
% 这时输入信号信噪比为SNR = -33.4679 dB 检测门限为 -25.2288 dB
% 噪声方差为0.08的情况下 周期策动系数临界值a = 0.8002 输入信号幅度敏感值a_sig = 0.008
% 这时输入信号信噪比为SNR = -30.9691 dB 检测门限为 -20.9691 dB
clear;
clc;
% 采样相关参数
h = 0.01; % 采样步长
t_0 = 0;
t_n = 1000;
t = t_0 : h : t_n;
n = length(t);
y_duffing = zeros( n, 2 );
U_weaksig = zeros( n, 1 );
% 混沌系统相关参数
k = 0.5; % 这里再手动设置一下阻尼比，用于计算R-K参数
a = 0.8002; % 周期策动系数
omega = 1; % 周期策动角频率
% 待测信号与噪声信号相关参数
a_sig = 0.008; % 待测信号幅度
omega_sig = 1; % 待测信号角频率
var_noi = 0.08; % 噪声方差
% 为了在信号中加入噪声，在迭代的模块前提前生成待测信号
for j = 1 : n
t_1 = t_0 + h;
U_weaksig(j) = a_sig * cos( omega_sig * t_0 ); % 生成待测的微弱信号
t_0 = t_1;
end
% 在待测信号里加入噪声
U_noise = sqrt( var_noi ) * randn( n, 1 );
U_weaksig_SNR = U_weaksig + U_noise;
% 计算信噪比和检测门限
SNR = 10 * log10( ( 0.5 * a_sig^2 ) / ( 0.5 * var_noi ) );
ThresholdDetection = 10 * log10( a_sig );
% 设置混沌系统初始值
x_0 = 0;
y_0 = 10^-10;
y_duffing( 1, : ) = [ x_0, y_0 ];
% 生成待测信号完成之后，重新初始化t的值，方便下面的迭代
t_0 = 0;
t_1 = 0;
% 开始迭代
for i = 1 : n
t_1 = t_0 + h;
% 当前x值：y( i , 1 )；
% 当前y值：y( i , 2 )；
U = a * cos( omega * t_0 ) + U_weaksig_SNR(i);
% U = a * cos( omega * t_0 ) + U_weaksig(i); % 不加噪声的测试句
% 计算K值
K11 = y_duffing( i , 2 );
K12 = ( 1 + h/2 ) * y_duffing( i , 2 );
K13 = ( 1 + h/2 + h^2/4 ) * y_duffing( i , 2 );
K14 = ( 1 + h + h^2 / 2 + h^3/4 ) * y_duffing( i , 2 );
K21 = -1 * k * y_duffing( i , 2 ) + y_duffing( i , 1 ) - y_duffing( i , 1 )^3 + U;
K22 = ( -1 * k + 1/2 * h * k^2 ) * y_duffing( i , 2 ) ...
+ ( 1 - 3/4 * h^2 - 1/2 * h * k ) * y_duffing( i , 1 ) ...
- 3/2 * h * y_duffing( i , 1 )^2 ...
+ ( -1 + 1/2 * h * k ) * y_duffing( i , 1 )^3 ...
+ ( 1 - 1/2 * h * k ) * U ...
+ ( 1/2 * h - 1/8 * h^3 );
K23 = ( -1 * k + 1/2 * h * k^2 - 1/4 * h^2 * k^3 ) * y_duffing( i , 2 ) ...
+ ( 1 - 3/4 * h^2 - 1/2 * h * k + 1/4 * h^2 * k^2 + 3/8 * h^3 * k ) * y_duffing( i , 1 ) ...
+ ( - 3/2 * h + 3/4 * h^2 * k ) * y_duffing( i , 1 )^2 ...
+ ( -1 + 1/2 * h * k - 1/4 * h^2 * k^2 ) * y_duffing( i , 1 )^3 ...
+ ( 1 - 1/2 * h * k + 1/4 * h^2 * k^2 ) * U ...
+ ( 1/2 * h - 1/8 * h^3 - 1/4 * h * k + 1/16 * h^4 * k );
K24 = ( -1 * k + h * k^2 - 1/2 * h^2 * k^3 + 1/4 * h^3 * k^4 ) * y_duffing( i , 2 )...
+ ( 1 - 3 * h^2 - h * k + 1/2 * h^2 * k^2 + 3/4 * h^3 * k - 1/4 * h^3 * k^3 - 3/8 * h^4 * h^2 ) * y_duffing( i , 1 )...
+ ( -3 * h + 3/2 * h^2 * k - 3/4 * h^3 * k^4 ) * y_duffing( i , 1 )^2 ...
+ ( -1 + h * k - 1/2 * h^2 * k^2 + 1/4 * h^3 * k^3 ) * y_duffing( i , 1 )^3 ...
+ ( 1 - h * k + 1/2 * h^2 * k^2 - 1/4 * h^3 * k^3 ) * U ...
+ ( h - h^3 - 1/2 * h^2 * k + 1/4 * h^2 * k^2 + 1/8 * h^4 * k - 1/16 * h^5 * k^2 );
% 对函数值进行迭代
y_duffing( i + 1 , 1 ) = y_duffing( i , 1 ) + h/6 * ( K11 + 2 * K12 + 2 * K13 + K14 );
y_duffing( i + 1 , 2 ) = y_duffing( i , 2 ) + h/6 * ( K21 + 2 * K22 + 2 * K23 + K24 );
t_0 = t_1;
end
% 计算Melnikov函数
Ms = zeros( 1 , n );
dMs = ( y_duffing( ( 2 : ( n + 1 ) ) , 2 ) )' .* ( a * cos( omega * t ) - 0.5 * omega * ( y_duffing( ( 2 : ( n + 1 ) ) , 2 ) )' );
NN = ( round( 2 * pi * 100 ) / 100 ) / h;
% NN = 6.28 / h;
for p = 0 : ( n - 1 )
if( ( p - NN ) > 0 )
for q = ( p - NN + 1 ) : p
Ms(p) = Ms(p) + dMs(q) * h;
end
else
for q = 1 : p
Ms(p) = Ms(p) + dMs(q) * h;
end
end
end
% 基于间接法对输入信号进行功率谱估计
y_corr_in = xcorr( U_weaksig_SNR, 'unbiased' ); % 输入信号功率谱分析的测试句
y_fft_corr_in = fft( y_corr_in, n );
P_in = abs( y_fft_corr_in );
index_in = 0 : round( n / 2 - 1 );
k_P_in = index_in * 1 / ( n * h );
P_dB_in = 10 * log10( P_in( index_in + 1 ) );
% 基于间接法对输出信号进行功率谱估计
y_corr = xcorr( y_duffing( 2 : ( n + 1 ) , 1 ), 'unbiased' );
y_fft_corr = fft( y_corr, n );
P = abs( y_fft_corr );
index = 0 : round( n / 2 - 1 );
k_P = index * 1 / ( n * h );
P_dB = 10 * log10( P( index + 1 ) );
% 画图
subplot( 2, 2, 1 );
plot ( y_duffing( : , 1 ), y_duffing( : , 2 ) );
title('基于4阶Runge Kutta离散化的Duffing系统相图');
subplot( 2, 2, 2 );
plot ( y_duffing( : , 1 ) );
title('Duffing系统时域图');
subplot( 2, 2, 3 );
plot ( U_weaksig_SNR );
hold on;
plot ( U_weaksig );
legend('带噪声的待测信号','不带噪声的待测信号');
title('待测信号（带噪声）时域图');
subplot( 2, 2, 4 );
plot ( t, Ms );
ylim ( [ -1 1 ] )
title('Melnikov函数图');
figure;
subplot( 2, 1, 1 );
plot( k_P_in, P_dB_in );
title('输入信号功率谱分布图');
subplot( 2, 1, 2 );
plot( k_P, P_dB );
title('输出信号功率谱分布图');
⛄ 运行结果
⛄ 参考文献
[1] 赵文礼, 黄振强, 赵景晓. 基于Duffing振子的微弱信号检测方法及其电路的实现[J]. 电路与系统学报, 2011, 16(6):6.
[2] 路鹏. 基于Duffing振子的微弱信号检测[D]. 吉林大学.
[3] 孙玉胜, 詹小霞, 石军,等. 基于Duffing振子的微弱信号检测[J]. 郑州轻工业学院学报：自然科学版, 2012, 27(5):4.
[4] 韩振榕. 基于Duffing振子的微弱信号检测研究[D]. 苏州大学.
