参考文章：http://www.cnblogs.com/adgk07/p/9314892.html
首先根据他这个代码和我之前手上已经拥有的那个代码，编写了一个适合自己的代码。
首先模仿他的代码，测试成功。
思路：
短时傅里叶变换，其实还是傅里叶变换，只不过把一段长信号按信号长度（nsc）、重叠点数（nov）重新采样。
% 结合之前两个版本的stft，实现自己的周期图法，力求通俗易懂，代码分明。
% 该代码写的时候是按照输入信号为实数的思路写的，在每个片段fft时进行前一半行的转置存储。后续代码思路已修改。 clear; clc; close all % %---------------------Chirp 信号生成(start)--------------------- fa = [ 50 800 ]; % 扫描频率上下限 fs = 6400; % 采样频率 % 分辨率相关设定参数 te = 1; % [s] 扫频时间长度 fre = 8; % [s] 频率分辨率 tre = 0.002; % [Hz] 时间分辨率 t = 0:1/fs:te; % [s] 时间序列 sc = chirp(t,fa(1),te,fa(2)); % 信号生成 %待传参(实参) signal = sc; nsc = floor(fs/fre); nff = 2^nextpow2(nsc);%每个窗口进行fft的长度 nov = floor(nsc*0.9); % %---------------------Chirp 信号生成(end)--------------------- %% 使用fft实现周期图法 % Input % signal - Signal vector 输入信号(行向量) % nsc - Number SeCtion 每个小分段的信号长度 % nov - Number OverLap 重叠点数 % fs - Frequency of Sample 采样率 % Output % S - A matrix that each colum is a FFT for time of nsc % 如果nfft为偶数，则S的行数为(nfft/2+1)，如果nfft为奇数，则行数为(nfft+1)/2 % 每一列是一个小窗口的FFT结果，因为matlab的FFT结果是对称的，只需要一半 % W - A vector labeling frequency % T - A vector labeling time % Signal Preprocessing h = hamming(nsc, 'periodic'); % Hamming weight function 海明窗加权，窗口大小为nsc L = length(signal); % Length of Signal 信号长度 nstep = nsc-nov; % Number of STep per colum 每个窗口的步进 ncol = fix( (L-nsc)/nstep ) + 1; % Number of CoLum 信号被分成了多少个片段 nfft = 2^nextpow2(nsc); % Number of Fast Fourier Transformation 每个窗口FFT长度 nrow = nfft/2+1; %Symmetric results of FFT STFT_X = zeros(nrow,ncol); %STFT_X is a matrix,初始化最终结果 index=1;%当前片段第一个信号位置在原始信号中的索引 for i=1:ncol %提取当前片段信号值,并用海明窗进行加权(均为长度为nsc的行向量) temp_S=signal(index:index+nsc-1).*h'; %对长度为nsc的片段进行nfft点FFT变换 temp_X=fft(temp_S,nfft); %从长度为nfft点(行向量)的fft变换中取一半(转换为列向量)，存储到矩阵的列向量 STFT_X(:,i)=temp_X(1:nrow)'; %将索引后移 index=index + nstep; end % Turn into DFT in dB STFT1 = abs(STFT_X/nfft); STFT1(2:end-1,:) = 2*STFT1(2:end-1,:); % Add the value of the other half STFT3 = 20*log10(STFT1); % Turn sound pressure into dB level % Axis Generating fax =(0:(nfft/2)) * fs/nfft; % Frequency axis setting tax = (nsc/2 : nstep : nstep*(ncol-1)+nsc/2)/fs; % Time axis generating [ffax,ttax] = meshgrid(tax,fax); % Generating grid of figure % Output W = ffax; T = ttax; S = STFT3; %% 画频谱图 subplot(1,3,1) % 绘制自编函数结果 my_pcolor(W,T,S) caxis([-130.86,-13.667]); title('自编函数'); xlabel('时间 second'); ylabel('频率 Hz'); subplot(1,3,2) % 绘制 Matlab 函数结果 s = spectrogram(signal,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 实现与比较');
内部调用了my_pcolor.m
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', 8) ylabel(handl, 'Magnitude, dB') % 坐标尺名称 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', 8) ylabel(handl, 'Magnitude, dB') % 坐标尺名称 end
接下来在主体不变的情况下，修改输入信号和函数返回结果，使其和自己之前的代码效果相同。
其实只要搞清楚了理论，实现上很简单。
一、首先分析Matlab周期图法API。
[spec_X,spec_f,spec_t]=spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);
如果函数没有返回结果的调用，MATLAB直接帮我们绘图了。
1）当输入信号为复数时的代码：
clear clc close all Fs = 1000; % Sampling frequency T = 1/Fs; % Sampling period L = 1000; % Length of signal t = (0:L-1)*T; % Time vector S = 20*sin(150*2*pi*t)+40*cos(250*2*pi*t)*j; nsc=100;%海明窗的长度,即每个窗口的长度 nov=0;%重叠率 nff= max(256,2^nextpow2(nsc));%N点采样长度 %% matlab自带函数 figure(1) spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs); title('spectrogram函数画图') [spec_X,spec_f,spec_t]=spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);
变量：
其中：
spec_X矩阵行数为nff,列数是根据信号signal的长度和每个片段nsc分割决定的，共计col列。
spec_f 为nff*1的列向量。
spec_t为 1*ncol的行向量。
此时自己实现的Matlab代码为：
1 % 结合之前两个版本的stft，实现自己的周期图法，力求通俗易懂，代码分明。 2 clear; clc; close all 3 %% 信号输入基本参数 4 Fs = 1000; % Sampling frequency 5 T = 1/Fs; % Sampling period 6 L = 1000; % Length of signal 7 t = (0:L-1)*T; % Time vector 8 S = 20*sin(150*2*pi*t)+40*cos(250*2*pi*t)*j; 9 10 11 %传参 12 signal = S; 13 nsc = 100; %每个窗口的长度,也即海明窗的长度 14 nff = max(256,2^nextpow2(nsc));%每个窗口进行fft的长度 15 nov = 0; %重叠率 16 fs = Fs; %采样率 17 18 %% 使用fft实现周期图法 19 %后续可封装为函数: 20 % function [ S , W , T ] = mf_spectrogram... 21 % ( signal , nsc , nov , fs ) 22 % Input 23 % signal - Signal vector 输入信号(行向量) 24 % nsc - Number SeCtion 每个小分段的信号长度 25 % nov - Number OverLap 重叠点数 26 % fs - Frequency of Sample 采样率 27 % Output 28 % S - A matrix that each colum is a FFT for time of nsc 29 % 如果nfft为偶数，则S的行数为(nfft/2+1)，如果nfft为奇数，则行数为(nfft+1)/2 30 % 每一列是一个小窗口的FFT结果，因为matlab的FFT结果是对称的，只需要一半 31 % W - A vector labeling frequency 频率轴 32 % T - A vector labeling time 时间轴 33 % Signal Preprocessing 34 h = hamming(nsc, 'periodic'); % Hamming weight function 海明窗加权，窗口大小为nsc 35 L = length(signal); % Length of Signal 信号长度 36 nstep = nsc-nov; % Number of STep per colum 每个窗口的步进 37 ncol = fix( (L-nsc)/nstep ) + 1; % Number of CoLum 信号被分成了多少个片段 38 nfft = max(256,2^nextpow2(nsc)); % Number of Fast Fourier Transformation 每个窗口FFT长度 39 nrow = nfft/2+1; 40 % 41 % 42 STFT1 = zeros(nfft,ncol); %STFT1 is a matrix 初始化最终结果 43 index = 1;%当前片段第一个信号位置在原始信号中的索引 44 for i=1:ncol 45 %提取当前片段信号值,并用海明窗进行加权(均为长度为nsc的行向量) 46 temp_S=signal(index:index+nsc-1).*h'; 47 %对长度为nsc的片段进行nfft点FFT变换 48 temp_X=fft(temp_S,nfft); 49 %将长度为nfft点(行向量)的fft变换转换为列向量，存储到矩阵的列向量 50 STFT1(:,i)=temp_X'; 51 %将索引后移 52 index=index + nstep; 53 end 54 55 56 %如果是为了和标准周期图法进行误差比较，则后续计算(abs,*2,log(+1e-6))不需要做，只需 57 %plot(abs(spec_X)-abs(STFT1))比较即可 58 59 STFT2 = abs(STFT1/nfft); 60 %nfft是偶数，说明首尾是奈奎斯特点，只需对2:end-1的数据乘以2 61 STFT2(2:end-1,:) = 2*STFT2(2:end-1,:); % Add the value of the other half 62 %STFT3 = 20*log10(STFT1); % Turn sound pressure into dB level 63 STFT3 = 20*log10(STFT2 + 1e-6); % convert amplitude spectrum to dB (min = -120 dB) 64 65 % Axis Generating 66 fax =(0:(nfft-1)) * fs./nfft; % Frequency axis setting 67 tax = (nsc/2 : nstep : nstep*(ncol-1)+nsc/2)/fs; % Time axis generating 68 69 % Output 70 W = fax; 71 T = tax; 72 S = STFT3; 73 74 75 %% matlab自带函数 76 figure(); 77 spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs); 78 title('spectrogram函数画图') 79 [spec_X,spec_f,spec_t]=spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs); 80 %减法，看看差距 81 figure() 82 plot(abs(spec_X)-abs(STFT1)) 83 84 %% 画频谱图 85 figure 86 %利用了SIFT3 87 surf(W, T, S') 88 shading interp 89 axis tight 90 box on 91 view(0, 90) 92 set(gca, 'FontName', 'Times New Roman', 'FontSize', 14) 93 xlabel('Frequency, Hz') 94 ylabel('Time, s') 95 title('Amplitude spectrogram of the signal') 96 handl = colorbar; 97 set(handl, 'FontName', 'Times New Roman', 'FontSize', 14) 98 ylabel(handl, 'Magnitude, dB') 99 100 %跳频图 101 figure(); 102 surf(T,W,S) 103 shading interp 104 axis tight 105 box on 106 view(0, 90) 107 set(gca, 'FontName', 'Times New Roman', 'FontSize', 14) 108 ylabel('Frequency, Hz') 109 xlabel('Time, s') 110 title('Amplitude spectrogram of the signal') 111 handl = colorbar; 112 set(handl, 'FontName', 'Times New Roman', 'FontSize', 14) 113 ylabel(handl, 'Magnitude, dB')
重点关注8行、59行、66行、69行、82行
其实此时只要牢牢记住结果集中行数为nfft就行了。
2）当输入信号为实数时的代码：
clear
clc
close all
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sampling period
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
S = 20*sin(150*2*pi*t)+40*cos(250*2*pi*t);
nsc=100;%海明窗的长度,即每个窗口的长度
nov=0;%重叠率
nff= max(256,2^nextpow2(nsc));%N点采样长度
%% matlab自带函数
figure(1)
spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);
title('spectrogram函数画图')
[spec_X,spec_f,spec_t]=spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);
其中：
spec_X矩阵行数为nff/2+1,列数是根据信号signal的长度和每个片段nsc分割决定的，共计col列。
spec_f 为nff/2+1*1的列向量。
spec_t为 1*ncol的行向量。
主要是因为输入信号为实数时，fft变换的结果有一半是对称的，不必考虑。
自己实现的Matlab代码为：
1 % 结合之前两个版本的stft，实现自己的周期图法，力求通俗易懂，代码分明。 2 clear; clc; close all 3 %% 信号输入基本参数 4 Fs = 1000; % Sampling frequency 5 T = 1/Fs; % Sampling period 6 L = 1000; % Length of signal 7 t = (0:L-1)*T; % Time vector 8 S = 20*sin(150*2*pi*t)+40*cos(250*2*pi*t); 9 10 11 %传参 12 signal = S; 13 nsc = 100; %每个窗口的长度,也即海明窗的长度 14 nff = max(256,2^nextpow2(nsc));%每个窗口进行fft的长度 15 nov = 0; %重叠率 16 fs = Fs; %采样率 17 18 %% 使用fft实现周期图法 19 %后续可封装为函数: 20 % function [ S , W , T ] = mf_spectrogram... 21 % ( signal , nsc , nov , fs ) 22 % Input 23 % signal - Signal vector 输入信号(行向量) 24 % nsc - Number SeCtion 每个小分段的信号长度 25 % nov - Number OverLap 重叠点数 26 % fs - Frequency of Sample 采样率 27 % Output 28 % S - A matrix that each colum is a FFT for time of nsc 29 % 如果nfft为偶数，则S的行数为(nfft/2+1)，如果nfft为奇数，则行数为(nfft+1)/2 30 % 每一列是一个小窗口的FFT结果，因为matlab的FFT结果是对称的，只需要一半 31 % W - A vector labeling frequency 频率轴 32 % T - A vector labeling time 时间轴 33 % Signal Preprocessing 34 h = hamming(nsc, 'periodic'); % Hamming weight function 海明窗加权，窗口大小为nsc 35 L = length(signal); % Length of Signal 信号长度 36 nstep = nsc-nov; % Number of STep per colum 每个窗口的步进 37 ncol = fix( (L-nsc)/nstep ) + 1; % Number of CoLum 信号被分成了多少个片段 38 nfft = max(256,2^nextpow2(nsc)); % Number of Fast Fourier Transformation 每个窗口FFT长度 39 nrow = nfft/2+1; 40 % 41 % 42 STFT1 = zeros(nfft,ncol); %STFT1 is a matrix 初始化最终结果 43 index = 1;%当前片段第一个信号位置在原始信号中的索引 44 for i=1:ncol 45 %提取当前片段信号值,并用海明窗进行加权(均为长度为nsc的行向量) 46 temp_S=signal(index:index+nsc-1).*h'; 47 %对长度为nsc的片段进行nfft点FFT变换 48 temp_X=fft(temp_S,nfft); 49 %将长度为nfft点(行向量)的fft变换转换为列向量，存储到矩阵的列向量 50 STFT1(:,i)=temp_X'; 51 %将索引后移 52 index=index + nstep; 53 end 54 55 % -----当输入信号为实数时，对其的操作（也就是说只考虑信号的前一半） 56 % 由于实数FFT的对称性,提取STFT1的nrow行 57 STFT_X = STFT1(1:nrow,:); % Symmetric results of FFT 58 59 %如果是为了和标准周期图法进行误差比较，则后续计算(abs,*2,log(+1e-6))不需要做，只需 60 %plot(abs(spec_X)-abs(STFT_X))比较即可 61 62 STFT2 = abs(STFT_X/nfft); 63 %nfft是偶数，说明首尾是奈奎斯特点，只需对2:end-1的数据乘以2 64 STFT2(2:end-1,:) = 2*STFT2(2:end-1,:); % Add the value of the other half 65 %STFT3 = 20*log10(STFT1); % Turn sound pressure into dB level 66 STFT3 = 20*log10(STFT2 + 1e-6); % convert amplitude spectrum to dB (min = -120 dB) 67 68 % Axis Generating 69 fax =(0:(nrow-1)) * fs./nfft; % Frequency axis setting 70 tax = (nsc/2 : nstep : nstep*(ncol-1)+nsc/2)/fs; % Time axis generating 71 72 % Output 73 W = fax; 74 T = tax; 75 S = STFT3; 76 77 78 %% matlab自带函数 79 figure(); 80 spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs); 81 title('spectrogram函数画图') 82 [spec_X,spec_f,spec_t]=spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs); 83 %减法，看看差距 84 figure() 85 plot(abs(spec_X)-abs(STFT_X)) 86 87 %% 画频谱图 88 figure 89 %利用了SIFT3 90 surf(W, T, S') 91 shading interp 92 axis tight 93 box on 94 view(0, 90) 95 set(gca, 'FontName', 'Times New Roman', 'FontSize', 14) 96 xlabel('Frequency, Hz') 97 ylabel('Time, s') 98 title('Amplitude spectrogram of the signal') 99 handl = colorbar; 100 set(handl, 'FontName', 'Times New Roman', 'FontSize', 14) 101 ylabel(handl, 'Magnitude, dB') 102 103 %跳频图 104 figure(); 105 surf(T,W,S) 106 shading interp 107 axis tight 108 box on 109 view(0, 90) 110 set(gca, 'FontName', 'Times New Roman', 'FontSize', 14) 111 ylabel('Frequency, Hz') 112 xlabel('Time, s') 113 title('Amplitude spectrogram of the signal') 114 handl = colorbar; 115 set(handl, 'FontName', 'Times New Roman', 'FontSize', 14) 116 ylabel(handl, 'Magnitude, dB')
主要记到行数为nrow=nfft/2+1行。
重点关注代码8行，57行，62行，85行。
二、分析思路
1）自己在使用fft实现时，通过一个for循环，对每个nsc长度的片段进行加窗、nfft点的fft变换。
2）
1. 当输入信号为实数时，由于fft的对称性，从结果SIFT1仅取fft结果行数的一半放到结果矩阵SIFT_X中去。
2. 当输入信号为复数时，该矩阵的行数仍然为nff行，因此结果集即为SIFT1。
3. 此时每一列即是一个小片段的fft变换的结果，该列结果是复数。 可以和Matlab API得到的结果进行差值比较，发现结果完全一样。
这3处分析的不同，也正是两个自己实现代码的修改处，可参考上述代码。
3）后续由于考虑到绘图、幅值、频谱的影响，因此又在复数矩阵的结果上进行了一些运算，因此在后续封装为函数时根据需要进行改进。