# v_stftw

## PURPOSE

V_STFTW converts a time-domain signal into the time-frequency domain with the Short-time Fourier Transform [Y,SO,T,F]=(X,NW,M,OV)

## SYNOPSIS

function [y,so,tax,fax,gd]=v_stftw(x,nw,m,ov,nt)

## DESCRIPTION

    V_STFTW converts a time-domain signal into the time-frequency domain with the Short-time Fourier Transform [Y,SO,T,F]=(X,NW,M,OV)
      Usage:  (1) %%%%%% square-root-Hamming window with 50% overlap [Default options] %%%%%%
                  [y,so]=v_stftw(x,nw);                % default parameters: m='rM', ov=2
                  ...                                  % time-frequency domain processing goes here
                  ...                                  % y only contains positive frequencies, so size(y,2)=floor(1+nw/2)
                  z=v_istftw(y,so);                    % z is the same as x except for the first and last half-frames

              (2) %%%%%% Hann window with 75% overlap, zero-pad start/end of signal %%%%%%
                  [y,so]=v_stftw(x,nw,'en',4);         % Hanning window with ov=4
                  ...                                  % time-frequency domain processing goes here
                  z=v_istftw(y,so);                    % z is the same as x everywhere due to 'e' option

              (3) %%%%%% Plot Spectrogram if called without outputs %%%%%%
                  v_stftw(x,nw);                       % Plot spectrogram of first channel

              (4) %%%%%% Real-time STFT-domain processing %%%%%%
                  [y,so]=v_stftw([],nw,'pn',4);          % initialise v_stftw state variable, so
                  io=[];                               % initialise v_istftw state variable, io
                  while true                           % loop forever
                    % ... [acquire input chunk x]         % acquire new chunk of input samples
                    [y,so]=v_stftw(x,so,'p');           % 'p' option means more chunks to come
                    % ... [process in T-F domain]      % process signal in time-frequency domain
                    [z,io]=v_istftw(y,so,io);             % convert back to time domain
                    % ... [output chunk z];            % output chunk of processed samples
                  end

## CROSS-REFERENCE INFORMATION

This function calls:

- v_axisenlarge V_AXISENLARGE - enlarge the axes of a figure (f,h)
- v_cblabel V_CBLABEL add a label to a colorbar C=(L,H)
- v_rfft V_RFFT Calculate the DFT of real data Y=(X,N,D)
- v_windows V_WINDOWS Generate a standard windowing function (TYPE,N,MODE,P,OV)

This function is called by:

