clc
clear all
close all
%% 产生信号(跳频信号)
%频率集
f1 = 25e3;f2 = 100e3; f3= 125e3; f4= 75e3; f5= 175e3; f6= 150e3; f7= 200e3; f8= 50e3;
fs = 500e3; %采样频率
fh= 2000 ; %跳速
Th= 1/fh; %跳周期
t1 = 0: 1/fs : Th- 1/fs; %第一跳周期内的采样时间点
t=0:1/fs:8*Th-1/fs; %整段信号的采样时间点
x1 = cos(2*pi*f1*t1); x2 = cos(2*pi*f2*t1); x3 = -1*cos(2*pi*f3*t1); x4 = -1*cos(2*pi*f4*t1);
x5 = cos(2*pi*f5*t1); x6 = cos(2*pi*f6*t1); x7 = -1*cos(2*pi*f7*t1); x8 = -1*cos(2*pi*f8*t1);
s= [x1 x2 x3 x4 x5 x6 x7 x8 ]; %整段接收信号
x=s;
x=awgn(x,5,'measured');%加高斯白噪声，信噪比为5dB
figure(1);
plot(t,x); %接收信号的时域图
x = x';
x = hilbert(x);
N = length(x);
f = 0 :fs/length(x):(fs/2);
y = fft(x);
figure(2);
xlabel('频率 Hz');
ylabel('幅度');
plot(f,abs(y(1:N/2+1))) %接收信号的频域图
xlabel('Time/s'),ylabel('Amplitude');
title('frequency spectrum')
grid on;
%% STFT
H = window(@hamming,201); %加窗（汉明窗，@hamming。127个点）
[TFR1,T,F]=tfrstft(x,1:length(x),N,H); %短时傅里叶变换函数tfrstft
%输入为：信号x
% 信号长度
% 所加窗函数
%输出为：短时傅里叶变换得到的时频矩阵
% 时频矩阵对应的时间序列
% 时频矩阵对应的频率序列
figure(3);
mesh(T(1:N),F(1:N/2),abs(TFR1(1:N/2,1:N))); %画三维图
xlabel('Time'),ylabel('Frequency'),zlabel('Amplitude');
grid on;
figure(4);
contour(T,F,abs(TFR1)); %画等高线
xlabel('Time'),ylabel('Normalized frequency');
axis([0 N 0 0.5])
figure(5);
imagesc(abs(TFR1)); %imagesc(A) 将矩阵A中的元素数值按大小转化为不同颜色，
%并在坐标轴对应位置处以这种颜色染色
title('STFT imagesc 汉明窗（127）');
figure(10);
[fmax,ff]=max(abs(TFR1));%%%%%%%%%%%% 时频脊线
plot(ff/N)
xlabel('时间t');
ylabel('频率f');
title('STFT变换时频脊线')
ylim([0 0.5])
%% 提取时间/频率平面
stft_t=max(abs(TFR1),[],1); %每个时间点所有频率对应的最大频谱幅值
stft_t=stft_t./max(stft_t); %幅度归一化
stft_t=stft_t-mean(abs(stft_t)); %去直流（mean取平均值）
figure(6);
plot((0:length(stft_t)-1)*max(t)/length(stft_t),stft_t./max(stft_t));
title('Relevance of STFT and Time')
xlabel('Time/s'),ylabel('Normalized |STFT(t)|')
y=fftshift(abs(fft(stft_t)));
f=(0:(length(y)-1))*fs/length(y) - 0.5*fs;
index=find(y==max(y));
fhd=(index-1)*fs/length(y) - 0.5*fs;%频率是从0开始的，数组下标从1开始
fhdd=fhd(find(fhd>0));
Thd=1/fhdd;
figure(7);
plot(f,y);
xlabel('Frequency/Hz'),ylabel('Amplitude');
serch=x(1:Thd*N/max(t));
xd1=fftshift(abs(fft(serch)));
figure;plot((0:length(xd1)-1)*fs/length(xd1),abs(xd1));
fd1=find(xd1==max(xd1))*fs/length(xd1);
f1_erro=abs(fd1-f1)/f1;
% stft_f=max(abs(TFR1),[],2); %每个频率点在所有时间内对应的最大频谱幅值
% figure;
% % f1=((0:1/2*(length(stft_f)-1))*0.5*fs/(1/2*(length(stft_f))));
% % plot(f1,stft_f(1:1/2*(length(stft_f)))./max(stft_f(1:1/2*(length(stft_f)))));
% f1=((0:(length(stft_f)-1))*fs/((length(stft_f))));
% plot(0.5*f1,stft_f(1:(length(stft_f)))./max(stft_f(1:(length(stft_f)))));
% %
% title('Relevance of STFT and Frequency')
% xlabel('Frequency/Hz'),ylabel('Normalized |STFT(f)|')
%% 参数估计
% %跳周期估计
%
% I=mean(stft_t);
% I_max=find(stft_t>0.79*I);
% I_min=find(stft_t<=0.79*I);
% stft_t(I_max)=0;
% stft_t(I_min)=1;
% figure;
% plot((0:length(stft_t)-1)*max(t)/length(stft_t),stft_t);
% stft_t1=stft_t(101:end);
% IndMin=find(diff(sign(diff(stft_t)))>0)+1; %寻找极小值
% figure;
% plot(stft_t);
% title('Relevance of STFT and Time')
% xlabel('Time'),ylabel('Normalized |STFT(t)|')
% hold on;
% plot(IndMin,stft_t(IndMin),'r^')
%
% max_t=max(t)-max(t1);
% a=min(stft_t1(IndMin));
% I=find(((stft_t1(IndMin)-a)./a)<0.6);%寻找波谷
% ht=IndMin(I)*max_t/length(stft_t1);%波谷对应的时刻为跳时刻
% Th_t=mean(abs(diff(ht)));
% fh_t=1/Th_t;
% eff=abs(Th_t-Th)/Th
% % eff2=abs(fh_t-fh)/fh
% A=find(stft_t1(IndMin)==a);
% T1=mod(A*max_t/length(stft_t1)+max(t1),Th_t)
% T1_err=abs(T1-max(t1))/(max(t1))
% % figure;
% % plot(((0:(length(stft_t)-1))*fs/(length(stft_t))),fftshift(abs(fft(stft_t))));
% % L=N/8-1; %STFT滑窗长度
% % h=mod(L*sum(IndMin(I)),fs/fh_t)/(fs/fh_t);
% % p=length(ht);
% % h=(sum(ht)-(p-1)*p*fs*Th_t/2)/p-Th*fs;
% %------------跳频频率估计----------------%
% serch=x(1:ht(1)*N/max(t));
% % serch=x(1:Th*N/max(t));
% % s1=pmusic(serch,1);
% s1=fft(serch);
% figure;
% plot((0:length(s1)-1)*fs/length(s1),abs(s1));
% f_t(1)=find(abs(s1)==max(abs(s1)))*fs/length(s1);
% f_err(1)=abs(f_t(1)-f1)/f1;
%
% for ii=1:length(ht)
% if ii==length(ht)
% serch=x(ht(ii)*N/max(t):end);
% else
% serch=x(ht(ii)*N/max(t):round((ht(ii)+Th_t)*N/max(t)));%把跳频信号按时序分成单频信号
% end
% s=fft(serch);
% f_t(ii+1)=find(abs(s)==max(abs(s)))*fs/length(s);
% f_t(ii+1)=f_czt(serch,f_t(ii+1),fs);%czt频谱细化
% f_err(ii+1)=abs(f_t(ii+1)-f(ii+1))/f(ii+1);
% end
评论10