傅里葉變換

參照網址短時傅里葉變換

參照上述網址，嘗試編寫短時傅里葉變換得程序，爲了加深對短時傅里葉計算過程的理解，也爲後續的其他算法先打下基礎。編程能力較差，雖然仿照，但可能依然由很多考慮不到。

function [s,f,t] = sftf(x,window,overlap,nfft,fs)
    % x： 輸入的一維數組
    % window： 窗的類型，當輸入爲一整數時，默認使用hanning窗
    % overlap：重疊的長度
    % nfft：做一次傅里葉變換得長度
    % fs：採樣頻率

%     默認參數設置
    if nargin < 2 || isempty(window);   window = hanning(256); end
    if length(window) == 1;             window = hanning(window); end
    if nargin < 3 || isempty(overlap);  overlap = floor(length(window)/2); end
    if nargin < 4 || isempty(nfft);     nfft = length(window); end
    if nargin < 5 || isempty(fs);       fs = 2; end
    if length(window) > nfft;           nfft = length(window); end
    
    n = length(window);
    step = n - overlap;                                     %每次移動的步長
    step_nums = ceil(abs(length(x)-overlap)/abs(step));    %需要移動多少步
    %     需要將x補全,最後不足以領補全
        [x_row, x_col] = size(x);
        if x_row ==1; x = x'; end                               % x變換爲一列
        if x_col ~= 1 && x_row ~=  1; error('輸入的x不是一維數組'); end
        x_len = step * (step_nums-1) + n;
        mm = x_len-length(x);
        x_bu = zeros(mm, 1);
        x = [x; x_bu];
    %     將x進行拆分
        S = zeros(n, step_nums);
        t = 1:step:x_len-n+1;
        for i = 1:length(t)
           S(:,i) = x(t(i):t(i)+n-1).*window;
        end
    
%     進行快速傅里葉變換

    SFTF = fft(S,nfft);
    if rem(n,2)==1; x_fre = (n+1)/2; end
    if rem(n,2)==0; x_fre = n/2+1; end
    SFTF = abs(SFTF/n);
    s = SFTF(1:x_fre,:);
    s(2:end-1,:) = 2*s(2:end-1,:);
    f = (0:x_fre-1)*fs/n;
    t = t/fs;
    
% 畫圖部分    
%     figure
%     surf(t,f,10*log10(abs(s)),'EdgeColor','none');
%     axis xy;
%     axis tight;
%     colormap(jet);
%     view(0,90);
%     xlabel('Time');

    
end

使用數據進行測試結果如下

clc
Fs = 1000;
N = 2000;
t = 0:1/Fs:(N-1)/Fs;
y = 1.5*sin(2*pi*80*t)+1.5*sin(2*pi*160*t)+1.5*sin(2*pi*240*t)+2*sin(2*pi*320*t)+2*sin(2*pi*120*t);
[s,f,t] = sftf(y, 300,100,300,Fs);
figure
surf(t,f,10*log10(abs(s)),'EdgeColor','none');
axis xy;
axis tight;
colormap(jet);
view(0,90);
xlabel('Time');
ylabel('Frequency (Hz)');