## SOURCE CODE

    0001 function [y,so,tax,fax,gd]=v_stftw(x,nw,m,ov,nt)
    0002 %V_STFTW converts a time-domain signal into the time-frequency domain with the Short-time Fourier Transform [Y,SO,T,F]=(X,NW,M,OV)
    0003 %  Usage:  (1) %%%%%% square-root-Hamming window with 50% overlap [Default options] %%%%%%
    0004 %              [y,so]=v_stftw(x,nw);                % default parameters: m='rM', ov=2
    0005 %              ...                                  % time-frequency domain processing goes here
    0006 %              ...                                  % y only contains positive frequencies, so size(y,2)=floor(1+nw/2)
    0007 %              z=v_istftw(y,so);                    % z is the same as x except for the first and last half-frames
    0008 %
    0009 %          (2) %%%%%% Hann window with 75% overlap, zero-pad start/end of signal %%%%%%
    0010 %              [y,so]=v_stftw(x,nw,'en',4);         % Hanning window with ov=4
    0011 %              ...                                  % time-frequency domain processing goes here
    0012 %              z=v_istftw(y,so);                    % z is the same as x everywhere due to 'e' option
    0013 %
    0014 %          (3) %%%%%% Plot Spectrogram if called without outputs %%%%%%
    0015 %              v_stftw(x,nw);                       % Plot spectrogram of first channel
    0016 %
    0017 %          (4) %%%%%% Real-time STFT-domain processing %%%%%%
    0018 %              [y,so]=v_stftw([],nw,'pn',4);          % initialise v_stftw state variable, so
    0019 %              io=[];                               % initialise v_istftw state variable, io
    0020 %              while true                           % loop forever
    0021 %                % ... [acquire input chunk x]         % acquire new chunk of input samples
    0022 %                [y,so]=v_stftw(x,so,'p');           % 'p' option means more chunks to come
    0023 %                % ... [process in T-F domain]      % process signal in time-frequency domain
    0024 %                [z,io]=v_istftw(y,so,io);             % convert back to time domain
    0025 %                % ... [output chunk z];            % output chunk of processed samples
    0026 %              end
    0027
    0028 %
    0029 %          (5) %%%%%% Optionally process a long signal in chunks %%%%%%
    0030 %              [y1,so]=v_stftw(x(1:1000),nw,'ep');      % convert chunk 1 of 3
    0031 %              [y2,so]=v_stftw(x(1001:2000),so,'p'); % Options other than 'p' will be remembered
    0032 %              [y3,so]=v_stftw(x(2001:end),so);      % Omit 'p' for final chunk (which can have x=[] if desired)
    0033 %              y=[y1; y2; y3];                       % Concatenate to get complete T-F data
    0034 %              [z1,io]=v_istftw(y(1:5,:),so);          % convert back to time domain
    0035 %              [z2,io]=v_istftw(y(6:10,:),so,io);    % convert back to time domain
    0036 %              z3=v_istftw(y(11:end,:),so,io);         % omit 'io' output on final chunk
    0037 %              z=[z1; z2; z3];                       % z is the same as x everywhere due to 'e' option
    0038 %
    0039 %  Inputs: x(t,...)    input signal (one signal per column)
    0040 %          nw       Window length (will be rounded up to a multiple of ov)
    0041 %                       alternatively, the so output from the previous chunk's call to v_stftw
    0042 %          m        mode string including window code
    0043 %          ov       integer overlap factor; the frame hop is nw/ov. [2]
    0044 %          nt       Optional DFT length (normally >=nw). [Default is nw as modified by 'i','I','u' options]
    0045 %
    0046 % Outputs: y(tf,k,...)  output STFT (frame=tf, freq=k). Number of frequencies is normally floor(1+nw/2) from DC to Nyquist
    0047 %          so           structure (see below for details) used in the inverse transformation, v_istftw(), and also as the
    0048 %                       nw argument in a subsequent call to v_stftw() when processing a long signal in multiple chunks.
    0049 %          tax          sample numbers corresponding to the centre of each frame; divide by fs to get times
    0050 %          fax          normalized centre frequencies of each bin (multiply by fs for actual frequencies)
    0051 %          gd(tf,k,...) group delay in samples. A impulse at x(i,n+1) will give gd(i,:)=n. Will be infinite wherever y()=0.
    0052 %
    0053 % The mode string, m, is a character string containing a combination of the following options and a window-code (see below):
    0054 %
    0055 %             'p'   x is a partial signal which will be followed by another chunk (prevents paddding of the final frame)
    0056 %
    0057 %             'r'   pad final frame with reflected data if necessary [default]
    0058 %             'z'   pad final frame with zeros if necessary
    0059 %             't'   truncate data to an exact number of frames
    0060 %             'e'   pad with samples at beginning and end of signal to include all frames with at least one signal sample
    0061 %
    0062 %             'i'   zero-pad by a factor of 2 to increase apparent frequency resolution
    0063 %             'I'   zero-pad by a factor of 4 to increase apparent frequency resolution
    0064 %             'u'   zero-pad to round transform length up to a power of 2 to speed computation
    0065 %
    0066 %             's'   plot spectrogram of first channel (default if no outputs)
    0067 %             'S'   plot mean spectra and mean power waveforms of all channels
    0068 %
    0069 %            window-code   Window   mode-string  overlaps   Sidelobe  3dB-BW  6dB-BW  Falloff    Comment
    0070 %
    0071 %                'c'      cos(1)        s         2,3,5      -23dB     1.2     1.6  -12dB/oct   used in MP3
    0072 %                'k'      kaiser(5)     boq         2        -23dB     1.2     1.7   -6dB/oct   used in AAC
    0073 %                'm'      hamming       s         3,4,5      -43dB     1.3     1.8   -6dB/oct   low sidelobes [default if ov>2]
    0074 %                'M'      hamming       sq        2,3,5      -24dB     1.1     1.5   -6dB/oct   sqrt Hamming  [default if ov=2]
    0075 %                'n'      hann          s         3,4,5      -31dB     1.4     2.0  -18dB/oct   rapid falloff
    0076 %                'R'      rectangle                 1        -13dB     0.9     1.2   -6dB/oct   narrow bandwidth [default if ov=1]
    0077 %                'v'      vorbis        s           2        -21dB     1.3     1.8  -18dB/oct   used in Vorbis
    0078 %                'V'      rsqvorbis     sq          2        -26dB     1.1     1.5   -6dB/oct   sqrt of raised and squared Vorbis
    0079 %
    0080 %                 where mode-string is the mode input to the v_windows call. For perfect reconstruction,
    0081 %                 the ov input must be a multiple of one of the numbers listed in the overlaps column.
    0082 %                 To see more details of a specific window, execute:
    0083 %                      v_windows(window-code,55440,mode-string,param);   % Note: use window-code='w' for rsqvorbis
    0084 %
    0085 % The contents of the so output structure are:
    0086 %
    0087 %           so.mx           cumulative number of input samples
    0088 %           so.mf           cumulative number of full frames output so far (excluding 'e' option prefix frames)
    0089 %           so.nt           length of Fourier transform
    0090 %           so.nh           frame hop in samples
    0091 %           so.wa(1,nw)     analysis window (e.g. Rectangular or Hamming)
    0092 %           so.po           index of final-frame padding option: 1='z', 2='r', 3='t'
    0093 %           so.me           true iff 'e' option given i.e. signal is zero-padded at start and end
    0094 %           so.ff           true iff no frames have yet been output
    0095 %           so.mp           true iff 'p' option given i.e. there are additional chunks to follow
    0096 %
    0097 % Notes:
    0098 %
    0099 %  (1) The scaling preserves power in a 2-sided transform so that
    0100 %
    0101 %        p = mean(x.^2) ~=~ mean(sum(abs(yy).^2,2)) where yy=[y conj(y(:,nw-size(y,2)+1:-1:2))]
    0102 %
    0103 %      where yy is the double-sided spectrum. This equivalence is stochastic rather than exact.
    0104 %  (2) For perfect reconstruction (see [3]), the overlap factor, ov, must be an integer multiple
    0105 %      of one of the values listed above in the ov column for the selected window.
    0106 %  (3) Using the default parameters, perfect reconstruction will not be obtained for the first and
    0107 %      the last frames of the signal. Using the 'e' option will add padding frames to the start and
    0108 %      end of the signal so that the entire signal will be reconstructed perfectly.
    0109 %  (4) By default, padding at the start and end of the signal will avoid introducing a discontinuity
    0110 %      by using a time-reversed reflection of the input signal. Alternatively, the 'z' option will
    0111 %      pad with zeros and the 't' option will truncate the signal to fit into an exact number of frames.
    0112 %  (5) To avoid temporal aliasing in the time-frequency domain , an overlap factor of ov>=4 is required for
    0113 %      most windows. Using an overlap factor of ov=2 (the default) halves the number of frames to process
    0114 %      but may result in aliasing for some signals although this is normally slight (see [1]).
    0115 %
    0116 % Refs: [1] J. B. Allen. Short term spectral analysis, synthesis, and modification by discrete Fourier transform.
    0117 %           IEEE Trans. Acoustics, Speech and Signal Processing, 25 (3): 235–238, June 1977.
    0118 %       [2]    R. Crochiere. A weighted overlap-add method of short-time fourier analysis/synthesis.
    0119 %           IEEE Trans. Acoustics, Speech and Signal Processing, 28 (1): 99–102, 1980.
    0120 %       [3]    J. Princen and A. Bradley. Analysis/synthesis filter bank design based on time domain aliasing cancellation.
    0121 %           IEEE Trans. Acoustics, Speech and Signal Processing, 34 (5): 1153–1161, 1986. doi: 10.1109/TASSP.1986.1164954.
    0122 %
    0123 % Revision History
    0124 %    2022-04-10 First version
    0125 persistent nw0 nt0 wi0 w0 wt wm wp
    0126 if isempty(nw0)
    0127     nw0=0;
    0128     wt={'rectangle','hamming','hamming','cos','kaiser','hann','rsqvorbis','vorbis'};
    0129     wm={'','sq','s','s','boq','s','sq','s'};
    0130     wp=[0 0 0 1 5 0 0 0];
    0131 end
    0132 sx=size(x);
    0133 if sx(1)==1 && length(sx)==2 % transpose x if it is a row vector
    0134     x=x.';
    0135     sx=size(x);
    0136 end
    0137 nx=sx(1);                   % number of samples
    0138 nc=prod(sx(2:end));         % total number of channels
    0139 x=reshape(x,nx,nc);         % make input data 2-dimensional
    0140 if nargin<3 || ~numel(m)    % no mode options specified
    0141     m='';                   % use default modes
    0142 end
    0143 if isstruct(nw)             % nw is structure from previous call to v_stftw
    0144     so=nw;                  % replicate into so since most fields will be identical
    0145     wa=nw.wa;               % analysis window
    0146     nh=nw.nh;               % frame hop (samples)
    0147     po=nw.po;               % padding option
    0148     nt=nw.nt;               % overwrite nt with length of transform
    0149     me=nw.me;               % 'e' option given
    0150     ff=nw.ff;               % first frame flag
    0151     mf=nw.mf;
    0152     if ~isempty(nw.xx)        % there is a data tail from a previous call
    0153         x=cat(1,nw.xx,x);     % pre-append saved data
    0154     end
    0155     so.mx=so.mx+nx;          % cumulative signal length (using original value of nx)
    0156     sx=size(x);             % recalculate sx (especially necessary if input argument x=[]
    0157     nx=sx(1);               % recalculate the number of samples
    0158     nw=length(wa);          % overwrite nw with the window length
    0159 else
    0160     if nargin<4 || ~numel(ov)
    0161         ov=2;               % default overlap factor
    0162     else
    0163         ov=round(ov);         % force integer overlap factor
    0164     end
    0165     nh=ceil(nw/ov);         % force integer hop length
    0166     nw=ov*nh;               % recalculate DFT length
    0167     wi=find(sum(repmat(m(:),1,8)==repmat('RMmcknVv',length(m),1),1),1); % check if m specifies the window
    0168     if isempty(wi)          % if no window specified
    0169         wi=1+(ov>1)+(ov>2); % use 'R', 'M' or 'm'
    0170     end
    0171     % sort out zero-padding
    0172     if nargin<5             % DFT length is not explicitly specified
    0173         nt=nw;                  % default transform length
    0174         if any(m=='i')          % 'i' = multiply nt by 2 for increased frequency resolution
    0175             nt=2*nt;
    0176         end
    0177         if any(m=='I')          % 'I' = multiply nt by 4 for increased frequency resolution
    0178             nt=4*nt;
    0179         end
    0180         if any(m=='u')          % 'u' = round up nt to a power of 2 for increased frequency resolution
    0181             [fw,pw]=log2(nt);
    0182             nt=pow2(pw-(fw==0.5));
    0183         end
    0184     end
    0185     me=any(m=='e');
    0186     if nw==nw0 && wi==wi0 && nt==nt0    % check if window already cached
    0187         wa=w0;                          % use cached version if available
    0188     else
    0189         wa = v_windows(wt{wi},nw,wm{wi},wp(wi)); % only need the parameter for 'c' and 'k' windows
    0190         wa=wa/sqrt(sum(wa.^2)*nt);      % scale to preserve power
    0191         nw0=nw;                         % save parameters of cached window
    0192         nt0=nt;
    0193         wi0=wi;
    0194         w0=wa;                          % save window in cache
    0195     end
    0196     % sort out enframing details
    0197     ff=1;                               % flag to indicate the first frame
    0198     po=find(sum(repmat(m(:),1,3)==repmat('zrt',length(m),1),1),1); % find padding option
    0199     if isempty(po)
    0200         po=2;                           % default to po=2 ('r')
    0201     end
    0202     % create so structure
    0203     mf=0;
    0204     so.wa=wa;                       % save analysis window
    0205     so.po=po;                       % save padding option for final frame
    0206     so.nh=nh;                       % save frame hop
    0207     so.nt=nt;                       % length of transform
    0208     so.me=me;                       % 'e' option given
    0209     so.mx=nx;                       % cumulative signal length
    0210     so.mf=0;                        % zero cumulative frames output so far (excluding prefix frames)
    0211 end
    0212 % now do the STFT
    0213 mp=any(m=='p');                     % 'p' option given: this is not the last chunk
    0214 f0=(1-nw/nh)*(me && ff);            % initial frame: either 0 or (1-ov)
    0215 if nx>0                             % there are some samples
    0216     if me && ~mp                    % include all frames that contain at least one sample
    0217         f1=ceil(nx/nh)-1;           % last frame that include at least one sample
    0218     elseif po<3 && ~mp              % ensure all samples are in at least one frame
    0219         f1=max(ceil((nx-nw)/nh),0);
    0220     else                            % only include full frames
    0221         f1=floor((nx-nw)/nh);               % last full frame (might be negative and/or <f0)
    0222     end
    0223 else                                        % input signal is empty
    0224     f1=-1;
    0225 end
    0226 nf=max((f1-f0+1)*(f1>=0),0);                % number of frames to output
    0227 if nf>0                                     % if there are frames to output
    0228     ff=0;                                   % turn off first frame flag
    0229     indx=repmat(nh*(f0:f1).',nw,1)+reshape(repmat(0:nw-1,nf,1),nf*nw,1);    % index into x
    0230     mk=indx<0 | indx>nx-1;                  % mask for values outside 0:nx-1
    0231     indx=nx+0.5-abs(mod(indx,2*nx)+0.5-nx); % reflect values outside the range
    0232     v=x(indx,:);                            % create output
    0233     if po==1                                % padding option was 'z'
    0234         v(mk,:)=0;                          % set padding values to zero
    0235     end
    0236     y=v_rfft(reshape(v,[nf,nw,sx(2:end)]).*repmat(wa,[nf,1,sx(2:end)]),nt,2);  % perform windowed DFT
    0237 else
    0238     y=zeros([0,nw,sx(2:end)]);
    0239 end
    0240 if nargout>1                                % need additional outputs
    0241     if nf>0
    0242         tax=(f0+mf:f1+mf)*nh+0.5*(1+nw);    % frame centre times (where input samples are 1,2,...)
    0243     else
    0244         tax=[];
    0245     end
    0246     fax=(0:size(y,2)-1)/nt;                 % bin centre frequencies (normalized by sample frequency)
    0247     if nf>0
    0248         so.xx=x((f1+1)*nh+1:end,:);            % save data tail for next chunk (starts at sample nh+1 of the last frame in y)
    0249     else
    0250         so.xx=x;                            % if nf==0 then data tail is entire signal (including any previous tail)
    0251     end
    0252     so.ff=ff;                               % first frame flag
    0253     so.mp=mp;                               % flag indicating more chunks to follow
    0254     so.mf=mf+nf+f0;                         % add number of additional non-prefix frames output
    0255     if nargout>4                            % need group delay output
    0256         if nf>0                             % if there are frames to output
    0257             gd=real(v_rfft(reshape(v,[nf,nw,sx(2:end)]).*repmat(wa.*(0:nt-1),[nf,1,sx(2:end)]),nt,2)./y);  % calculate group delay
    0258         else
    0259             gd=zeros([0,nw,sx(2:end)]);
    0260         end
    0261     end
    0262 end
    0263 if nf>0 && (~nargout || any(lower(m)=='s')) % no outputs or 'sS' options so plot graphs
    0264     ta=(f0+mf:f1+mf)*nh+0.5*(1+nw);
    0265     nv=size(y,2);
    0266     fa=(0:nv-1)/nt;
    0267     clf;
    0268     if any(m=='S')
    0269         subplot(2,1,1);
    0270         yp=abs(reshape(y,[nf,nv,nc])).^2;   % calculate power spectrum
    0271         plot(ta,10*log10(reshape(sum(cat(2,yp,yp(:,2:floor((nt+1)/2),:)),2),nf,nc))); % add in power from negative frequencies
    0272         v_axisenlarge([-1 -1.05]);
    0273         xlabel('Sample Number');
    0274         ylabel('Frame power (dB)');
    0275         subplot(2,1,2);
    0276         plot(fa,10*log10(reshape(mean(yp,1),nv,nc)));
    0277         v_axisenlarge([-1 -1.05]);
    0278         xlabel('Normalized Frequency');
    0279         ylabel('Mean Power (dB)');
    0280     else
    0281         imagesc(ta,fa,20*log10(abs(y(:,:,1)))'); % plot only the first channel
    0282         colorbar;
    0283         axis 'xy';
    0284         ylabel('Normalized Frequency');
    0285         xlabel('Sample Number');
    0286         v_cblabel('Energy (dB)');
    0287         hold on;
    0288         clim=get(gca,'clim');
    0289         set(gca,'clim',[max(clim(1),clim(2)-40) clim(2)]);
    0290         ylim=get(gca,'ylim');
    0291         plot([0.5;0.5],ylim',':k');             % indicate start of signal
    0292         if ~mp
    0293             plot([1;1]*(so.mx+0.5),ylim',':k'); % indicate end of signal
    0294         end
    0295         hold off
    0296     end
    0297 end
    0298