Matlab周期图法使用FFT实现
生活随笔
收集整理的這篇文章主要介紹了
Matlab周期图法使用FFT实现
小編覺得挺不錯的,現(xiàn)在分享給大家,幫大家做個參考.
參考文章：http://www.cnblogs.com/adgk07/p/9314892.html
首先根據(jù)他這個代碼和我之前手上已經(jīng)擁有的那個代碼，編寫了一個適合自己的代碼。
首先模仿他的代碼，測試成功。
思路：
短時傅里葉變換，其實(shí)還是傅里葉變換，只不過把一段長信號按信號長度（nsc）、重疊點(diǎn)數(shù)（nov）重新采樣。
% 結(jié)合之前兩個版本的stft，實(shí)現(xiàn)自己的周期圖法，力求通俗易懂，代碼分明。% 該代碼寫的時候是按照輸入信號為實(shí)數(shù)的思路寫的，在每個片段fft時進(jìn)行前一半行的轉(zhuǎn)置存儲。后續(xù)代碼思路已修改。 clear; clc; close all % %---------------------Chirp 信號生成(start)--------------------- fa = [ 50 800 ]; % 掃描頻率上下限 fs = 6400; % 采樣頻率 % 分辨率相關(guān)設(shè)定參數(shù) te = 1; % [s] 掃頻時間長度 fre = 8; % [s] 頻率分辨率 tre = 0.002; % [Hz] 時間分辨率 t = 0:1/fs:te; % [s] 時間序列 sc = chirp(t,fa(1),te,fa(2)); % 信號生成%待傳參(實(shí)參) signal = sc; nsc = floor(fs/fre); nff = 2^nextpow2(nsc);%每個窗口進(jìn)行fft的長度 nov = floor(nsc*0.9); % %---------------------Chirp 信號生成(end)--------------------- %% 使用fft實(shí)現(xiàn)周期圖法 % Input % signal - Signal vector 輸入信號(行向量) % nsc - Number SeCtion 每個小分段的信號長度 % nov - Number OverLap 重疊點(diǎn)數(shù) % fs - Frequency of Sample 采樣率 % Output % S - A matrix that each colum is a FFT for time of nsc % 如果nfft為偶數(shù)，則S的行數(shù)為(nfft/2+1)，如果nfft為奇數(shù)，則行數(shù)為(nfft+1)/2 % 每一列是一個小窗口的FFT結(jié)果，因?yàn)閙atlab的FFT結(jié)果是對稱的，只需要一半 % W - A vector labeling frequency% T - A vector labeling time % Signal Preprocessing h = hamming(nsc, 'periodic'); % Hamming weight function 海明窗加權(quán)，窗口大小為nsc L = length(signal); % Length of Signal 信號長度 nstep = nsc-nov; % Number of STep per colum 每個窗口的步進(jìn) 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,初始化最終結(jié)果 index=1;%當(dāng)前片段第一個信號位置在原始信號中的索引 for i=1:ncol%提取當(dāng)前片段信號值,并用海明窗進(jìn)行加權(quán)(均為長度為nsc的行向量)temp_S=signal(index:index+nsc-1).*h';%對長度為nsc的片段進(jìn)行nfft點(diǎn)FFT變換temp_X=fft(temp_S,nfft);%從長度為nfft點(diǎn)(行向量)的fft變換中取一半(轉(zhuǎn)換為列向量)，存儲到矩陣的列向量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) % 繪制自編函數(shù)結(jié)果my_pcolor(W,T,S)caxis([-130.86,-13.667]);title('自編函數(shù)');xlabel('時間 second');ylabel('頻率 Hz');subplot(1,3,2) % 繪制 Matlab 函數(shù)結(jié)果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 自帶函數(shù)');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 實(shí)現(xiàn)與比較');
內(nèi)部調(diào)用了my_pcolor.m
function [ ] = my_pcolor( f , t , s ) %MY_PCOLOR 繪圖函數(shù) % Input % f - 頻率軸矩陣 % t - 時間軸矩陣 % s - 幅度軸矩陣gca = pcolor(f,t,s); % 繪制偽彩色圖set(gca, 'LineStyle','none'); % 取消網(wǎng)格，避免一片黑handl = colorbar; % 彩圖坐標(biāo)尺set(handl, 'FontName', 'Times New Roman', 'FontSize', 8)ylabel(handl, 'Magnitude, dB') % 坐標(biāo)尺名稱 end View Code?
function [ ] = my_pcolor( f , t , s ) %MY_PCOLOR 繪圖函數(shù) % Input % f - 頻率軸矩陣 % t - 時間軸矩陣 % s - 幅度軸矩陣gca = pcolor(f,t,s); % 繪制偽彩色圖set(gca, 'LineStyle','none'); % 取消網(wǎng)格，避免一片黑handl = colorbar; % 彩圖坐標(biāo)尺set(handl, 'FontName', 'Times New Roman', 'FontSize', 8)ylabel(handl, 'Magnitude, dB') % 坐標(biāo)尺名稱 end?
接下來在主體不變的情況下，修改輸入信號和函數(shù)返回結(jié)果，使其和自己之前的代碼效果相同。
?
其實(shí)只要搞清楚了理論，實(shí)現(xiàn)上很簡單。
一、首先分析Matlab周期圖法API。
[spec_X,spec_f,spec_t]=spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs);如果函數(shù)沒有返回結(jié)果的調(diào)用，MATLAB直接幫我們繪圖了。
1）當(dāng)輸入信號為復(fù)數(shù)時的代碼：
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點(diǎn)采樣長度 %% matlab自帶函數(shù) figure(1) spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs); title('spectrogram函數(shù)畫圖') [spec_X,spec_f,spec_t]=spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);變量：
其中：
spec_X矩陣行數(shù)為nff,列數(shù)是根據(jù)信號signal的長度和每個片段nsc分割決定的，共計col列。
spec_f 為nff*1的列向量。
spec_t為 1*ncol的行向量。
此時自己實(shí)現(xiàn)的Matlab代碼為：
1 % 結(jié)合之前兩個版本的stft，實(shí)現(xiàn)自己的周期圖法，力求通俗易懂，代碼分明。 2 clear; clc; close all 3 %% 信號輸入基本參數(shù) 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));%每個窗口進(jìn)行fft的長度 15 nov = 0; %重疊率 16 fs = Fs; %采樣率 17 18 %% 使用fft實(shí)現(xiàn)周期圖法 19 %后續(xù)可封裝為函數(shù): 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 重疊點(diǎn)數(shù) 26 % fs - Frequency of Sample 采樣率 27 % Output 28 % S - A matrix that each colum is a FFT for time of nsc 29 % 如果nfft為偶數(shù)，則S的行數(shù)為(nfft/2+1)，如果nfft為奇數(shù)，則行數(shù)為(nfft+1)/2 30 % 每一列是一個小窗口的FFT結(jié)果，因?yàn)閙atlab的FFT結(jié)果是對稱的，只需要一半 31 % W - A vector labeling frequency 頻率軸 32 % T - A vector labeling time 時間軸 33 % Signal Preprocessing 34 h = hamming(nsc, 'periodic'); % Hamming weight function 海明窗加權(quán)，窗口大小為nsc 35 L = length(signal); % Length of Signal 信號長度 36 nstep = nsc-nov; % Number of STep per colum 每個窗口的步進(jìn) 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 初始化最終結(jié)果 43 index = 1;%當(dāng)前片段第一個信號位置在原始信號中的索引 44 for i=1:ncol 45 %提取當(dāng)前片段信號值,并用海明窗進(jìn)行加權(quán)(均為長度為nsc的行向量) 46 temp_S=signal(index:index+nsc-1).*h'; 47 %對長度為nsc的片段進(jìn)行nfft點(diǎn)FFT變換 48 temp_X=fft(temp_S,nfft); 49 %將長度為nfft點(diǎn)(行向量)的fft變換轉(zhuǎn)換為列向量，存儲到矩陣的列向量 50 STFT1(:,i)=temp_X'; 51 %將索引后移 52 index=index + nstep; 53 end 54 55 56 %如果是為了和標(biāo)準(zhǔn)周期圖法進(jìn)行誤差比較，則后續(xù)計算(abs,*2,log(+1e-6))不需要做，只需 57 %plot(abs(spec_X)-abs(STFT1))比較即可 58 59 STFT2 = abs(STFT1/nfft); 60 %nfft是偶數(shù)，說明首尾是奈奎斯特點(diǎn)，只需對2:end-1的數(shù)據(jù)乘以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自帶函數(shù) 76 figure(); 77 spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs); 78 title('spectrogram函數(shù)畫圖') 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重點(diǎn)關(guān)注8行、59行、66行、69行、82行
其實(shí)此時只要牢牢記住結(jié)果集中行數(shù)為nfft就行了。
2）當(dāng)輸入信號為實(shí)數(shù)時的代碼：
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點(diǎn)采樣長度
%% matlab自帶函數(shù)
figure(1)
spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);
title('spectrogram函數(shù)畫圖')
[spec_X,spec_f,spec_t]=spectrogram(S,hamming(nsc, 'periodic'),nov,nff,Fs);
?
其中：
spec_X矩陣行數(shù)為nff/2+1,列數(shù)是根據(jù)信號signal的長度和每個片段nsc分割決定的，共計col列。
spec_f 為nff/2+1*1的列向量。
spec_t為 1*ncol的行向量。
主要是因?yàn)檩斎胄盘枮閷?shí)數(shù)時，fft變換的結(jié)果有一半是對稱的，不必考慮。
自己實(shí)現(xiàn)的Matlab代碼為：
1 % 結(jié)合之前兩個版本的stft，實(shí)現(xiàn)自己的周期圖法，力求通俗易懂，代碼分明。 2 clear; clc; close all 3 %% 信號輸入基本參數(shù) 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));%每個窗口進(jìn)行fft的長度 15 nov = 0; %重疊率 16 fs = Fs; %采樣率 17 18 %% 使用fft實(shí)現(xiàn)周期圖法 19 %后續(xù)可封裝為函數(shù): 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 重疊點(diǎn)數(shù) 26 % fs - Frequency of Sample 采樣率 27 % Output 28 % S - A matrix that each colum is a FFT for time of nsc 29 % 如果nfft為偶數(shù)，則S的行數(shù)為(nfft/2+1)，如果nfft為奇數(shù)，則行數(shù)為(nfft+1)/2 30 % 每一列是一個小窗口的FFT結(jié)果，因?yàn)閙atlab的FFT結(jié)果是對稱的，只需要一半 31 % W - A vector labeling frequency 頻率軸 32 % T - A vector labeling time 時間軸 33 % Signal Preprocessing 34 h = hamming(nsc, 'periodic'); % Hamming weight function 海明窗加權(quán)，窗口大小為nsc 35 L = length(signal); % Length of Signal 信號長度 36 nstep = nsc-nov; % Number of STep per colum 每個窗口的步進(jìn) 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 初始化最終結(jié)果 43 index = 1;%當(dāng)前片段第一個信號位置在原始信號中的索引 44 for i=1:ncol 45 %提取當(dāng)前片段信號值,并用海明窗進(jìn)行加權(quán)(均為長度為nsc的行向量) 46 temp_S=signal(index:index+nsc-1).*h'; 47 %對長度為nsc的片段進(jìn)行nfft點(diǎn)FFT變換 48 temp_X=fft(temp_S,nfft); 49 %將長度為nfft點(diǎn)(行向量)的fft變換轉(zhuǎn)換為列向量，存儲到矩陣的列向量 50 STFT1(:,i)=temp_X'; 51 %將索引后移 52 index=index + nstep; 53 end 54 55 % -----當(dāng)輸入信號為實(shí)數(shù)時，對其的操作（也就是說只考慮信號的前一半） 56 % 由于實(shí)數(shù)FFT的對稱性,提取STFT1的nrow行 57 STFT_X = STFT1(1:nrow,:); % Symmetric results of FFT 58 59 %如果是為了和標(biāo)準(zhǔn)周期圖法進(jìn)行誤差比較，則后續(xù)計算(abs,*2,log(+1e-6))不需要做，只需 60 %plot(abs(spec_X)-abs(STFT_X))比較即可 61 62 STFT2 = abs(STFT_X/nfft); 63 %nfft是偶數(shù)，說明首尾是奈奎斯特點(diǎn)，只需對2:end-1的數(shù)據(jù)乘以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自帶函數(shù) 79 figure(); 80 spectrogram(signal,hamming(nsc, 'periodic'),nov,nff,Fs); 81 title('spectrogram函數(shù)畫圖') 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主要記到行數(shù)為nrow=nfft/2+1行。
重點(diǎn)關(guān)注代碼8行，57行，62行，85行。
二、分析思路
1）自己在使用fft實(shí)現(xiàn)時，通過一個for循環(huán)，對每個nsc長度的片段進(jìn)行加窗、nfft點(diǎn)的fft變換。
2）
1. 當(dāng)輸入信號為實(shí)數(shù)時，由于fft的對稱性，從結(jié)果SIFT1僅取fft結(jié)果行數(shù)的一半放到結(jié)果矩陣SIFT_X中去。
2. 當(dāng)輸入信號為復(fù)數(shù)時，該矩陣的行數(shù)仍然為nff行，因此結(jié)果集即為SIFT1。
3. 此時每一列即是一個小片段的fft變換的結(jié)果，該列結(jié)果是復(fù)數(shù)。 可以和Matlab API得到的結(jié)果進(jìn)行差值比較，發(fā)現(xiàn)結(jié)果完全一樣。
這3處分析的不同，也正是兩個自己實(shí)現(xiàn)代碼的修改處，可參考上述代碼。
3）后續(xù)由于考慮到繪圖、幅值、頻譜的影響，因此又在復(fù)數(shù)矩陣的結(jié)果上進(jìn)行了一些運(yùn)算，因此在后續(xù)封裝為函數(shù)時根據(jù)需要進(jìn)行改進(jìn)。
?
轉(zhuǎn)載于:https://www.cnblogs.com/shuqingstudy/p/10155980.html
總結(jié)
以上是生活随笔為你收集整理的Matlab周期图法使用FFT实现的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： banner手动切换效果
- 下一篇： windows下多进程加协程并发模式