Manual Implementation of STFT of an audio signal.

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.

Answer 1

[Wayne King]

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.

Answer 2

[Wayne King]

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.

[UJJWAL]

Hi,

Thanks for your response.

The problem is this :-

The total length of the samples is 117424.

For each sample I will find the 65536-point FFT as its sampling frequency is 44100 Hz.

This amounts to a 65536 x 117424 matrix which is an extremely large matrix and gives me an out of memory error.

Hence I decided to take 1000 samples at a time which would need a 65536 x 1000 matrix , plot it, then clear the matrix and take the next 1000 segments and so on, while holding the plot. But this idea did not work. I am not overlapping the windows because in my work , I am required to find the STFT without the overlapping of the windows.

[Walter Roberson]

Each of your 65536 x 117424 matrices is 500 megabytes. How exactly are you plotting that matrix? Unless you are somehow rendering it down to a *much* smaller matrix, the graphics for each of those plots is going to require at least 500 megabytes: handle graphics hold on to the original data values until the values are replaced or the graphic is removed. You are thus attempting to use the graphics as temporary storage at 500+ megabytes each, 117 1/2 times. Where is that 59 gigabytes of graphics memory going to come from? MATLAB does not have access to memory beyond what is normally available to your program.

[UJJWAL]

ok Walter.. so given the specifications of the audio file .. how do I find its spectrogram without using the in-built function ?? Any suggestions ?

[Walter Roberson]

What output are you hoping to generate out of the aggregate data? For example were you hoping to "line up" the spectrograms one after another on the same baseline? If so, how wide did you say your screen was?

Could you post a link to the graph you are creating out of a single group of 1000 samples?

http://www.mathworks.com/matlabcentral/answers/7924-where-can-i-upload-images-and-files-for-use-on-matlab-answers

[UJJWAL]

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]

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.

Answer 3

[Daniel Shub]

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)

Answer 4

[Wayne King]

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)');

Answer 5

[UJJWAL]

Actually I do not want to use the spectrogram function. I want to write my own code for finding the spectrogram

[Daniel Shub]

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]

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.

Answer 6

[raj]

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