matlab代码

%%%%随机共振参数优化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% close all;
clear;
clc;
% tic
N =200;                        % 每个码元的持续时间
fs=10;                          % 随机共振，龙格库塔算法频率
% n = 10;
%rr = zeros(n,100000);
%M = 1;
SNR_index = -13;                 % 输入信号的信噪比 dB
%load('ns_test');            % 读取码元文件
load('bits_2000')
%%%%load('sim_test')
% bit_num = 100;                  %码元个数
% bits = rand_seed(bit_num);
% save('bits');
%noise_amp =1            % 加性方式添加噪声
 
% AN= 4;          % 噪声范围参数
% noise_ampx=0:0.01:AN;
% AN_x =  noise_ampx.*5;      % 噪声的分辨率

%%  选择调制方式，由数字码元生成波形
signal_set = 'BNRZ'; % either 'BNRZ', 'bpsk2', 'ask', or 'fsk'

if strcmp(signal_set, 'BNRZ')
    signal1 =0.3*ones(1,N);
    signal0 = -signal1;
elseif strcmp(signal_set, 'bpsk1')
    signal1 = ones(1,N);
    signal0 = -signal1;
elseif strcmp(signal_set, 'bpsk2')
    signal1 = sqrt(2)*sin(2*pi*[0:N-1]/N);
    signal0 = -signal1;
elseif strcmp(signal_set,'fsk')
    signal1 = sqrt(2)*sin(2*pi*2*[0:N-1]/N);
    signal0 = sqrt(2)*sin(2*pi*3*[0:N-1]/N);
elseif strcmp(signal_set, 'ask')
    signal1 = 0.3*ones(1,N);
    signal0 = zeros(1,N);

    
else
    perror(sprintf('Unknown signal set %s\n',signal_set));
end
% 
 x = [];    % 信号波形的数组
 y = [];
