tomwalters@0: % support file for 'aim-mat'
tomwalters@0: %
tomwalters@0: % This external file is included as part of the 'aim-mat' distribution package
bleeck@3: % (c) 2011, University of Southampton
bleeck@3: % Maintained by Stefan Bleeck (bleeck@gmail.com)
bleeck@3: % download of current version is on the soundsoftware site: 
bleeck@3: % http://code.soundsoftware.ac.uk/projects/aimmat
bleeck@3: % documentation and everything is on http://www.acousticscale.org
bleeck@3: 
tomwalters@0: 
tomwalters@0: function f=spectrogram(sig,nr_f,maxfre)
tomwalters@0: 
tomwalters@0: nsig=changesr(sig,maxfre*2);
tomwalters@0: sr=getsr(nsig);
tomwalters@0: vals=getvalues(nsig);
tomwalters@0: 
tomwalters@0: % sig_len=getlength(sig);
tomwalters@0: % time_step=sig_len/nr_t;
tomwalters@0: % window=[0:time_step:sig_len];
tomwalters@0: % win_len=getnrpoints(nsig)/nr_t;
tomwalters@0: % window = hann(round(win_len));
tomwalters@0: % noverlap=round(win_len/2);
tomwalters@0: noverlap=nr_f*2-2;
tomwalters@0: [Bb,freqs,ts]=specgram(vals,nr_f*2,sr*2,[],noverlap);
tomwalters@0: 
tomwalters@0: B=Bb(2:nr_f+1,:);
tomwalters@0: 
tomwalters@0: nr_t=size(B,2);
tomwalters@0: nr_f=size(B,1);
tomwalters@0: 
tomwalters@0: srf=nr_t/getlength(sig);
tomwalters@0: f=field(nr_t,nr_f,srf);
tomwalters@0: f=setvalues(f,abs(B));
tomwalters@0: f=setmaxfre(f,maxfre);
tomwalters@0: 
tomwalters@0: return
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % function f=spectrogram(sig,tstart,tstop,nr_t,nr_f,maxfre)
tomwalters@0: % alle Zeiten in Sekunden
tomwalters@0: % macht ein Sektrogramm aus dem Signal sig, das zu nr_t Zeiten gesampelt
tomwalters@0: % wird und nr_f Frequenzen hat
tomwalters@0: % das Signal wird aus sig in der Zeit von start bis stop genommen
tomwalters@0: % das Signal sig wird erst so auf eine Samplerate gebracht, dass 
tomwalters@0: % das Resultat der FFT eine Auflösung bis maxfre hat.
tomwalters@0: 
tomwalters@0: if nargin<6
tomwalters@0:     maxfre=GetSR(sig)/2;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % die alte (orginale) Samplerate
tomwalters@0: sr_old=GetSR(sig);
tomwalters@0: 
tomwalters@0: %wenn die maxfre anders ist als die, die das Signal mitbringt, dann muss das Signal heruntergesamplet werden:
tomwalters@0: if sr_old~=maxfre*2
tomwalters@0:     sr_neu=maxfre*2;
tomwalters@0:     sig=changesr(sig,sr_neu);   % änder das ganze Signal zu der neuen SR
tomwalters@0: else 
tomwalters@0:     sr_neu=sr_old;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % Punkte, die die FFT benötigt, damit nr_f punkte hinten raus kommen
tomwalters@0: nr_fft=nr_f*2;
tomwalters@0: % und das ist ein Signal von dieser Länge: (auch mit neuer Samplerate)
tomwalters@0: stime=bin2time(nr_fft,sr_neu);
tomwalters@0: % Dauer des zu untersuchenden Signals
tomwalters@0: duration=tstop-tstart;
tomwalters@0: % wenn die gewünschte Samplingrate nicht hoch genug ist, weil zuviel Punkte erforderlich sind
tomwalters@0: if duration < stime
tomwalters@0:     error('spectrogram: desired Samplingrate too low: take less points or higher freqeuncy!');
tomwalters@0:     return;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % meine Notation ist anders als Matlab: Bei mir ist x die Zeit
tomwalters@0: % und y die Frequenz
tomwalters@0: sr_field=nr_t/duration; % die "SampleRate" des Feldes (sehr klein, da nur wenig Punkte)
tomwalters@0: f=field(nr_t,nr_f,sr_field);
tomwalters@0: f=setmaxfre(f,maxfre);  % Für die Grafik: Wie hoch die Frequenzen geben
tomwalters@0: 
tomwalters@0: % damit die fft immer von einem richtigen vollen Signal gemacht werden kann, werden die 
tomwalters@0: % FFTs an Stellen berechenet, die immer ein volles signal haben
tomwalters@0: time_step=(tstop-tstart-stime)/nr_t;
tomwalters@0: 
tomwalters@0: count=0;
tomwalters@0: w = hann(nr_fft); % for later: the window
tomwalters@0: % for t=56*time_step:time_step:tstop-stime-time_step
tomwalters@0: 
tomwalters@0: sig=setstarttime(sig,0);
tomwalters@0: for t=tstart:time_step:tstop-stime-time_step
tomwalters@0:     count=count+1;
tomwalters@0:     s1=t;
tomwalters@0:     s2=t+stime;
tomwalters@0:     if s2<tstop
tomwalters@0:         s=getpart(sig,t,t+stime);    % von diesem Teil soll die FFT gemacht werden
tomwalters@0:     else % hier sollte er nicht hinkommen
tomwalters@0:         s=getpart(sig,t,tstop);
tomwalters@0:         s=expand(s,stime,0);    % fülle den Rest mit Nullen auf, damit das Signal auf jeden Fall lang genug ist
tomwalters@0:     end
tomwalters@0:     s=setstarttime(s,0);
tomwalters@0:         % damit die FFT nr_f Punkte hat, muss das Signal, womit sie gemacht wird, so lang sein:
tomwalters@0: %     s2=changesr(s,sr_neu);
tomwalters@0:     % ein HanningFenster drüberlegen:
tomwalters@0:     s=s*w;
tomwalters@0:     ps=powerspectrum(s,nr_f*2);
tomwalters@0:     % the resultig powesprect is one point too large, because of the zerovalue->strip it
tomwalters@0:     ps=strippowerspectrum(ps); % cut the leading zero
tomwalters@0: %     figure(1)
tomwalters@0: %     plot(s,'.-');
tomwalters@0: %     figure(2)
tomwalters@0: %     plot(ps,'.-');
tomwalters@0:     % normiere das Spektrum auf die maximale Amplitude des Signals
tomwalters@0:     ms=max(s);
tomwalters@0:     ps=ps*ms;
tomwalters@0:     f=setcolumn(f,count,ps);
tomwalters@0: end