文章目录

前言程序介绍1.模糊函数计算2.模拟退火算法3.代价函数4.结果图象 总结

前言

本文仿真了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_indexcase 1% 1. pcwx = exp(1i*2*pi*fc*t);case 2% 2. lfmk = B/T;f0 = fc-B/2;x = exp(1i*2*pi*(f0*t+k/2*t.^2));case 3% 3. hfmf0 = 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) == 1s((i-1)*length(tbark)+1:i*length(tbark))=exp(1j*2*pi*fc*tbark);elses((i-1)*length(tbark)+1:i*length(tbark))=exp(1j*(2*pi*fc*tbark+pi));endendx = [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);endx = [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:Ns((i-1)*length(tsubpulse)+1:i*length(tsubpulse)) = exp(1j*(2*pi*fc*tsubpulse+pi*k*tsubpulse.^2+route_best(i)/M*2*pi));endx = [s,zeros(1,length(t)-length(s))];endre_fs = 0.9*fs:2:1.1*fs;alpha = re_fs/fs; % Doppler ratio, alpha = 1-2*v/cdoppler = (1-alpha)*fc;% Doppler = 2v/c*fc = (1-alpha)*fcN_a = length( resample(x,fs,min(re_fs)) );N = N_a + length(x)-1;afmag = zeros(length(alpha),N);ticparfor i = 1:length(alpha)if ceil(length(x)/alpha(i)) ~= N_ax_alpha = [resample(x,fs,re_fs(i)),zeros(1, N_a-ceil(length(x)/alpha(i)))];elsex_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));endtocdelay = ([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)))

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 > resEnergy1=Energy_best; % 用于控制循环的结束条件for i = 1: Markov_length% 产生新解(对当前解添加扰动)if rand >0.5% 两点交换a = 0;b = 0;while (a==b)a = randi(N);b = randi(N);endroute_new(a) = randi(M);route_new(b) = randi(M);elseroute_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_bestEnergy_best=Energy_new;route_best = route_new;endelseif rand<exp(-(Energy_new-Energy_current)/temperature)Energy_current = Energy_new;route_current = route_new;elseroute_new = route_current; % 否则路线不更新，保存更改之前的路线endendendhold on;stem(route_new,'b')drawnow;F=getframe(gcf);I=frame2im(F);[I,map]=rgb2ind(I,256);if pic_num == 1imwrite(I,map,'TSP.gif','gif', 'Loopcount',inf,'DelayTime',0.2);elseimwrite(I,map,'TSP.gif','gif','WriteMode','append','DelayTime',0.2);endclfhold offtitle(['SA Algorithm for TSP problem: iteration = ',num2str(pic_num)])xlabel('N')%ylabel('')pic_num = pic_num + 1;if Energy1==Energy_best&&Energy_current==Energy_bestbreakelsetemperature = temperature*ratio; %降温过程endend

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

4.结果图象

总结

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