Mercurial > hg > aimmat
diff aim-mat/tools/spectrogram.m @ 0:74dedb26614d
Initial checkin of AIM-MAT version 1.5 (6.4.2011).
author | tomwalters |
---|---|
date | Fri, 20 May 2011 12:32:31 +0100 |
parents | |
children | 20ada0af3d7d |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/aim-mat/tools/spectrogram.m Fri May 20 12:32:31 2011 +0100 @@ -0,0 +1,115 @@ +% support file for 'aim-mat' +% +% This external file is included as part of the 'aim-mat' distribution package +% http://www.pdn.cam.ac.uk/cnbh/aim2006 +% $Date: 2008-06-10 18:00:16 +0100 (Tue, 10 Jun 2008) $ +% $Revision: 585 $ + +function f=spectrogram(sig,nr_f,maxfre) + +nsig=changesr(sig,maxfre*2); +sr=getsr(nsig); +vals=getvalues(nsig); + +% sig_len=getlength(sig); +% time_step=sig_len/nr_t; +% window=[0:time_step:sig_len]; +% win_len=getnrpoints(nsig)/nr_t; +% window = hann(round(win_len)); +% noverlap=round(win_len/2); +noverlap=nr_f*2-2; +[Bb,freqs,ts]=specgram(vals,nr_f*2,sr*2,[],noverlap); + +B=Bb(2:nr_f+1,:); + +nr_t=size(B,2); +nr_f=size(B,1); + +srf=nr_t/getlength(sig); +f=field(nr_t,nr_f,srf); +f=setvalues(f,abs(B)); +f=setmaxfre(f,maxfre); + +return + + + + +% function f=spectrogram(sig,tstart,tstop,nr_t,nr_f,maxfre) +% alle Zeiten in Sekunden +% macht ein Sektrogramm aus dem Signal sig, das zu nr_t Zeiten gesampelt +% wird und nr_f Frequenzen hat +% das Signal wird aus sig in der Zeit von start bis stop genommen +% das Signal sig wird erst so auf eine Samplerate gebracht, dass +% das Resultat der FFT eine Auflösung bis maxfre hat. + +if nargin<6 + maxfre=GetSR(sig)/2; +end + +% die alte (orginale) Samplerate +sr_old=GetSR(sig); + +%wenn die maxfre anders ist als die, die das Signal mitbringt, dann muss das Signal heruntergesamplet werden: +if sr_old~=maxfre*2 + sr_neu=maxfre*2; + sig=changesr(sig,sr_neu); % änder das ganze Signal zu der neuen SR +else + sr_neu=sr_old; +end + +% Punkte, die die FFT benötigt, damit nr_f punkte hinten raus kommen +nr_fft=nr_f*2; +% und das ist ein Signal von dieser Länge: (auch mit neuer Samplerate) +stime=bin2time(nr_fft,sr_neu); +% Dauer des zu untersuchenden Signals +duration=tstop-tstart; +% wenn die gewünschte Samplingrate nicht hoch genug ist, weil zuviel Punkte erforderlich sind +if duration < stime + error('spectrogram: desired Samplingrate too low: take less points or higher freqeuncy!'); + return; +end + + +% meine Notation ist anders als Matlab: Bei mir ist x die Zeit +% und y die Frequenz +sr_field=nr_t/duration; % die "SampleRate" des Feldes (sehr klein, da nur wenig Punkte) +f=field(nr_t,nr_f,sr_field); +f=setmaxfre(f,maxfre); % Für die Grafik: Wie hoch die Frequenzen geben + +% damit die fft immer von einem richtigen vollen Signal gemacht werden kann, werden die +% FFTs an Stellen berechenet, die immer ein volles signal haben +time_step=(tstop-tstart-stime)/nr_t; + +count=0; +w = hann(nr_fft); % for later: the window +% for t=56*time_step:time_step:tstop-stime-time_step + +sig=setstarttime(sig,0); +for t=tstart:time_step:tstop-stime-time_step + count=count+1; + s1=t; + s2=t+stime; + if s2<tstop + s=getpart(sig,t,t+stime); % von diesem Teil soll die FFT gemacht werden + else % hier sollte er nicht hinkommen + s=getpart(sig,t,tstop); + s=expand(s,stime,0); % fülle den Rest mit Nullen auf, damit das Signal auf jeden Fall lang genug ist + end + s=setstarttime(s,0); + % damit die FFT nr_f Punkte hat, muss das Signal, womit sie gemacht wird, so lang sein: +% s2=changesr(s,sr_neu); + % ein HanningFenster drüberlegen: + s=s*w; + ps=powerspectrum(s,nr_f*2); + % the resultig powesprect is one point too large, because of the zerovalue->strip it + ps=strippowerspectrum(ps); % cut the leading zero +% figure(1) +% plot(s,'.-'); +% figure(2) +% plot(ps,'.-'); + % normiere das Spektrum auf die maximale Amplitude des Signals + ms=max(s); + ps=ps*ms; + f=setcolumn(f,count,ps); +end