tomwalters@0: % method of class @signal
tomwalters@0: % 
tomwalters@0: %   INPUT VALUES:
tomwalters@0: %  
tomwalters@0: %   RETURN VALUE:
tomwalters@0: %
tomwalters@0: % 
bleeck@3: % 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 fsig=powerspectrum(sig,nr_fft)
tomwalters@0: % usage: sig=powerspectrum(a)
tomwalters@0: % returns a fsignal class containing the powerspectrum of the signal class a
tomwalters@0: % a should have a lenght of 2^n the result however, is one point longer!
tomwalters@0: % in the first bin is additionally the zero value
tomwalters@0: % the result is normed, so that the highest value is 0dB
tomwalters@0: % berechnet das Power Spektrum aus dem Signal a.
tomwalters@0: % die Phasen werden weggeworfen
tomwalters@0: 
tomwalters@0: if nargin < 2
tomwalters@0:     nr_fft=1024;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: Fs = getsr(sig);
tomwalters@0: vals=getvalues(sig);
tomwalters@0: y=fft(vals,nr_fft);
tomwalters@0: py=y.*conj(y);
tomwalters@0: % 
tomwalters@0: % % Verlängerung des Signals um einen Punkt, damit der Nullanteil noch dabei ist
tomwalters@0: nr=round(size(py,1)/2+1); % eines mehr (der Nullanteil)
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % use the powerspectrum from the toolbox
tomwalters@0: % [ppy,w]=periodogram(vals,[],'onesided',nr_fft,Fs); % only real values
tomwalters@0: 
tomwalters@0: % otherwise calculate it yourself:
tomwalters@0: ppy=py(1:nr)/nr*2;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: %normierung auf Energie
tomwalters@0: %s=sum(abs(ppy));
tomwalters@0: %energy=s*s;
tomwalters@0: energy = sum(ppy.^2);
tomwalters@0: if energy==0
tomwalters@0:     ppy=1;
tomwalters@0: else
tomwalters@0:     ppy=ppy/energy; 
tomwalters@0:     ppy=20*log(ppy);
tomwalters@0:     ppy=ppy-max(ppy);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: fsig=fsignal(ppy,nr); % Signal mit der richtigen Samplerate
tomwalters@0: fsig=setdf(fsig,Fs/(nr-1)/2); % kleinester Frequenzabstand
tomwalters@0: 
tomwalters@0: fsig=setmaxfre(fsig,Fs/2);
tomwalters@0: 
tomwalters@0: % sig=setsr(sig,nr); %sr bedeutet für fsignals was anderes, nämlich die Zahl der Punkte
tomwalters@0: fsig=setname(fsig,sprintf('Power Spectrum of Signal \n%s',sig.name));
tomwalters@0: fsig=setunit_x(fsig,'Frequency (Hz)');
tomwalters@0: fsig=setunit_y(fsig,'Power Spectral Density (dB/Hz)');