view aim-mat/tools/spectrogram.m @ 0:74dedb26614d

Initial checkin of AIM-MAT version 1.5 (6.4.2011).

% 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