Matlab周期图法使用FFT实现
生活随笔
收集整理的這篇文章主要介紹了
Matlab周期图法使用FFT实现
小編覺得挺不錯的,現在分享給大家,幫大家做個參考.
參考文章：http://www.cnblogs.com/adgk07/p/9314892.html
首先根據他這個代碼和我之前手上已經擁有的那個代碼，編寫了一個適合自己的代碼。
首先模仿他的代碼，測試成功。
思路：
短時傅里葉變換，其實還是傅里葉變換，只不過把一段長信號按信號長度（nsc）、重疊點數（nov）重新采樣。
% 結合之前兩個版本的stft，實現自己的周期圖法，力求通俗易懂，代碼分明。% 該代碼寫的時候是按照輸入信號為實數的思路寫的，在每個片段fft時進行前一半行的轉置存儲。后續(xù)代碼思路已修改。 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 dBSTFT1 = abs(STFT_X/nfft);STFT1(2:end-1,:) = 2*STFT1(2:end-1,:); % Add the value of the other halfSTFT3 = 20*log10(STFT1); % Turn sound pressure into dB level% Axis Generatingfax =(0:(nfft/2)) * fs/nfft; % Frequency axis settingtax = (nsc/2 : nstep : nstep*(ncol-1)+nsc/2)/fs; % Time axis generating[ffax,ttax] = meshgrid(tax,fax); % Generating grid of figure% OutputW = 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 dBs1 = abs(s/nff);s2 = s1(1:nff/2+1,:); % Symmetric results of FFTs2(2:end-1,:) = 2*s2(2:end-1,:); % Add the value of the other halfs3 = 20*log10(s2); % Turn sound pressure into dB levelmy_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 %后續(xù)可封裝為函數: 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 %如果是為了和標準周期圖法進行誤差比較，則后續(xù)計算(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 %后續(xù)可封裝為函數: 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 %如果是為了和標準周期圖法進行誤差比較，則后續(xù)計算(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循環(huán)，對每個nsc長度的片段進行加窗、nfft點的fft變換。
2）
1. 當輸入信號為實數時，由于fft的對稱性，從結果SIFT1僅取fft結果行數的一半放到結果矩陣SIFT_X中去。
2. 當輸入信號為復數時，該矩陣的行數仍然為nff行，因此結果集即為SIFT1。
3. 此時每一列即是一個小片段的fft變換的結果，該列結果是復數。 可以和Matlab API得到的結果進行差值比較，發(fā)現結果完全一樣。
這3處分析的不同，也正是兩個自己實現代碼的修改處，可參考上述代碼。
3）后續(xù)由于考慮到繪圖、幅值、頻譜的影響，因此又在復數矩陣的結果上進行了一些運算，因此在后續(xù)封裝為函數時根據需要進行改進。
?
轉載于:https://www.cnblogs.com/shuqingstudy/p/10155980.html
總結
以上是生活随笔為你收集整理的Matlab周期图法使用FFT实现的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： banner手动切换效果
- 下一篇： windows下多进程加协程并发模式