飞码网-免费源码博客分享网站

点击这里给我发消息

matlab如何对音频文件进行短时傅里叶变换(STFT)(自定义函数)

飞码网-免费源码博客分享网站 爱上飞码网—https://www.codefrees.com— 飞码网-matlab-python-C++ 爱上飞码网—https://www.codefrees.com— 飞码网-免费源码博客分享网站
本代码是一个Matlab函数,提供给定信号x[n]的短时傅里叶变换(STFT)。该函数是Matlab命令 "spectrogram "的替代方案。该函数的输出是
1) 一个矩阵,包含复数STFT系数,列间为时间,行间为频率。
2)一个频率向量。
3)一个时间向量。

为了明确函数的用法,现举例说明。为方便起见,在函数的开头给出了输入和输出参数。

案例是将一个wave文件进行STFT短时傅里叶变换

其中:输入和输出如下。如果需要请更换为自己需要的输入


% Input:
% x - signal in the time domain
% win - analysis window function
% hop - hop size
% nfft - number of FFT points
% fs - sampling frequency, Hz
%
% Output:
% STFT - STFT-matrix (only unique points, time 
%        across columns, frequency across rows)
% f - frequency vector, Hz
% t - time vector, s

本程序读入的wave文件在文章末尾给出,测试的时候将其和运行.m文件放到同一个文件夹即可

clear, clc, close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              Short-Time Fourier Transform            %
%               with MATLAB Implementation             %
%                                                      %
%        飞码网—codefrees.com-版权所有        20/10/13 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load an audio file
[x, fs] = audioread('track.wav');   % load an audio file
x = x(:, 1);                        % get the first channel

% define analysis parameters
wlen = 1024;                        % window length (recomended to be power of 2)
hop = wlen/4;                       % hop size (recomended to be power of 2)
nfft = 4096;                        % number of fft points (recomended to be power of 2)

% perform STFT
win = blackman(wlen, 'periodic');
[S, f, t] = stft(x, win, hop, nfft, fs);

% calculate the coherent amplification of the window
C = sum(win)/wlen;

% take the amplitude of fft(x) and scale it, so not to be a
% function of the length of the window and its coherent amplification
S = abs(S)/wlen/C;

% correction of the DC & Nyquist component
if rem(nfft, 2)                     % odd nfft excludes Nyquist point
    S(2:end, :) = S(2:end, :).*2;
else                                % even nfft includes Nyquist point
    S(2:end-1, :) = S(2:end-1, :).*2;
end

% convert amplitude spectrum to dB (min = -120 dB)
S = 20*log10(S + 1e-6);

% plot the spectrogram
figure(1)
surf(t, f, S)
shading interp
axis tight
view(0, 90)
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
xlabel('Time, s')
ylabel('Frequency, Hz')
title('Amplitude spectrogram of the signal')

hcol = colorbar;
set(hcol, 'FontName', 'Times New Roman', 'FontSize', 14)
ylabel(hcol, 'Magnitude, dB')



function [STFT, f, t] = stft(x, win, hop, nfft, fs)

% function: [STFT, f, t] = stft(x, win, hop, nfft, fs)
%
% Input:
% x - signal in the time domain
% win - analysis window function
% hop - hop size
% nfft - number of FFT points
% fs - sampling frequency, Hz
%
% Output:
% STFT - STFT-matrix (only unique points, time 
%        across columns, frequency across rows)
% f - frequency vector, Hz
% t - time vector, s

% representation of the signal as column-vector
x = x(:);

% determination of the signal length 
xlen = length(x);

% determination of the window length
wlen = length(win);

% stft matrix size estimation and preallocation
NUP = ceil((1+nfft)/2);     % calculate the number of unique fft points
L = 1+fix((xlen-wlen)/hop); % calculate the number of signal frames
STFT = zeros(NUP, L);       % preallocate the stft matrix

% STFT (via time-localized FFT)
for l = 0:L-1
    % windowing
    xw = x(1+l*hop : wlen+l*hop).*win;
    
    % FFT
    X = fft(xw, nfft);
    
    % update of the stft matrix
    STFT(:, 1+l) = X(1:NUP);
end

% calculation of the time and frequency vectors
t = (wlen/2:hop:wlen/2+(L-1)*hop)/fs;
f = (0:NUP-1)*fs/nfft;

end


 
track.zip
c552608aa9ba50eed841a8391c1061bd.zip (333.37 KB)
飞码网-免费源码博客分享网站 爱上飞码网—https://www.codefrees.com— 飞码网-matlab-python-C++ 爱上飞码网—https://www.codefrees.com— 飞码网-免费源码博客分享网站
赞 ()
内容页底部广告位3
留言与评论(共有 0 条评论)
   
验证码: