部分声纳波形模糊函数MATLAB仿真_模糊函数的matlab代码

---

前言

本文仿真了CW、LFM、HFM、Bark二相编码、Costas频率编码信号的模糊函数，此外还利用模拟退火算法寻找N相编码的最优编码并绘制了其模糊函数。
模糊函数的计算和模拟退火算法有参考网络资源。

---

程序介绍

1.模糊函数计算

代码如下：

close all; clear all; clc;

T = 0.26;
fs = 20e3;
fc = 1e3;
B = 1000;
t = 0:1/fs:T-1/fs;

type_index = 6; 
type_name = {'PCW', 'LFM', 'HFM', 'Bark', 'Costas', 'PolyPhase'};

switch type_index
    case 1
        % 1. pcw
        x = exp(1i*2*pi*fc*t);
        
    case 2
        % 2. lfm
        k = B/T;
        f0 = fc-B/2;
        x = exp(1i*2*pi*(f0*t+k/2*t.^2));
        
    case 3
        % 3. hfm
        f0 = fc+B/2;
        beta = B/f0/(fc-B/2)/T;
        x = exp(1i*2*pi/beta*log(1+beta*f0*t));
        
    case 4
        % 4.bark信号
        % 4. hfm + pcw
        %     f0 = fc+B/2;
        %     beta = B/f0/(fc-B/2)/T;
        %     x = sin(2*pi/beta*log(1+beta*f0*t)) + sin(2*pi*fc*t);
%         bark = [1,1,1,1,1,-1,-1,1,1,-1,1,-1,1];
        bark = (randi(2,1,100)-1)*2-1;
        Tbark = T/length(bark);
        tbark = 0:1/fs:Tbark-1/fs;
        s = zeros(1,length(bark)*length(tbark));
        for i = 1:length(bark)
            if bark(i) == 1
                s((i-1)*length(tbark)+1:i*length(tbark))=exp(1j*2*pi*fc*tbark);
            else
                s((i-1)*length(tbark)+1:i*length(tbark))=exp(1j*(2*pi*fc*tbark+pi));
            end
        end
        x = [s, zeros(1,length(t)-length(s))];
        
    case 5
        % 5.costas信号
        costas = [2,4,8,5,10,9,7,3,6,1];
        f = fc-B/2+(costas-1)*B/(length(costas)-1);
        Tcostas = T/length(costas);
        tcostas = 0:1/fs:Tcostas-1/fs;
        s = zeros(1,length(costas)*length(tcostas));
        for i = 1:length(costas)
            s((i-1)*length(tcostas)+1:i*length(tcostas)) = exp(1j*2*pi*f(i)*tcostas);
        end
        x = [s,zeros(1,length(t)-length(s))];
        
    case 6
        % 6.PolyPhase信号
        N = 20;
        M = 2;
        [route_best,Energy_best] = SA_TSP(N, M);
        Tsubpulse = T/N;
        k = B/Tsubpulse;
        tsubpulse = 0:1/fs: Tsubpulse-1/fs;
        s = zeros(1,N*length(tsubpulse));
        for i=1:N
            s((i-1)*length(tsubpulse)+1:i*length(tsubpulse)) = exp(1j*(2*pi*fc*tsubpulse+pi*k*tsubpulse.^2+route_best(i)/M*2*pi));
        end
        x = [s,zeros(1,length(t)-length(s))];
end

re_fs = 0.9*fs:2:1.1*fs;
alpha = re_fs/fs;       % Doppler ratio, alpha = 1-2*v/c
doppler = (1-alpha)*fc;	% Doppler = 2v/c*fc = (1-alpha)*fc

N_a = length( resample(x,fs,min(re_fs)) );
N = N_a + length(x)-1;
afmag = zeros(length(alpha),N);

tic
parfor i = 1:length(alpha)
    if ceil(length(x)/alpha(i)) ~= N_a
        x_alpha = [resample(x,fs,re_fs(i)),zeros(1, N_a-ceil(length(x)/alpha(i)))];
    else
        x_alpha = resample(x,fs,re_fs(i));
    end
    %     x_alpha = Doppler(x,1/re_fs(i),1/fs);
    x_temp = zeros(1,length(x_alpha));
    for j = 1:length(x_alpha)
        x_temp(j) = conj(x_alpha(length(x_alpha)-j+1));
    end
    %     disp(num2str(i))
    afmag_temp = conv(x_temp,x);
    M = length(afmag_temp);
    
    afmag(i,:) = afmag_temp*sqrt(alpha(i));
    
end
toc

delay = ([1:N]-N_a)/fs;

tau = 0.2;
fd = 100;
indext = find(delay>=-tau & delay<=tau);
indexf = find(doppler>=-fd & doppler<=fd);
delay1 = delay(indext);
doppler1 = doppler(indexf);

mag = abs(afmag);
mag = mag/max(max(mag));
% mag = 10*log10(mag);
mag1 = mag(:,indext);
mag1 = mag1(indexf,:);
[row,column] = find(mag1<-100);
mag1(row,column)=-60;

