stft_timefreq_complexsignal.m

短时傅立叶变换测频

function ff=stft_timefreq_complexSignal(inputdata,fs)
% 利用FFT进行时频分析，即短时傅立叶变换的特例
% 适用于复信号的瞬时频率估计，信号频率可以是负的
% fuxiongjun 
% Copyright (c) fuxiongjun 2008-2010
% input:inputdata---row data
%       fs  --sampling rate        

close all
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
IfSimSin = 0;
IfSimChirp = 1;
IfSimBark = 0;
BarkLen = 13;   %巴克码长度

fs = 1.5e9;
Fo = 0e6;     %载频
Wo = 2*pi*Fo;  
Tp = 200e-6;    %脉宽
A = 1;          %脉冲幅度
B = 1000e6;      %chirp带宽
K = B / Tp;     %调频斜率
Ts = 1 / fs;            %采样时间间隔
T = -Tp/2 : Ts : Tp/2;        %时间轴
Len = length(T);        %采样点数量

SNR_input = 10;            % 输入信噪比dB
SNR = 10^(SNR_input/10);    % 信噪比转化为单位1
noisePower = A^2/(SNR);     % 噪声功率
power = 10*log10(noisePower);   % 噪声功率 in dB

%列举各种巴克码
Bark1 = -1;
Bark2 = [-1 1];
Bark3 = [-1 -1 1];
Bark4 = [-1 -1 1 -1];
Bark5 = [-1 -1 -1 1 -1];
Bark7 = [-1 -1 -1 1 1 -1 1];
Bark11 = [-1 -1 -1 1 1 1 -1 1 1 -1 1];
Bark13 = [-1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1];

if IfSimSin == 1
    inputdata = exp( j* (Wo*T)  );
    [a,b] = size(inputdata);
    Noise = wgn(a,b,power,'complex');  % 高斯白噪声信号
    inputdata = inputdata + Noise;
end

if IfSimChirp == 1
    inputdata = exp( j* (Wo*T) + j* pi*K*T.^2  );
    [a,b] = size(inputdata);
    Noise = wgn(a,b,power,'complex');  % 高斯白噪声信号
    inputdata = inputdata + Noise;
end

if IfSimBark == 1
    switch BarkLen
        case 1
            Bark = Bark1;
        case 2
            Bark = Bark2;
        case 3
            Bark = Bark3;
        case 4
            Bark = Bark4;
        case 5
            Bark = Bark5;
        case 7
            Bark = Bark7;
        case 11
            Bark = Bark11;
        case 13
            Bark = Bark13;
        otherwise
            disp('Wrong Barker length!');
            return;
    end % end switch
    
    CodeLen = floor( Len / BarkLen );
    RealSigLen = CodeLen * BarkLen;
    TempSig = repmat(Bark,CodeLen,1);
    TempSig = reshape(TempSig, 1, RealSigLen );
    T1 = (1:RealSigLen)*Ts;
    inputdata = exp( j* Wo*T1) .* TempSig;
    [a,b] = size(inputdata);
    Noise = wgn(a,b,power,'complex');  % 高斯白噪声信号
    inputdata = inputdata + Noise;
    SigSpec = fft(inputdata);
    figure;
    plot( abs((SigSpec)));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


count=length(inputdata);
prompt={'谱估计窗长为:（数据点数）--控制频率分辨率'};
freqtitle='给出谱估计数据窗长度(偶数）';
default_freq={'1024'}; 
freq=inputdlg(prompt,freqtitle,1,default_freq,'on');
if isempty(freq)
   warndlg('您取消了本次操作',freqtitle);
end
window_length=str2num(freq{1});
window_length=pow2(nextpow2(window_length)); %如果用户输入的不是2的幂，则进行如此转化

		point_interval=1000000*1/fs;  %转化成微秒,两个数据间的时间间隔
		stra='频率值的输出时间间隔为:';
		stra2='(输入的点间隔)';
		strb=num2str(point_interval);
		strc='(us)X  ';
		strabc=strcat(stra,stra2,strb,strc);
		prompt={strabc};
		inter_time_title='频率值的输出时间间隔—给出输入点间隔的倍数';
		default_inter_time={'256'};
		inter_time=inputdlg(prompt,inter_time_title,1,default_inter_time);
 		if isempty(inter_time)
   		warndlg('您取消了本次操作',inter_time_title);
		end
		step=str2num(inter_time{1});


