请写出matlab的mfcc函数

function [MFCCs, FBEs, frames] = mfcc(samples, fs, Tw, Ts, alpha, window, R, M, N, L)

% Ensure that we have a row vector
if size(samples, 1) > 1
samples = samples';
end

% Explode samples into frames, apply hamming window and FFT each frame
frames = enframe(samples, hamming(round(Tw*fs)), round(Ts*fs));
complex_spectrum = fft(frames, N, 2);
power_spectrum = (1/N) * abs(complex_spectrum).^2;

% Prepare filterbank
[f, ~] = get_filter_banks(fs, N, R, M, L);

% Apply filterbank and log-compress
FBEs = log10(f * power_spectrum(:, 1:(N/2+1))');

% Compute DCT
MFCCs = dct(alpha * FBEs')';

% Keep only first L coefficients
MFCCs = MFCCs(:, 1:L);

end

% Helper function to get filter banks
function [f, fn] = get_filter_banks(fs, N, R, M, L)

% Convert frequencies to mel scale
fmin = 0;
fmax = fs / 2;
melmin = 2595 * log10(1 + (fmin / 700));
melmax = 2595 * log10(1 + (fmax / 700));
melcenters = linspace(melmin, melmax, M+2);

% Convert mel scale back to frequency
fcenters = 700 * (10.^(melcenters / 2595) - 1);

% Round frequencies to nearest FFT bin
bin = floor((N + 1) * fcenters / fs);

% Prepare filterbank matrix
f = zeros(M, L);
for m = 1:M
f(m, :) = triangular_filter(bin(m), bin(m+1), bin(m+2), L);
end

% Normalize filterbanks
fn = zeros(M, L);
for i = 1:M
for j = 1:L
if (i == 1)
sum_low = 0;
else
sum_low = sum(f(1:i-1, j));
end
sum_high = sum(f(i, j));
fn(i, j) = sum_high / (sum_low + sum_high);
end
end

end

% Helper function to create triangular filter
function f = triangular_filter(x_low, x_peak, x_high, L)
k = 1:L;
f = zeros(1, L);
for i = 1:L
if (k(i) < x_low || k(i) > x_high)
f(i) = 0;
elseif (k(i) >= x_low && k(i) < x_peak)
f(i) = 2 * (k(i) - x_low) / ((x_high - x_low) * (x_peak - x_low));
elseif (k(i) >= x_peak && k(i) <= x_high)
f(i) = 2 * (x_high - k(i)) / ((x_high - x_low) * (x_high - x_peak));
end
end
end

% Helper function to explode samples into frames
function frames = enframe(samples, window, shift)
nframes = floor((length(samples)-length(window))/shift) + 1;
frames = zeros(nframes, length(window));
for i = 1:nframes
frames(i, :) = samples((i-1)*shift+1:(i-1)*shift+length(window)) .* window;
end
end