【限免】短时傅里叶变换时频分析【附MATLAB代码】

来源：微信公众号：EW Frontier

简介

一种能够同时对信号时域和频域分析的方法——短时傅里叶变换（STFT），可以在时频二维角度准确地描述信号 的时间、频域的局部特性，与其他算法不同，通过该算法可以得到信号的瞬时频率随时间变化的变化规律，其在雷达信号的脉内 特征分析中的效果明显。本文根据仿真结果，对不同类型信号经短时傅里叶变换后的结果进行统计，形成了基于短时傅里叶变换 的雷达信号脉内特征自动识别流程，对电子侦察情报的获取及应用有重要的意义。

电子侦察面对的电磁环境越来越复杂，这意味着雷达 信号越来越复杂，也就是说雷达信号环境中常规脉冲雷达信 号使用的比例逐渐减小。因此，对线性调频、相位编码、非 线性调频、频率捷变等雷达信号的识别已成为不可忽视重要 环节，同时，高速地对这些信号进行自动识别也是发展的需 求。本文利用短时傅里叶变换对各种复杂信号进行脉内特征 分析，以探索脉内特征自动识别流程。

STFT主要原理

短时傅里叶变换公式为：

式（1）中，m（τ-t）为窗函数在时域窗函数取截信号， 窗函数滤波出来的局部时域数据作FFT，就是在τ时刻得 到该时域窗函数对应信号的傅里叶变换，设置不同的τ值， 窗函数的中心位置会不断移动，这样就得到了不同τ下该信 号的傅里叶变换，这些不同时刻傅里叶变换的集合就是Sx （ω，t）。

从式（1）可以发现，从形式上来看，傅里叶变换和短时傅里叶变换的区别仅在于时间域上少了一个窗函 数。短时傅里叶变换通过采用滑窗处理来弥补传统傅里叶 变换的不足之处。它能够对某一小段时间滑窗内的信号做 傅里叶变换，反映该信号的频域随时间变换的大致规律。 同时，短时傅里叶变换还保留了传统傅里叶变换较好的抗 干扰能力。仿真研究表明，若时域的滑窗时间越短，则变 换后的频率分辨率会越低；若滑窗时间延长，则时域的分 辨率就会降低，这是短时傅里叶变换在时域分辨率、频域 分辨率方面存在的矛盾。

MATLAB代码

function [t, f, S] = stft1(x,N,M,Nfft,Fs,win_type)
%STFT1 Short-Time Fourier Transform (STFT) - Method I.
%
%   [t, f, S] = stft1(x,N,M,Nfft,Fs,win_type) calculates the Short-Time Fourier Transform 
%   (STFT) or in other words the spectrogram of the input signal x.
%
%   Inputs:
%            x is a row vector that contains the signal to be examined.
%            N is the selected length of the window in samples.
%            M is the selected amount of overlap between successive windows in samples.
%            Nfft is the selected number of DFT calculation points (Nfft>=N). 
%            Fs is the signal sampling frequency in Hz.
%            win_type is a string containing one of the windows shown
%                      below. The default value of this variable corresponds to the rectangular window.
%
% Outputs: S is a matrix that contains the complex spectrogram of x.    
%          i.e. S(:,k) = fft(k-th frame) computed on Nfft points.
%          f is a vector of frequency points in Hz where the spectrogram is calculated.
%          t is a vector of time points in sec. Each point corresponds to a
%            signal frame.
%
%   Copyright 2020-2030, Ilias S. Konsoulas.
 
%% Window selection and contruction.
switch win_type 
    
    case  'cheby'
         win = chebwin(N).';
 
    case 'blackman'
         win = blackman(N).';
 
    case 'hamm'
         win = hamming(N).';
        
    case 'hann'
         win = hanning(N).';  
        
    case 'kaiser'
         beta = 5;
         win = kaiser(N,beta).';  
    
    case 'gauss'
         win = gausswin(N).';
         
    otherwise  % otherwise use the rectangular window
         win = ones(1,N);
end
 
%% Input Signal Segmentation Params.
x = x(:).';
L = length(x);
% Number of segments (frames) the signal is divided to.
K = floor((L-M)/(N-M)); 
 
%% DFT Matrix Construction
n = 0:Nfft-1;
w = exp(-1i*2*pi*n/Nfft);
W = zeros(Nfft,Nfft);
for k=0:Nfft-1
    W(:,k+1) = w.^k;
end
 
%% Number of Unique FFT Points.
NUPs = Nfft;   
if isreal(x)
   if mod(Nfft,2)   % if N is odd.
        NUPs = (Nfft+1)/2;
   else             % if N is even.
        NUPs = Nfft/2+1; 
   end
end 
 
%% Input Signal Segmentation (Framing)
 
X = zeros(N,K);
 
for k=1:K
    X(:,k) = x((k-1)*(N-M)+1:k*N - (k-1)*M).*win;
end
 
X = vertcat(X,zeros(Nfft-N,K));  % Zero-Padding each frame. 
 
%% DFT Calculation of all frames as a single matrix product.
S = W*X;  
 
S = S(1:NUPs,:);
 
%% Frame Time Points
t = (N-1)/Fs + (0:K-1)*(N-M)/Fs; % Frame Time Points.
 
%% Frequency Points in Hz.
f = (0:NUPs-1)*Fs/Nfft;  
 
%% NOTES:
% When K is an integer then the following equation holds:
% (N-1)*Ts + (K-1)*(N-M)*Ts = (L-1)*Ts = total signal duration in sec.
% or (N-1) + (K-1)*(N-M) = (L-1).

% Short-Time Fourier Transform - Method I (Spectrogram) Demo.
% This demo shows how the custom function stft1.m and is used to
% to produce the spectrogram of an input signal.
 
%   Copyright 2020 - 2030, Ilias S. Konsoulas.
 
%% Workspace Initialization.
 
clc; clear; close all;
 
%% Select and load the signal to be analyzed.
% load('chirp','Fs','y'); x = y;
% load('gong', 'Fs','y'); x = y;
% load timit2.asc;  Fs = 8000; x = timit2; 
% load('train','Fs','y'); x = y;
load('splat','Fs','y'); x = y;
% load('laughter','Fs','y'); x = y;
%  [x Fs] = audioread('andean-flute.wav');
 
%% Play Back the selected audio signal:
soundsc(x,Fs,24);
 
%% Signal Normalization.
x = x.'/max(abs(x));  
 
%% STFT Parameters.
L    = length(x);
N    = 512;   % Selected window size.
M    = 450;   % Selected overlap between successive segments in samples.
Nfft = 512;   % Selected number of FFT points.
 
[t,f,S] = stft1(x,N,M,Nfft,Fs,'hamm');
 
%% Plot the Spectrogram.
h = figure('Name','STFT - Method I Demo');
colormap('jet');
 
[T,F] = meshgrid(t,f/1000); % f in KHz.
surface(T,F,10*log10(abs(S.^2) + eps),'EdgeColor','none');
 
axis tight;
grid on;
title(['Signal Length: ',num2str(L),', Window Length: ', num2str(N),', Overlap: ', num2str(M), ' samples.']);
xlabel('Time (sec)');
ylabel('Frequency (KHz)');
colorbar('Limits',[-80, 40]);
cbar_handle = findobj(h,'tag','Colorbar');
set(get(cbar_handle,'YLabel'),'String','(dB)','Rotation',0);
zlim([-80 40]);

MATLAB仿真结果