%提示用户给出保存文件名及其路径名
[yysname,yyspath]=uiputfile('*.dat','将频率估计结果保存为:');
save_name = fullfile(yyspath,yysname);  


%给出等效时间间隔提示
window_length_time=1000000*window_length/fs;  %单位为us
window_length_string=num2str(window_length_time);  %数据转化为字符串实现数字的显示
str11='谱估计数据窗长的等效时间为(us):';
str12=window_length_string;
resolution_freq=1000/window_length_time; %单位为KHz
resolution_freq_str=num2str(resolution_freq);
str21='频率分辨率为：(KHz)';
str22=resolution_freq_str;
step_time=1000000*step/fs;
step_string=num2str(step_time);
str31='数据窗移动步进的等效时间为(us):';
str32=step_string;
msg1=strcat(str11,str12,str21,str22,str31,str32);
msgbox(msg1);

 
 leng_f=floor((count-1-window_length/2)/step)+1;  %定义输出频率的个数,也是作FFT的循环次数
 if leng_f<=0
   errordlg('待处理的数据量太少，请把谱估计数据窗长度设小点');
end

 
ff=zeros(1,leng_f);    %定义输出频率向量
f=(fs/window_length)*(0:(window_length)-1);
h1=waitbar(0,'正在作时频分析...');
for start_point=0:leng_f-1     %最外层循环
not_use_number=count-start_point*step;     %剩余数据长度   
   if not_use_number<window_length
   ye=inputdata(1+start_point*step:start_point*step+not_use_number);
  else
   ye=inputdata(1+start_point*step:start_point*step+window_length);
  end
 
yyy = fft(ye,window_length);
Pyy = yyy.*conj(yyy)/window_length;
Pyy2=Pyy;
locate_f=find(Pyy2==max(Pyy2));  

addr_f=min(locate_f);


ff(start_point+1)=f(addr_f);
if ff(start_point+1)>fs/2
    ff(start_point+1)=ff(start_point+1)-fs;
end
time(start_point+1)=start_point*step/fs; %反映真实时刻

waitbar(start_point/leng_f,h1);
end
close(h1);

if length(time)==leng_f-1
   ff(leng_f)=[];
end

len_end=min(length(time),length(ff));
time=time(1:len_end-1);  %即使多一样长，也删除最后一次的频率估计，因为最后一次数据量少，不准
ff=ff(1:len_end-1);

%给出输出频率个数的提示
time_length_string=num2str(length(time));  %数据转化为字符串实现数字的显示
str11='输出频率值个数为:';
str12=time_length_string;
msg1=strcat(str11,str12);
msgbox(msg1);   


figure(50);
plot(time,ff);
 xlabel('                                     time(s)','HorizontalAlignment','left');
 ylabel('frequency(Hz)');
title('time-frequency analysis based on STFT');

fo = ff(1);
fh = ff( length( ff ) );
BW = (fh - fo)/1000000;
RF = ff( floor(length( ff )/2) )/1000000;

str11='信噪比:';
SNR_str = num2str(SNR_input);
str12 = SNR_str ;
str13='(dB)  ';
str21='信号载频估计值为:';
RF_str = num2str(RF);
str22 = RF_str;
str23='(MHz)  ';
BW_str = num2str(BW);
str31 = '带宽估计值：';
str32 = BW_str;
str33='(MHz)  ';
msg1=strcat(str11,str12,str13,str21,str22,str23,str31,str32,str33);
msgbox(msg1);


[fidw,message] = fopen(save_name,'wb','l'); 
  if fidw==-1							 
     msgbox('文件保存失败')
  end
   
count_write=fwrite(fidw,ff ,'uint32');  %读出用little_endian读取。（对应这里的fopen，也是默认方式）
fclose(fidw);

%文件保存好后给出提示
str11='频率估计结果已经保存为';
str12=save_name;
msg1=strcat(str11,str12);
msgbox(msg1);