Matlab_spectrogram_短时傅里叶分析_实现与讨论 [未完成]

function [ S , W , T ] = mf_spectrogram...
    ( signal , nsc , nov , fs )
%MF_SPECTROGRAM Short-time Fourier transform realization
% Input         
%       signal      - Signal vector
%       nsc         - Abb. of Number SeCtion
%       nov         - Abb. of Number OverLap
%       fs          - Abb. of Frequency of Sample
% Output
%       S           - A matrix that each colum is a FFT for time of nsc
%       W           - A vector labeling frequency 
%       T           - A vector labeling time

  % Signal Preprocessing
    h = hamming(nsc, ‘periodic‘);   % Hamming weight function
    L = length(signal);             % Length of Signal
    nst = nsc-nov;                  % Number of STep per colum
    ncl = fix( (L-nsc)/nst ) + 1;   % Number of CoLum
    nff = 2^nextpow2(nsc);          % Number of Fast Fourier Transformation
    Ps = zeros(nsc,ncl);
    for n = 1:ncl                   % Ps means Processed Signal
        Ps(:,n) = signal( (n-1)*nst+1 : (n-1)*nst+nsc ).*h‘;
    end                             % Ps is a matrix

  % Short-time Fourier transform
    STFT0 = fft(Ps,nff);

  % Turn into DFT in dB
    STFT1 = abs(STFT0/nff);
    STFT2 = STFT1(1:nff/2+1,:);             % Symmetric results of FFT
    STFT2(2:end-1,:) = 2*STFT2(2:end-1,:);  % Add the value of the other half 
    STFT3 = 20*log10(STFT2);                % Turn sound pressure into dB level

  % Axis Generating
    fax = fs*(0:(nff/2))/nff;                           % Frequency axis setting
    tax = ( .5*nsc : nst : nst*(ncl-1)+.5*nsc ) / fs;   % Time axis generating
    [ffax,ttax] = meshgrid(tax,fax);                    % Generating grid of figure
    
  % Output
    W = ffax;
    T = ttax;
    S = STFT3;
end

function [  ] = my_pcolor( f , t , s )
%MY_PCOLOR 绘图函数
% Input             
%       f               - 频率轴矩阵
%       t               - 时间轴矩阵
%       s               - 幅度轴矩阵
    gca = pcolor(f,t,s);                      % 绘制伪彩色图
        set(gca, ‘LineStyle‘,‘none‘);         % 取消网格，避免一片黑
    handl = colorbar;                         % 彩图坐标尺
        set(handl, ‘FontName‘, ‘Times New Roman‘, ‘FontSize‘, 14)
        ylabel(handl, ‘Magnitude, dB‘)        % 坐标尺名称
end

clc 
clear
close all

% 基本参数
fa = [ 50 800 ];    % 扫描频率上下限
fs = 6400;          % 采样频率

% 分辨率相关设定参数
te = 1;             % [s] 扫频时间长度
fre = 8;            % [s] 频率分辨率
tre = 0.002;        % [Hz] 时间分辨率

% Chirp 信号生成
t = 0:1/fs:te;                  % [s] 时间序列
sc = chirp(t,fa(1),te,fa(2));   % 信号生成

% 分辨率相关输入参数
nsc = floor(fs/fre);
% nov = floor(nsc-(fs*tre));
nov = floor(nsc*0.9);
nff = max(256,2^nextpow2(nsc));

% 计算与绘制结果
subplot(1,3,1)                            % 绘制自编函数结果
    [S,W,T] = mf_spectrogram(sc,nsc,nov,fs);
    my_pcolor(W,T,S)
    caxis([-130.86,-13.667]);
    title(‘自编函数‘);
    xlabel(‘时间 second‘);
    ylabel(‘频率 Hz‘);
subplot(1,3,2)                            % 绘制 Matlab 函数结果
    s = spectrogram(sc,hamming(nsc),nov,nff,fs,‘yaxis‘);
        % Turn into DFT in dB
        s1 = abs(s/nff);
        s2 = s1(1:nff/2+1,:);             % Symmetric results of FFT
        s2(2:end-1,:) = 2*s2(2:end-1,:);  % Add the value of the other half 
        s3 = 20*log10(s2);                % Turn sound pressure into dB level
    my_pcolor(W,T,s3 )
    caxis([-130.86,-13.667]);
    title(‘Matlab 自带函数‘);
    xlabel(‘时间 second‘);
    ylabel(‘频率 Hz‘);
subplot(1,3,3)                            % 两者误差
    my_pcolor(W,T,20*log10(abs(10.^(s3/20)-10.^(S/20))))
    caxis([-180,-13.667]);
    title(‘error‘);
    ylabel(‘频率 Hz‘);
    xlabel(‘时间 second‘);
suptitle(‘Spectrogram 实现与比较‘);