MATLAB: How to get surface data from spectrogram

3d printingspectrogramstlsurface data

I currently have a spectrogram that is generated from an audio file and I am struggling to figure out how to get the surface data.

I would like to get the surface data of the spectrogram in order to convert it to an STL file using STLwrite or surf2STL. Is there anyway to get the X,Y,Z data from the spectrogram? Keep in mind, the X,Y,Z data would change when the program uses another audio file to generate another spectrogram. Any help is appreciated.

Here is the audio file used in the program. https://drive.google.com/file/d/1d__c9yvfezmhJakezxnVsj7OoDoXkFa_/view?usp=sharing

% import Matlab's Hallelujah sound file
% load handel
% filename = 'Halleluja.wav';
% audiowrite(filename,y,Fs);
% clear y Fs;
% [sample_data, sample_rate] = audioread(filename);
%uncomment out the next two lines if you want to import your own sound file
filename = 'Halleluja.wav'
[sample_data, sample_rate] = audioread(filename');
info = audioinfo(filename)
%listen to unfiltered sound file
sound(sample_data, sample_rate)
%stop sound
clear sound
%calculate the sampling times
sample_period = 1/sample_rate;                              %calculate the time between samples
%samples = [1,2*sample_rate];                                %modify number of samples
%[sample_data, sample_rate] = audioread(filename,samples);   %re-import sound file using modified rate
%remove the right channel of audio data
sample_data_width = size(sample_data);
if(sample_data_width(2) == 2)
    sample_data(:,1) = [];  
end
%Generate the sampling times
time = 0:sample_period:(length(sample_data)-1)/sample_rate;
%Plot the Time Domain Representation of Unfiltered Sound
plot(time,sample_data)
height_multiplier = 1.6;
axis([0,11,-1/height_multiplier,1/height_multiplier])
title('Time Domain Representation - Unfiltered Sound')
xlabel('Time (seconds)')
ylabel('Amplitude')
xlim([0 time(end)])
%Plot the Frequency Domain Representation of Unfiltered Sound
m = length(sample_data);                                    %Calculate the original sample length
n = pow2(nextpow2(m));                                      %Transform the length so the number of samples is a power of 2
y = fft(sample_data, n);                                    %Calculate the Fast Fourier Transform
frequencies = (0:n-1)*(sample_rate/n);                      %Generate the frequencies 
amplitudes = abs(y)/n;                                      %Generate the amplitudes
%Plot the Frequency Domain
plot(frequencies(1:floor(n/2)), amplitudes(1:floor(n/2)))
title('Frequency Domain Representation - Unfiltered Sound')
xlabel('Frequency')
ylabel('Amplitude')
width_slider2 = 1;
xlim([0 frequencies(end)/width_slider2])
%Filtered Sound
%Filter out the noise
order = 7;
[b,a] = butter(order,1000/(sample_rate/2),'low');
filtered_sound = filter(b,a,sample_data);
sound(filtered_sound, sample_rate)                          %Listen the unfiltered sound file
clear sound                                                 %Stop sound
%Plot the Time Domain Representation - Filtered Sound 
%Generate the new sampling times
time2 = (0:sample_period:(length(filtered_sound)-1)/sample_rate);
%Plot the Time Domain Representation
plot(time2,filtered_sound)
height_multiplier2 = 1;
axis([0,11,-1/height_multiplier2,1/height_multiplier2])
title('Time Domain Representation - Filtered Sound')
xlabel('Time (seconds)')
ylabel('Amplitude')
xlim([0 time2(end)])
%Plot the Frequency Domain Representation - Filtered Sound 
m1 = length(sample_data);                                   % Calculate the original sample length
n1 = pow2(nextpow2(m1));                                    % Transform the length so the number of sampels is a power of 2.
y1 = fft(filtered_sound, n1);                               % Calculate the Fast Fourier Transform  
frequencies2 = (0:n1-1)*(sample_rate/n1);                   % Generate the frequencies
amplitudes2 = abs(y1)/n1;                                   % Generate the amplitudes
%Plot the Frequency Domain Representation
plot(frequencies2(1:floor(n1/2)),amplitudes2(1:floor(n1/2)))
title('Frequency Domain Representation - Filtered Sound')
xlabel('Frequency')
ylabel('Amplitude')
width_slider2 = 11;
xlim([0 frequencies2(end)/width_slider2])
%3D Spectrogram - Filtered Sound
win = 128;                                                  % Window length in samples                                                          
hop = win/2;                                                % Number of samples between overlapping windows
nfft = win;                                                 % Width of each frequency bin
spectrogram(filtered_sound,win,hop,nfft,sample_rate,'yaxis');
view(-45,65)

Best Answer

You have to learn to read the help and documentation of the functions used. If you try:

help spectrogram

This is what will be returned:

 SPECTROGRAM Spectrogram using a Short-Time Fourier Transform (STFT).
    S = SPECTROGRAM(X) returns the short-time Fourier transform of the
    signal specified by vector X in the matrix S. By default, X is divided
    into eight segments with 50% overlap, and each segment is windowed with
    a Hamming window. The number of frequency points used to calculate the
    discrete Fourier transforms is equal to the larger of 256 or the next
    power of two greater than the segment length.
 
    If X cannot be divided exactly into eight segments, X will be
    truncated.
 
    S = SPECTROGRAM(X,WINDOW), when WINDOW is a vector, divides X into
    segments of the same length as WINDOW, and then windows each segment
    with the vector specified in WINDOW.  If WINDOW is an integer, the
    function divides X into segments of length equal to that integer value
    and windows each segment with a Hamming window.  If WINDOW is not
    specified, the default is used.
 
    S = SPECTROGRAM(X,WINDOW,NOVERLAP) specifies NOVERLAP samples of
    overlap between adjoining segments. NOVERLAP must be an integer smaller
    than WINDOW if WINDOW is an integer.  NOVERLAP must be an integer
    smaller than the length of WINDOW if WINDOW is a vector.  If NOVERLAP
    is not specified, the default value is used to obtain a 50% overlap.
 
    S = SPECTROGRAM(X,WINDOW,NOVERLAP,NFFT) specifies the number of
    frequency points used to calculate the discrete Fourier transforms.
    If NFFT is not specified, the default NFFT is used.
 
    S = SPECTROGRAM(X,WINDOW,NOVERLAP,NFFT,Fs) specifies the sample rate,
    Fs, in Hz. If Fs is specified as empty, it defaults to 1 Hz. If it is
    not specified, normalized frequency is used.
 
    Each column of S contains an estimate of the short-term, time-localized
    frequency content of X.  Time increases across the columns of S, from
    left to right.  Frequency increases down the rows, starting at 0.  If X
    is a length NX complex signal, S is a complex matrix with NFFT rows and
    k = fix((NX-NOVERLAP)/(length(WINDOW)-NOVERLAP)) columns. For real X, S
    has (NFFT/2+1) rows if NFFT is even and (NFFT+1)/2 rows if NFFT is
    odd.
 
    [S,F,T] = SPECTROGRAM(...) returns a vector of frequencies, F, and a
    vector of times, T, at which the spectrogram is computed. F has length
    equal to the number of rows of S. T has length k (defined above) and
    its value corresponds to the center of each segment. If a sample
    rate is not provided, F contains normalized frequencies.
 
    [S,F,T] = SPECTROGRAM(X,WINDOW,NOVERLAP,F) computes the two-sided
    spectrogram at the normalized frequencies specified in the vector F. F
    must have at least two elements.
 
    [S,F,T] = SPECTROGRAM(X,WINDOW,NOVERLAP,F,Fs) computes the two-sided
    spectrogram at the frequencies specified in vector F.  F must be
    expressed in Hz and have at least two elements.
 
    [S,F,T,P] = SPECTROGRAM(...) P is a matrix representing the Power
    Spectral Density (PSD) of each segment. For real signals, SPECTROGRAM
    returns the one-sided modified periodogram estimate of the PSD of each
    segment; for complex signals and in the case when a vector of
    frequencies is specified, it returns the two-sided PSD.
 
    [S,F,T,P] = SPECTROGRAM(...,'MinThreshold',THRESH) sets the elements of
    P to zero when the corresponding elements of 10*log10(P) are less than
    THRESH. Specify THRESH in decibels. The default value of THRESH is
    -Inf.
 
    [S,F,T,P] = SPECTROGRAM(...,'reassigned') reassigns each PSD estimate
    to the location of its center of gravity.  The reassignment is done
    in-place and returned in P.
 
    [S,F,T,P,Fc,Tc] = SPECTROGRAM(...) returns the locations in frequency
    and time of the center of gravity of each estimate in the spectrogram.
    The frequencies and times are returned in matrices, Fc, and Tc,
    respectively.  Fc and Tc have the same dimensions as the spectrogram, S.
 
    [...]  = SPECTROGRAM(...,SPECTRUMTYPE) uses the window scaling
    algorithm specified by SPECTRUMTYPE when computing the power spectral
    density matrix P.
    SPECTRUMTYPE can be set to 'psd' or 'power':
       'psd'   - returns the power spectral density.
       'power' - scales each estimate of the PSD by the equivalent noise
                bandwidth of the window (in Hz).  Use this option to
                obtain an estimate of the power at each frequency.
    The default value for SPECTRUMTYPE is 'psd'.
 
    [...] = SPECTROGRAM(...,FREQRANGE)  returns the PSD over the specified
    range of frequencies based upon the value of FREQRANGE:
 
       'onesided' - returns the one-sided matrix P of a real input signal X.
          If NFFT is even, P has NFFT/2+1 rows and is computed over the
          interval [0,pi].  If NFFT is odd, then P has (NFFT+1)/2 rows
          and is computed over the interval [0,pi). When Fs is specified,
          the intervals become [0,Fs/2] and [0,Fs/2) for even and odd NFFT,
          respectively.
 
       'twosided' - returns the two-sided matrix P for either real or complex
          input X.  P has NFFT rows and is computed over the interval
          [0,2*pi). When Fs is specified, the interval becomes [0,Fs).
 
       'centered' - returns the centered two-sided matrix P for either real
          or complex X.  P has NFFT rows and is computed over the interval
          (-pi,pi] for even length NFFT and (-pi,pi) for odd length NFFT.
          When Fs is specified, the intervals become (-Fs/2, Fs/2] and
          (-Fs/2,Fs/2) for even and odd NFFT, respectively.
 
       FREQRANGE may be placed in any position in the input argument list
       after NOVERLAP.  The default value of FREQRANGE is 'onesided' when X
       is real and 'twosided' when X is complex.
 
    [...] = SPECTROGRAM(...,'OutputTimeDimension',TIMEDIMENSION) specifies
    the orientation of S, T, P, Fc, and Tc according to the location of
    the time dimension. If TIMEDIMENSION is set to 'downrows', the time
    dimension of S, P, Fc, and Tc is down the rows and the frequency
    dimension is across the columns. T is returned as a column vector. If
    TIMEDIMENSION is set to 'acrosscolumns', the time dimension of S, P,
    Fc, and Tc is across the columns and the frequency dimension is down
    the rows. T is returned as a row vector. This argument is ignored if
    this function is called with no output arguments. The default value is
    'acrosscolumns'.
 
    SPECTROGRAM(...) with no output arguments plots the PSD estimate for
    each segment on a surface in the current figure. 
 
    SPECTROGRAM(...,FREQLOCATION) controls where MATLAB displays the
    frequency axis on the plot. This string can be either 'xaxis' or
    'yaxis'.  Setting this FREQLOCATION to 'yaxis' displays frequency on
    the y-axis and time on the x-axis.  The default is 'xaxis' which
    displays the frequency on the x-axis. If FREQLOCATION is specified when
    output arguments are requested, it is ignored.
 
    EXAMPLE 1: Spectrogram of quadratic chirp
      t=0:0.001:2;                    % 2 secs @ 1kHz sample rate


      y=chirp(t,100,1,200,'q');       % Start @ 100Hz, cross 200Hz at t=1sec


      spectrogram(y,kaiser(128,18),120,128,1E3,'yaxis');
      title('Quadratic Chirp: start at 100Hz and cross 200Hz at t=1sec');
 
    EXAMPLE 2: Reassigned spectrogram of quadratic chirp
      t=0:0.001:2;                    % 2 secs @ 1kHz sample rate
      y=chirp(t,100,1,200,'q');       % Start @ 100Hz, cross 200Hz at t=1sec
      spectrogram(y,kaiser(128,18),120,128,1E3,'reassigned','yaxis');
      title('Quadratic Chirp: start at 100Hz and cross 200Hz at t=1sec');
 
    EXAMPLE 3:  Plot instantaneous frequency of quadratic chirp
      t=0:0.001:2;                    % 2 secs @ 1kHz sample rate
      y=chirp(t,100,1,200,'q');       % Start @ 100Hz, cross 200Hz at t=1sec
      % remove estimates less than -30 dB
      [~,~,~,P,Fc,Tc] = spectrogram(y,kaiser(128,18),120,128,1E3,'minthreshold',-30);
      plot(Tc(P>0),Fc(P>0),'. ')
      title('Quadratic Chirp: start at 100Hz and cross 200Hz at t=1sec');
      xlabel('Time (s)')
      ylabel('Frequency (Hz)')
 
    EXAMPLE 4: Waterfall display of the PSD of each segment of a VCO
      Fs = 10e3;
      t = 0:1/Fs:2;
      x1 = vco(sawtooth(2*pi*t,0.5),[0.1 0.4]*Fs,Fs);
      spectrogram(x1,kaiser(256,5),220,512,Fs);
      view(-45,65)
      colormap bone
 
    See also PERIODOGRAM, PWELCH, GOERTZEL.
    Documentation for spectrogram
       doc spectrogram
    Other functions named spectrogram
       gpuArray/spectrogram    tall/spectrogram

There you see what to do and what variables to use for plotting your surface-plot of the spectrogramme:

[S,F,T,P] = spectrogram(X,WINDOW,NOVERLAP,NFFT,Fs);
surf(T,F,log10(P)),shading flat

Adjust caxis and other niceties to make the plot the way you want it.

HTH

Related Solutions

MATLAB: Frame by frame freqeuncy analysis of voice.

Have a look at the documentation for the spectrogram function, that should do what you want. 0.025 is a rather short time so perhaps you can get something from calculating the instantaneous frequency - so have a look at that (and first the Hilbert transform too), this is likely a rough suggestion since we expect frequencies way higher than 20 Hz.

HTH

MATLAB: Manual spectrogram creation without using spectrogram command

hello

see below (myspecgram)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
[data,Fs] = audioread('myWAVaudiofile.wav');
channel = 1;
signal = data(:,channel);
samples = length(signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 4096;    % 
OVERLAP = 0.75;
% spectrogram dB scale
spectrogram_dB_scale = 100;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 2;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [sg,fsg,tsg] = specgram(signal,NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  
[sg,fsg,tsg] = myspecgram(signal,NFFT,Fs,OVERLAP);  
% FFT normalisation and conversion amplitude from linear to dB (peak)
sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only
% saturation of the dB range : 
% saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;
% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([' Spectrogram / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
function  [sg,freq_vector,time] = myspecgram(signal,nfft,Fs,Overlap)
% [sg,freq_vector,time] = myspecgram(signal,NFFT,Fs,hanning(NFFT),floor(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)
dt = 1/Fs;
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
    sg = [];
    for i=1:spectnum
        start = (i-1)*offset;
        sw = signal((1+start):(start+nfft)).*window;
        sg = [sg (abs(fft(sw))*4/nfft)];     % X=fft(x.*hanning(N))*4/N; % hanning only 
        time(i) = (start+nfft/2)*dt; % time defined as in the middle of the data buffer
    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
    
sg = sg(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

Best Answer

Related Solutions

MATLAB: Frame by frame freqeuncy analysis of voice.

MATLAB: Manual spectrogram creation without using spectrogram command

Related Question