|
tomwalters@0
|
1 % method of class @signal
|
|
tomwalters@0
|
2 %
|
|
tomwalters@0
|
3 % INPUT VALUES:
|
|
tomwalters@0
|
4 %
|
|
tomwalters@0
|
5 % RETURN VALUE:
|
|
tomwalters@0
|
6 %
|
|
tomwalters@0
|
7 %
|
|
bleeck@3
|
8 % This external file is included as part of the 'aim-mat' distribution package
|
|
bleeck@3
|
9 % (c) 2011, University of Southampton
|
|
bleeck@3
|
10 % Maintained by Stefan Bleeck (bleeck@gmail.com)
|
|
bleeck@3
|
11 % download of current version is on the soundsoftware site:
|
|
bleeck@3
|
12 % http://code.soundsoftware.ac.uk/projects/aimmat
|
|
bleeck@3
|
13 % documentation and everything is on http://www.acousticscale.org
|
|
bleeck@3
|
14
|
|
tomwalters@0
|
15
|
|
tomwalters@0
|
16 function fsig=powerspectrum(sig,nr_fft)
|
|
tomwalters@0
|
17 % usage: sig=powerspectrum(a)
|
|
tomwalters@0
|
18 % returns a fsignal class containing the powerspectrum of the signal class a
|
|
tomwalters@0
|
19 % a should have a lenght of 2^n the result however, is one point longer!
|
|
tomwalters@0
|
20 % in the first bin is additionally the zero value
|
|
tomwalters@0
|
21 % the result is normed, so that the highest value is 0dB
|
|
tomwalters@0
|
22 % berechnet das Power Spektrum aus dem Signal a.
|
|
tomwalters@0
|
23 % die Phasen werden weggeworfen
|
|
tomwalters@0
|
24
|
|
tomwalters@0
|
25 if nargin < 2
|
|
tomwalters@0
|
26 nr_fft=1024;
|
|
tomwalters@0
|
27 end
|
|
tomwalters@0
|
28
|
|
tomwalters@0
|
29
|
|
tomwalters@0
|
30
|
|
tomwalters@0
|
31 Fs = getsr(sig);
|
|
tomwalters@0
|
32 vals=getvalues(sig);
|
|
tomwalters@0
|
33 y=fft(vals,nr_fft);
|
|
tomwalters@0
|
34 py=y.*conj(y);
|
|
tomwalters@0
|
35 %
|
|
tomwalters@0
|
36 % % Verlängerung des Signals um einen Punkt, damit der Nullanteil noch dabei ist
|
|
tomwalters@0
|
37 nr=round(size(py,1)/2+1); % eines mehr (der Nullanteil)
|
|
tomwalters@0
|
38
|
|
tomwalters@0
|
39
|
|
tomwalters@0
|
40 % use the powerspectrum from the toolbox
|
|
tomwalters@0
|
41 % [ppy,w]=periodogram(vals,[],'onesided',nr_fft,Fs); % only real values
|
|
tomwalters@0
|
42
|
|
tomwalters@0
|
43 % otherwise calculate it yourself:
|
|
tomwalters@0
|
44 ppy=py(1:nr)/nr*2;
|
|
tomwalters@0
|
45
|
|
tomwalters@0
|
46
|
|
tomwalters@0
|
47
|
|
tomwalters@0
|
48 %normierung auf Energie
|
|
tomwalters@0
|
49 %s=sum(abs(ppy));
|
|
tomwalters@0
|
50 %energy=s*s;
|
|
tomwalters@0
|
51 energy = sum(ppy.^2);
|
|
tomwalters@0
|
52 if energy==0
|
|
tomwalters@0
|
53 ppy=1;
|
|
tomwalters@0
|
54 else
|
|
tomwalters@0
|
55 ppy=ppy/energy;
|
|
tomwalters@0
|
56 ppy=20*log(ppy);
|
|
tomwalters@0
|
57 ppy=ppy-max(ppy);
|
|
tomwalters@0
|
58 end
|
|
tomwalters@0
|
59
|
|
tomwalters@0
|
60 fsig=fsignal(ppy,nr); % Signal mit der richtigen Samplerate
|
|
tomwalters@0
|
61 fsig=setdf(fsig,Fs/(nr-1)/2); % kleinester Frequenzabstand
|
|
tomwalters@0
|
62
|
|
tomwalters@0
|
63 fsig=setmaxfre(fsig,Fs/2);
|
|
tomwalters@0
|
64
|
|
tomwalters@0
|
65 % sig=setsr(sig,nr); %sr bedeutet für fsignals was anderes, nämlich die Zahl der Punkte
|
|
tomwalters@0
|
66 fsig=setname(fsig,sprintf('Power Spectrum of Signal \n%s',sig.name));
|
|
tomwalters@0
|
67 fsig=setunit_x(fsig,'Frequency (Hz)');
|
|
tomwalters@0
|
68 fsig=setunit_y(fsig,'Power Spectral Density (dB/Hz)');
|
