手工实现Matlab spectogram函数

Question

我正在尝试实现我自己的函数，它给出了与Matlab spectogram函数相同的结果。到目前为止，我已经完成了如下功能：

function out = manulaSpectogram(x, win, noverlap, nfft)
x = x(:);
n = length(x);
wlen = length(win);
nUnique = ceil((1+nfft)/2); % number of uniqure points
L = fix((n-noverlap)/(wlen-noverlap)); % number of signal frames
out = zeros(L, nUnique);
index = 1:wlen;
    for i = 0:L-1
        xw = win.*x(index);
        X = fft(xw, nfft);
        out(i+1, :) = X(1:nUnique);
        index = index + (wlen - noverlap);
    end
end

在我的测试中，它工作得很好，当参数nfft大于或等于窗口长度时，它会给出与spectogram函数相同的结果。

% first test (nnft = window length):
A = [1,2,3,4,5,6];
window = 6;
overlap = 2;
nfft = 6;
s = spectrogram(A, hamming(window), overlap, nfft)'
s2 = manulaSpectogram(A, hamming(window), overlap, nfft)
% results:
s =
   9.7300 + 0.0000i  -5.2936 + 0.9205i   0.7279 - 0.3737i  -0.1186 + 0.0000i
s2 =
   9.7300 + 0.0000i  -5.2936 - 0.9205i   0.7279 + 0.3737i  -0.1186 + 0.0000i
   
% second test (nfft > window length):
A = [1,2,3,4,5,6];
window = 3;
overlap = 2;
nfft = 6;
s = spectrogram(A, hamming(window), overlap, nfft)'
s2 = manulaSpectogram(A, hamming(window), overlap, nfft)
% results:
s =
   2.3200 + 0.0000i   0.9600 + 1.9399i  -1.0400 + 1.5242i  -1.6800 + 0.0000i
   3.4800 + 0.0000i   1.5000 + 2.8752i  -1.5000 + 2.3209i  -2.5200 + 0.0000i
   4.6400 + 0.0000i   2.0400 + 3.8105i  -1.9600 + 3.1177i  -3.3600 + 0.0000i
   5.8000 + 0.0000i   2.5800 + 4.7458i  -2.4200 + 3.9144i  -4.2000 + 0.0000i
   
s2 =
   2.3200 + 0.0000i   0.9600 - 1.9399i  -1.0400 - 1.5242i  -1.6800 + 0.0000i
   3.4800 + 0.0000i   1.5000 - 2.8752i  -1.5000 - 2.3209i  -2.5200 + 0.0000i
   4.6400 + 0.0000i   2.0400 - 3.8105i  -1.9600 - 3.1177i  -3.3600 + 0.0000i
   5.8000 + 0.0000i   2.5800 - 4.7458i  -2.4200 - 3.9144i  -4.2000 + 0.0000i

在窗口长度小于nfft的情况下，结果不同。

% third test (nfft < window length):
A = [1,2,3,4,5,6];
window = 6;
overlap = 2;
nfft = 3;
s = spectrogram(A, hamming(window), overlap, nfft)'
s2 = manulaSpectogram(A, hamming(window), overlap, nfft)
% results:
s =
   9.7300 + 0.0000i   0.7279 - 0.3737i
 
s2 =
   3.6121 + 0.0000i  -1.6861 + 1.6807i

那么，我如何改进我的函数，即使在nnft小于窗口长度的情况下也能收到相同的结果？Matlab的spectogram是如何计算这种情况的？

我正在尝试实现我自己的函数，因为spectogram函数是一个大型算法的一部分，我需要从Matlab到C#语言来实现它，所以我想知道spectogram“黑盒”是做什么的。

Answer 1

我注意到，当窗口大小大于nfft标量数时，数据必须以某种方式进行转换。最后，我找到了一个内部Matlab函数，它可能是在原始的spectogram Matlab函数中调用的。它被命名为datawrap，包装了模数为nfft的输入数据。

所以在我的函数中，我必须在调用fft之前转换数据段(就像datawrap函数一样)。改进后的函数：

function out = manulaSpectogram(x, win, noverlap, nfft)
x = x(:);
n = length(x);
wlen = length(win);
nUnique = ceil((1+nfft)/2); % number of uniqure points
L = fix((n-noverlap)/(wlen-noverlap)); % number of signal frames
out = zeros(L, nUnique);
index = 1:wlen;
    for i = 0:L-1
        xw = win.*x(index);
        % added transformation
        if length(xw) > nfft
            xw = sum(buffer(xw, nfft), 2);
        end
        % end of added transformation
        X = fft(xw, nfft);
        out(i+1, :) = X(1:nUnique);
        index = index + (wlen - noverlap);
    end
end

我相信它可以正常工作，因为它给出了与Matlab spectogram函数相同的结果。