Matlab周期圖法使用FFT實現

本文轉載自查看原文 2018-12-21 17:20 1817 matlab

參考文章：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

View Code

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')

View Code

重點關注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')

View Code

主要記到行數為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）后續由於考慮到繪圖、幅值、頻譜的影響，因此又在復數矩陣的結果上進行了一些運算，因此在后續封裝為函數時根據需要進行改進。

免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。

猜您在找 MatLab實現FFT與功率譜 FFT-Matlab初步實現 Voronoi圖及matlab實現基2時抽8點FFT的matlab實現流程及FFT的內部機理點九切圖法使用 MATLAB中fft函數的正確使用方法 matlab 中fft的用法 Matlab進行FFT變換使用axes函數在matlab繪圖中實現圖中圖的繪制 MATLAB中fft