Manual Implementation of STFT of an audio signal.

34 次查看（过去 30 天）
UJJWAL
UJJWAL 2011-9-22
Hi ,
I have an Audio File whose spectrogram I want to find manually by calculating its STFT without using the inbuilt Spectrogram function. My steps can be mentioned as follows :-
a) I take samples out of the file in chunks of 1000 samples each and process them one by one. b) For each chunk , I find its STFT . c) I plot the magnitute squared of the STFT and hold that plot. d) I repeat the whole process again and again till all the samples of the audio file have been used.
However I face the following big problem :-
The process is never complete as everytime the computer runs out of memory. I try to use as few variables as possible and i clear the variables as soon as they are used. I also use the pack command before executing the program. However each time the computer almost hangs after some iterations. I think this is because the hold on command causes MATLAB to keep storing its previous plot data.
Is there any way around this issue ? Please enlighten. I am using an audio file of length 3.5 seconds at a sampling frequency of 44100 Hz.

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

采纳的回答

Wayne King
Wayne King 2011-9-26
Here's what you have to do, something like this:
dt = 1/44100;
t=0:dt:(117204*dt)-dt;
x = chirp(t,1500,1,8000);
S = zeros(501,117);
windowlength = 1e3;
k = 1;
for jj = 1:117
y = x(k:k+windowlength-1)';
ydft = fft(y).*gausswin(1e3);
S(:,jj) = ydft(1:501);
k = k+windowlength;
end
F = 0:44100/1000:44100/2;
T = 0:(1e3*dt):(117000*dt)-(1e3*dt);
surf(T,F,20*log10(abs(S)),'EdgeColor','none')
axis xy; axis tight; colormap(jet); view(0,90);
xlabel('Time');
ylabel('Frequency (Hz)');
Like I said, the segments do not overlap and I haven't been too careful with somethings, like whether I skip a sample between adjacent segments, but this is the general idea.

更多回答（5 个）

Wayne King
Wayne King 2011-9-22
Why are you plotting each segment and holding it? I think you should store the DFT of the individual segments as column vectors (or row vectors) in a matrix and then make a surf plot when you are done.
You will also want to overlap your segments and not block them.
  6 个评论
显示 4更早的评论
UJJWAL
UJJWAL 2011-9-26
Hey,
here is the problem :-
a) I have an audio file with sampling frequency 44100 Hz and with sample length of 117024. I need to find its spectrogram using STFT with Gaussian Window (which is known as the Gabor Transform).
b) To solve the problem, I had two ways to do :-
i) To take the whole file and and then try to process. It cannot be done because in that case the matrix i will get will be 65536 by 117024 where 65536 appears because 2^nextpow2(44100) is 65536.
ii) I tried to take the file in chunks that is 1000 samples at a time and process. In that case the processing can take place. However I have to store the Spectrogram in a matrix so that finally when I have processed all the samples ,I can display all of them together. However again the same memory problem turns up. So I used hold on thinking that it will keep plotting the spectrogram for each chunk of the audio file and finally I will be able to get the full spectrogram. However it turns out that when I use hold on then MATLAB again starts storing the whole data so memory is the problem.
I hope now this problem is clearer. Pls Help
Daniel Shub
Daniel Shub 2011-9-26
Have you looked at my answer? It seems to do what you want minus the windowing. Adding the window in is relatively easy (you need to multiple each column of my z with your window vector) or replace ones(1000,1) in the spectrogram method with a window vector.

请先登录，再进行评论。

Daniel Shub
Daniel Shub 2011-9-23
I am confused ...
Is this roughly what you are trying to do?
x = randn(117424,1); % Make some fake data of the correct size
y = [x; zeros((1000-mod(length(x), 1000)), 1)];
z = reshape(y, 1000, length(y)/1000);
Z = fft(z, 2^16);
mesh(abs(Z))
Then again, I think that is roughly
spectrogram(x, ones(1000, 1), 0, 2^16, 44.1e3)
Wayne King
Wayne King 2011-9-26
Hi, I think you are misunderstanding some things about the STFT. Why do you think the sampling frequency has anything to do with the length of the short-time Fourier transforms? The length of the short-time Fourier transforms depends on the length of the segments you choose to extract from your data.
Also, the number of columns is NOT equal to the number of points in your data vector.
Below I use a Gaussian window of length 1000 and overlap the segments by 900 points. Look at the size of S, the spectrogram matrix.
win = gausswin(1e3);
dt = 1/44100;
t=0:dt:(117204*dt)-dt;
x = chirp(t,1500,1,8000);
[S,F,T,P] = spectrogram(x,win,900,length(win),44100);
surf(T,F,10*log10(P),'EdgeColor','none')
axis xy; axis tight; colormap(jet); view(0,90);
xlabel('Time');
ylabel('Frequency (Hz)');
UJJWAL
UJJWAL 2011-9-26
Actually I do not want to use the spectrogram function. I want to write my own code for finding the spectrogram
  2 个评论
Daniel Shub
Daniel Shub 2011-9-26
Again, have you looked at my answer? I am not sure why you don't want to use spectrogram, but I gave you a solution that doesn't require it.
Wayne King
Wayne King 2011-9-26
Irrespective of whether you want to write your own STFT or not, you are not correctly understanding what the size of the output matrix is.
Let's say you don't overlap your time segments (which will result in a blocky-looking spectogram, but ok), if you want to take the data 1,000 samples points at a time, then you will only have about 117 columns in your matrix. If your data is real-valued, then each column will only have 501 points, the positive frequencies of the fft() of each of your 1,000 sample segments.

请先登录，再进行评论。

raj
raj 2011-11-6
hello friends i too have the same question, i have an audio file and i want to implement stft manually but i have a problem with the overlap i implemented an overlap using if elseif branch but the window size doesnt change ...please help me with it i am pasting my code here
clear all
close all
clc
fs = 50e3;
%t = 1/fs;
w1 = 28000;
feed = w1/2;
load srtm
x1 = (srtm(8,1:10000));
nx = length(x1);% size of signal
w = hamming(w1)';% hamming window
nw = length(w);
T = (1:nx/fs); %size of window
pos=1;
n = 1;
%Null Matrices
W_total1 =[];%nullmatrix
while (pos+nw <= nx)
% while enough signal left
if n == 1
w1start = 0;
w1(end) = w1;
elseif n==2
w1start = (n-1)*w1 - feed;
w1(end) = (w1 + feed);
else
w1start = (w1start + w1(end))/2;
w1(end) = w1start + w1;
end
z(n,:) = x1(pos:pos+nw-1).*w(:,28000);
n = n+1;
pos = pos + nw/2;
W = fft(z(n-1,:));
W = W(1:end/2);
W_total1 = [W_total1 W.'];
end
figure;
imagesc(T,1:fs/2,10*log10(abs(W_total1/2e-5)))
axis xy; axis tight; colormap(jet);
ylim([0 1000])
xlabel('Time');
ylabel('Frequency (Hz)');
colorbar

类别

在 Help Center 和 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