Manual implementation of spectogram function

2021 4 29
2 个回答
回答已采纳
更新时间：2021 4 29
6 次查看（30 天）

0 个投票

Hello,
I am trying to implement my own function that gives the same results as Matlab spectogram function.
So far I have accomplished function like this:
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
In my tests it works perfectly and gives the same results like spectogram function when parameter nfft is greater or equal to length of window.
% 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
In the case when length of window is less than nfft than the results are different.
% 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
So how can I improve my function to recieve the same results even in the case when nnft is less than window length? How Matlab's spectogram calculates this case?

请先登录，再回答此问题。

 采纳的回答

Pawel Krzyzak
Pawel Krzyzak 2021-4-29
编辑：Pawel Krzyzak 2021-4-29

0 个投票

I noticed that when window size is greater than nfft scalar number, the data has to be transformed somehow. Finally I found an inner Matlab function that probably is called in the original spectogram Matlab function. It is named datawrap and wraps input data modulo nfft.
So in my function I had to transform data segment (in the same way how datawrap function does it) before calling fft. Improved function:
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
I believe it works properly because it gives the same results as Matlab spectogram function.

更多回答（1 个）

Mathieu NOE
Mathieu NOE 2021-4-29

0 个投票

hello
I share with you my own version , hope it helps !
notice i always use hanning window and the correction factor for amplitude is related to this choice
fs=5000; % Sampling rate
N=1000; % Number of data points
t=[0:1:N-1]/fs;
y=sin(500*pi*t)...
+2*sin(500*1.5*pi*t)...
+3*sin(500*3*pi*t)...
+4*sin(500*4*pi*t);
figure(1);
windowsize = 200;
window = boxcar(windowsize);
nfft = windowsize;
noverlap = windowsize-1;
% [S,F,T] = spectrogram(y,window,noverlap,nfft,fs); % S = SPECTROGRAM(X,WINDOW,NOVERLAP,NFFT,Fs)
[S,F,T] = myspecgram(y, fs, nfft, 0.75); % overlap is 75% of nfft here
imagesc(T,F,(abs(S)))
colorbar('vert');
set(gca,'YDir','Normal')
xlabel('Time (secs)')
ylabel('Freq (Hz)')
title('Short-time Fourier Transform spectrum')
colormap('jet');
function [fft_specgram,freq_vector,time] = myspecgram(signal, Fs, nfft, Overlap)
% FFT peak spectrogram of signal (example sinus amplitude 1 = 0 dB after fft).
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
s_tmp = zeros(nfft,1);
s_tmp((1:samples),:) = signal;
signal = s_tmp;
samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
% compute fft with overlap
offset = fix((1-Overlap)*nfft);
spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
% % for info is equivalent to :
% noverlap = Overlap*nfft;
% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows
% main loop
fft_specgram = [];
for ci=1:spectnum
start = (ci-1)*offset;
sw = signal((1+start):(start+nfft)).*window;
fft_specgram = [fft_specgram abs(fft(sw))*4/nfft]; % X=fft(x.*hanning(N))*4/N; % hanning only
end
% one sidded fft spectrum % Select first half
if rem(nfft,2) % nfft odd
select = (1:(nfft+1)/2)';
else
select = (1:nfft/2+1)';
end
fft_specgram = fft_specgram(select,:);
freq_vector = (select - 1)*Fs/nfft;
% time vector
% time stamps are defined in the middle of the buffer
time = ((0:spectnum-1)*offset + round(offset/2))/Fs;
end

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Pawel Krzyzak
Pawel Krzyzak 2021-4-29
Thank your very much for your implementation but I still don't know how can I modify my function to recieve the same results as results given by matlab spectogram funcion.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Time-Frequency Analysis 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by