figure(1);
mesh(doppler1,delay1,mag1.');
% axis([-100,100,-tau,tau,-50,0]);
colorbar;
set(gcf,'color','w');
xlabel('Doppler (Hz)','FontName','Times New Roman','FontSize',10);
ylabel('Delay (sec)','FontName','Times New Roman','FontSize',10);
zlabel('Level (dB)','FontName','Times New Roman','FontSize',10);
title(['WAF of ',type_name{type_index} ,' Signal'],'FontName','Times New Roman','FontSize',10);

figure(2);
% v=[-3,-3];
% contour(delay1,doppler1,mag1,v,'ShowText','on');grid on;%模糊度图
contour(delay1,doppler1,mag1);grid on;%模糊度图
xlabel('Delay (Sec)','FontName','Times New Roman','FontSize',10);
ylabel('Doppler (Hz)','FontName','Times New Roman','FontSize',10);
title('Contour of AF','FontName','Times New Roman','FontSize',10)

figure(3)
subplot(211);
plot(doppler1,mag1(:,floor(length(indext)/2)),'b')
xlabel('Doppler (Hz)','FontName','Times New Roman','FontSize',10);
ylabel('Amp','FontName','Times New Roman','FontSize',10);
title('Zero Delay','FontName','Times New Roman','FontSize',10)
subplot(212)
plot(delay1,mag1(floor(length(indexf)/2),:),'b')
xlabel('Delay (sec)','FontName','Times New Roman','FontSize',10);
ylabel('Amp','FontName','Times New Roman','FontSize',10);
title('Zero Doppler','FontName','Times New Roman','FontSize',10)

%%
temp=clock;
temp=sum(temp(4:6))*sum(temp(2:3));
temp=round(temp/10);
rand('seed',temp);
%%
plot([0:length(x_alpha)-1]/length(x_alpha)*fs,abs(fft(x_alpha)));
hold on;
plot([0:length(x_alpha2)-1]/length(x_alpha2)*fs,abs(fft(x_alpha2)))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161

2.模拟退火算法

代码如下：

function  [route_best,Energy_best] = SA_TSP(N,M)
Initial_temp = 1000;
res = 1e-3;                             % 最低温限制
ratio = 0.9;                            % 降温参数，控制温度的下降
temperature = Initial_temp;
Markov_length = 20;                     % 改变解的次数
temp=clock;
temp=sum(temp(4:6))*sum(temp(2:3));
temp=round(temp/10);
rand('seed',temp);
route_new = randi(M, 1, N);                % 新产生的解路线
Energy_current = inf;                   % 当前解的能量
Energy_best = inf;                      % 最优解的能量

route_current = route_new;              % 当前解路线
route_best = route_new;                 % 最优解路线
pic_num = 1;

% 外层while循环控制降温过程，内层for循环控制新解的产生。
while temperature > res
    Energy1=Energy_best;                             % 用于控制循环的结束条件
    for i = 1: Markov_length
        
        % 产生新解(对当前解添加扰动)
        if rand >0.5
            % 两点交换
            a = 0;
            b = 0;
            while (a==b)
                a = randi(N);
                b = randi(N);
            end
            route_new(a) = randi(M);
            route_new(b) = randi(M);
        else
            route_new(randperm(N,3)) = randi(M,1,3);
        end
                
        [Energy_new, ~] = SeqCorr(route_new,N,M);
        % 按照Metroplis准则接收新解
        if Energy_new<Energy_current
            % 更新局部最优
            Energy_current = Energy_new;
            route_current = route_new;
            
            % 更新全局最优
            if Energy_new<Energy_best
                Energy_best=Energy_new;
                route_best = route_new;
            end
            
        else
            if rand<exp(-(Energy_new-Energy_current)/temperature)
                Energy_current = Energy_new;
                route_current = route_new;
            else
                route_new = route_current; % 否则路线不更新，保存更改之前的路线
            end
            
        end
        
    end
    
    hold on;
    stem(route_new,'b')
    drawnow;
    F=getframe(gcf);
    I=frame2im(F);
    [I,map]=rgb2ind(I,256);
    if pic_num == 1
        imwrite(I,map,'TSP.gif','gif', 'Loopcount',inf,'DelayTime',0.2);
    else
        imwrite(I,map,'TSP.gif','gif','WriteMode','append','DelayTime',0.2);
    end
    clf
    hold off
    title(['SA Algorithm for TSP problem: iteration = ',num2str(pic_num)])
    xlabel('N')
%     ylabel('')
    pic_num = pic_num + 1;
    
    if Energy1==Energy_best&&Energy_current==Energy_best
        break
    else
        temperature = temperature*ratio; %降温过程
    end
    
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88

3.代价函数

模拟退火算法所优化的目标函数为编码信号的积分旁瓣电平值（ISL）。
代码如下：

function [cost, Correlation] = SeqCorr(route_new, N, M)
code = exp(1i*[1:M]/M*2*pi);
Seq = code(route_new);
Correlation = abs(xcorr(Seq, Seq));
cost = abs(sum(Correlation)-Correlation(N));
end
1
2
3
4
5
6

4.结果图象

---

总结

以上是自己基于理解所做的仿真，若有错误，希望大家能指出并一起交流学习。