Matlab周期图法使用FFT实现

参考文章：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代码为：

% 结合之前两个版本的stft，实现自己的周期图法，力求通俗易懂，代码分明。
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;


%传参
signal = S;
nsc    = 100;                     %每个窗口的长度,也即海明窗的长度
nff    = max(256,2^nextpow2(nsc));%每个窗口进行fft的长度
nov    = 0;                       %重叠率
fs     = Fs;                      %采样率

%% 使用fft实现周期图法
%后续可封装为函数:
% function [ S , W , T ] = mf_spectrogram...
%     ( signal , nsc , nov , fs )
% 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   = max(256,2^nextpow2(nsc));   % Number of Fast Fourier Transformation  每个窗口FFT长度
nrow   = nfft/2+1;
%
%
STFT1 = zeros(nfft,ncol);  %STFT1 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变换转换为列向量，存储到矩阵的列向量
    STFT1(:,i)=temp_X';
    %将索引后移
    index=index + nstep;
end


%如果是为了和标准周期图法进行误差比较，则后续计算(abs,*2,log(+1e-6))不需要做，只需
%plot(abs(spec_X)-abs(STFT1))比较即可

STFT2 = abs(STFT1/nfft);
%nfft是偶数，说明首尾是奈奎斯特点，只需对2:end-1的数据乘以2
STFT2(2:end-1,:) = 2*STFT2(2:end-1,:);  % Add the value of the other half
%STFT3 = 20*log10(STFT1);                % Turn sound pressure into dB level
STFT3 = 20*log10(STFT2 + 1e-6);          % convert amplitude spectrum to dB (min = -120 dB)

% Axis Generating
fax =(0:(nfft-1)) * fs./nfft;                        % Frequency axis setting
tax = (nsc/2 : nstep : nstep*(ncol-1)+nsc/2)/fs;     % Time axis generating
 
 % Output
 W = fax;
 T = tax;
 S = STFT3;
 

%% matlab自带函数
figure();
spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);
title('spectrogram函数画图')
[spec_X,spec_f,spec_t]=spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);
%减法，看看差距
figure()
plot(abs(spec_X)-abs(STFT1))

%% 画频谱图
figure
%利用了SIFT3
surf(W, T,  S')
shading interp
axis tight
box on
view(0, 90)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
xlabel('Frequency, Hz')
ylabel('Time, s')
title('Amplitude spectrogram of the signal')
handl = colorbar;
set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(handl, 'Magnitude, dB')

%跳频图
figure();
surf(T,W,S)
shading interp
axis tight
box on
view(0, 90)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel('Frequency, Hz')
xlabel('Time, s')
title('Amplitude spectrogram of the signal')
handl = colorbar;
set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
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代码为：

% 结合之前两个版本的stft，实现自己的周期图法，力求通俗易懂，代码分明。
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);


%传参
signal = S;
nsc    = 100;                     %每个窗口的长度,也即海明窗的长度
nff    = max(256,2^nextpow2(nsc));%每个窗口进行fft的长度
nov    = 0;                       %重叠率
fs     = Fs;                      %采样率

%% 使用fft实现周期图法
%后续可封装为函数:
% function [ S , W , T ] = mf_spectrogram...
%     ( signal , nsc , nov , fs )
% 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   = max(256,2^nextpow2(nsc));   % Number of Fast Fourier Transformation  每个窗口FFT长度
nrow   = nfft/2+1;
%
%
STFT1 = zeros(nfft,ncol);  %STFT1 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变换转换为列向量，存储到矩阵的列向量
    STFT1(:,i)=temp_X';
    %将索引后移
    index=index + nstep;
end

% -----当输入信号为实数时，对其的操作（也就是说只考虑信号的前一半）
% 由于实数FFT的对称性,提取STFT1的nrow行
STFT_X = STFT1(1:nrow,:);             % Symmetric results of FFT

%如果是为了和标准周期图法进行误差比较，则后续计算(abs,*2,log(+1e-6))不需要做，只需
%plot(abs(spec_X)-abs(STFT_X))比较即可

STFT2 = abs(STFT_X/nfft);
%nfft是偶数，说明首尾是奈奎斯特点，只需对2:end-1的数据乘以2
STFT2(2:end-1,:) = 2*STFT2(2:end-1,:);  % Add the value of the other half
%STFT3 = 20*log10(STFT1);                % Turn sound pressure into dB level
STFT3 = 20*log10(STFT2 + 1e-6);          % convert amplitude spectrum to dB (min = -120 dB)

% Axis Generating
fax =(0:(nrow-1)) * fs./nfft;                        % Frequency axis setting
tax = (nsc/2 : nstep : nstep*(ncol-1)+nsc/2)/fs;     % Time axis generating
 
 % Output
 W = fax;
 T = tax;
 S = STFT3;
 

%% matlab自带函数
figure();
spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);
title('spectrogram函数画图')
[spec_X,spec_f,spec_t]=spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);
%减法，看看差距
figure()
plot(abs(spec_X)-abs(STFT_X))

%% 画频谱图
figure
%利用了SIFT3
surf(W, T,  S')
shading interp
axis tight
box on
view(0, 90)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
xlabel('Frequency, Hz')
ylabel('Time, s')
title('Amplitude spectrogram of the signal')
handl = colorbar;
set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(handl, 'Magnitude, dB')

%跳频图
figure();
surf(T,W,S)
shading interp
axis tight
box on
view(0, 90)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel('Frequency, Hz')
xlabel('Time, s')
title('Amplitude spectrogram of the signal')
handl = colorbar;
set(handl, 'FontName', 'Times New Roman', 'FontSize', 14)
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）后续由于考虑到绘图、幅值、频谱的影响，因此又在复数矩阵的结果上进行了一些运算，因此在后续封装为函数时根据需要进行改进。