% 
% 
%
    m=new_bits_two(:,536:536);
    for n=1:length(m)
        %%%%   for n=1:length(simple_test)
        %if bits(n)== 0
        
        if m(n)== 0
            %%%%       if simple_test(n)== 0
            x=[x signal0];
            %x = zeros(1,100);
        else x=[x signal1];
            %else x = ones(1,100);
        end
            
                                                                 
    end
  
    %%  波形添加噪声
    r=noisegen(x,SNR_index);                    % 加入指定信噪比的高斯白噪声
    %r = x + noise_amp*randn(1,length(x));     % 加性方式加噪声，混合信号r
    
    
    %%  随机共振参数
    h=1/fs;              % 时间步长
    % a_move=1;            % 可根据所考察参数，平移输出图形，使其更容易理解
    % b_move=1;
    %for M = 0:0.01:4
    % for j=1:1:10
    %
    %     r = x+ j*randn(1,length(x));
    %
    % for a =0.1:0.1:2       % 遍历参数a,b的取值，求误码率
    %     for b = 0.1:0.1:2
    for a =1:1:5      % 遍历参数a,b的取值，求误码率
        for b = 1:1:5
            %         for M = 1:1:10
            
            %r = x+ sqrt(2*M).*randn(1,length(x));
            %             r = x+ sqrt(2*M).*randn(1,length(x));
            %
            %             data{M} = r;
            c = fix(a);
            d = fix(b);
            %             e = fix(M);
            
            
            %         aa=a+a_move;
            %         bb=b+b_move;
            %         a=a./1;      % 参数的考察精度
            %         b=b./1;      % 参数的考察精度
            
            r1=sr(a,b,h,r);     % 随机共振处理
            %             sr_data{c,d,e} = r1;
            sr_data{c,d} = r1;
            %         a=a.*1;      % 参数的考察精度
            %         b=b.*1;      % 参数的考察精度
            
            r2=-r1;    % 码元判别的中间参数
            
            %for n=1:length(bits)                  %累加，用和做判决
            for n=1:length(m)
                %%%%   for n=1:length(simple_test)
                sum1=0; sum2=0;
                
                for k=N/4:N/1              % 积分判别码元，积分的范围
                    sum1=sum1+r1((n-1)*N+k);
                    sum2=sum2+r2((n-1)*N+k);
                end
                if sum1>= sum2
                    output(1,n)=1;
                    
                else
                    output(1,n)=0;
                end
            end
            %         save output.dat output;
            %  求出判决的值，输出为数列，便于比较计算误码率
            c = fix(a);
            d = fix(b);
            %             e = fix(M);
            error_num=0;
            %for i=1:length(bits)
            for i=1:length(m)
                %%%% for i=1:length(simple_test)
                %if bits(i)~=output(i)
                if m(i)~=output(i)
                    %%%%     if simple_test(i)~=output(i)
                    error_num=error_num+1;      % 统计误码个数
                end
            end
            %error_rate=error_num/length(bits);   % 计算误码率
            error_rate=error_num/length(m);
            %             error_ra(c,d,e)=error_rate;
            error_ra(c,d)=error_rate;
            %error_ra(a,b,j)=error_rate;
        end
    end
    
    
    %% 输出最佳参数及相应的最小误码率
    %BER_min=min(min(error_ra));          %输出最小误码率
    %BER_min=min(min(min(error_ra)));
    BER_min = min(error_ra(:));
    % s=size(error_ra);%计算三维维数组的大小
    % Lin=find(x<=BER_min);%计算最小值位置的单下标
    [min_val, position_min] = min(error_ra(:));
    % [c,d,e] = ind2sub(size(error_ra),position_min);
    [c,d] = ind2sub(size(error_ra),position_min);
    %Lax=find(x>=max_x);%计算最大值位置的单下标
    
    %[j,a,b]=ind2sub(s,Lin);%将最小值单下标转为三维多下标
    
    %[m,n,p]=ind2sub(s,Lax);%将最大值单下标转为三维多下标
    
    %Loc_in=[j,a,b];%最小值位置下标
    
    
    %[best_jj best_aa, best_bb]=find(error_ra==BER_min);      %定位最小误码率时的参数取值%[best_cc, best_dd, best_ee]=find(error_ra==BER_min);
    %[best_ee, best_cc, best_dd]=find(error_ra==BER_min);
    %[BER_min, position, index] = min(min(min(error_ra)));
    % [~ , ~]=find(error_ra==BER_min);
    % [best_aa, ~]=find(error_ra==BER_min);
    % [best_jj,  best_bb]=find(error_ra==BER_min);
    % best_a = a; %- a_move;                   %最小误码率时的参数a
    % best_b = b; %- b_move;                   %最小误码率时的参数b
    % best_j = j;
    best_c = c; %- a_move;                   %最小误码率时的参数a
    best_d = d; %- b_move;                   %最小误码率时的参数b
    % best_e = e;
    
    best_a = c;
    best_b = d;
    % best_M = e;
    % r_d = x + sqrt(2*e).*randn(1,length(x));
    % r1_d = sr(best_a,best_b,h,r_d);
    % r2_d = -r1_d;
    %            %for n=1:length(bits)                  %累加，用和做判决
    %             for n=1:length(m)
    %           %%%%   for n=1:length(simple_test)
    %                 sum1=0; sum2=0;
    %
    %                 for k=N/4:N/1              % 积分判别码元，积分的范围
    %                     sum1=sum1+r1_d((n-1)*N+k);
    %                     sum2=sum2+r2_d((n-1)*N+k);
    %                 end
    %                 if sum1>= sum2
    %                     output(n)=1;
    %
    %                 else
    %                     output(n)=0;
    %                 end
    %             end
    %             error_num1=0;
    %             %for i=1:length(bits)
    %             for i=1:length(m)
    %            %%%% for i=1:length(simple_test)
    %                 %if bits(i)~=output(i)
    %                 if x22(i,10)~=output(i)
    %            %%%%     if simple_test(i)~=output(i)
    %                     error_num1=error_num1+1;      % 统计误码个数
    %                 end
    %             end
    %            error_rate1=error_num1/length(m);
    % best_parameters_one = [best_j best_a]
    % best_parameters_two = [best_a best_b]
    % best_parameters_three =[best_j best_b]
    best_parameters = [best_a best_b]
    % best_parameters = [best_a best_b best_M]
    %best_parameters =[BER_min, position, index]
    %BER_min      % 输出最小误码率
    % best_r = data(1,best_M);
    % best_r1 = sr_data{best_c,best_d,best_M};
    best_r1 = sr_data{best_c,best_d};
    peak=max(abs(best_r1));
    rms=sqrt(sum(best_r1.^2)/length(best_r1));
    averageAmplitude=mean(abs(best_r1));
    
    %E = sum(r)/100000;%均值
    % E = mean(best_r{1,1});
    E = mean(best_r1);
    %D = (sum((r-E).*(r-E)))/100000;%方差
    
    S = sum(r);
    % D = var(best_r{1,1});
    D = var(best_r1);
    %P =(sum(r.*r.*r))/100000;%偏度
    % P = skewness(best_r{1,1});%偏度
    P = skewness(best_r1);%偏度
    best_a ;
    best_b;
    % best_M ;
    SNR = SNR_calculate(x,best_r1); %信噪比
    BER_min %误码率
    %C = (max(r1)-min(r1))/(sqrt(sum(r1.*r1)));%波峰指数
    %C_C = C*100000;
    C= peak/rms;%波峰指数
    %I = (max(r1)-min(r1))/sum(abs(r1));%脉冲指数
    %I_I = I*100000;
    I = peak/averageAmplitude;%脉冲指数
    D_net = [E D P best_a best_b SNR_index SNR BER_min C I];
    
    D_net = D_net';
    % D_nett(j) = D_net;
    % figure(2)
    % plot(best_r1)
%     save('dnet535.mat','D_net');


%% 输出图像
% figure(1) 
% surf(error_ra,'EdgeColor','None');%绘制误码率的3D图
% axis ([0 25 0 20 0 1])
% xlim([1 10]);
% colorbar;
% shading interp;
% title('Effects of SR parameters on BER');    %标题
% x1=xlabel('Parameter b');        %x轴标题
% x2=ylabel('Parameter a');        %y轴标题
% x3=zlabel('BER');        %z轴标题

%plot(D_net)



%%
% 如果参数的考察精度小于1，则输出的图像横纵坐标并不是随机共振参数，则可以下方的标记方式做标记
% title('随机共振参数与信号误比特率的关系');    %标题
% x1=xlabel('随机共振 参数b，减1为实际参数');        %x轴标题
% x2=ylabel('随机共振 参数a，减1为实际参数');        %y轴标题
% x3=zlabel('误比特率');        %z轴标题

%% 输出图像
xlim_max=10000;
ylim_max=2;

%% 输出图像
t=[0:length(r)-1];
figure(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(311)    % 原始波形
plot(t,x); axis([0 xlim_max -ylim_max ylim_max]); hold on;
% plot(t,r1,'k');   xlim([0 xlim_max]);ylim([-ylim_max ylim_max]);
title('Original Signal');
subplot(312)
plot(t,r); axis([0 xlim_max -5 5]); hold on;
subplot(313)
plot(t,best_r1); axis([0 xlim_max -ylim_max ylim_max]); hold on;
%% 计算信噪比增益
befor_SNR = SNR_calculate(x,r)
after_SNR = SNR_calculate(x,best_r1)
SNR_gain= after_SNR - befor_SNR