Mercurial > forge

comparison main/signal/pwelch.m @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children 13d7efaa47b8
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
1 ## Copyright (C) 1999-2001 Paul Kienzle
2 ##
3 ## This program is free software; you can redistribute it and/or modify
4 ## it under the terms of the GNU General Public License as published by
5 ## the Free Software Foundation; either version 2 of the License, or
6 ## (at your option) any later version.
7 ##
8 ## This program is distributed in the hope that it will be useful,
9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 ## GNU General Public License for more details.
12 ##
13 ## You should have received a copy of the GNU General Public License
14 ## along with this program; if not, write to the Free Software
15 ## Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16
17 ## usage: [Pxx, w] = pwelch(x,n,Fs,window,overlap,ci,range,units,trend)
18 ## [Pxx, Pci, w] = pwelch(x,n,Fs,window,overlap,ci,range,units,trend)
19 ##
20 ## Estimate power spectrum of a stationary signal. This chops the signal
21 ## into overlapping slices, windows each slice and applies a Fourier
22 ## transform to determine the frequency components at that slice. The
23 ## magnitudes of these slices are then averaged to produce the estimate Pxx.
24 ## The confidence interval around the estimate is returned in Pci.
25 ##
26 ## x: vector of samples
27 ## n: size of fourier transform window, or [] for default=256
28 ## Fs: sample rate, or [] for default=2 Hz
29 ## window: shape of the fourier transform window, or [] for default=hanning(n)
30 ## Note: window length can be specified instead, in which case
31 ## window=hanning(length)
32 ## overlap: overlap with previous window, or [] for default=length(window)/2
33 ## ci: confidence interval, or [] for default=0.95
34 ## ci must be between 0 and 1; if ci is not passed, or if it is
35 ## passed as 0, then no confidence intervals will be computed.
36 ## range: 'whole', or [] for default='half'
37 ## show all frequencies, or just half of the frequencies
38 ## units: 'squared', or [] for default='db'
39 ## show results as magnitude squared or as log magnitude squared
40 ## trend: 'mean', 'linear', or [] for default='none'
41 ## remove trends from the data slices before computing spectral estimates
42 ##
43 ## Example
44 ## [b,a] = cheby1(4,3,[0.2, 0.4]); ## define noise colour
45 ## pwelch(filter(b,a,randn(2^12,1))); ## estimate noise colour
46
47 ## 2001-04-02 Paul Kienzle
48 ## * return nfft/2+1 elements rather than nfft/2 for even nfft.
49 ## * use more accurate (and faster) computation of confidence intervals
50
51 ## TODO: Should be extended to accept a vector of frequencies at which to
52 ## TODO: evaluate the fourier transform (via filterbank or chirp
53 ## TODO: z-transform).
54 ## TODO: What should happen with the final window when it isn't full?
55 ## TODO: currently I dump it, but I should probably zero pad and add
56 ## TODO: it in.
57 ## TODO: Consider returning the whole gamit of Pxx, Pyy, Pxy, Cxy, Txy
58 ## TODO: as well as confidence intervals for each; if users tend
59 ## TODO: only to use one of these don't bother, but if they ever need
60 ## TODO: more than one, then it's free. Alternatively, break out the
61 ## TODO: compute engine into a function that the user can call directly.
62 ## TODO: Check if Cxy, Txy are computed frame-by-frame or on the average
63 ## TODO: of the frames. SpcTools and I do it on the average,
64 ## TODO: wdkirby@ix.netcom.com (1998-04-29 octave-sources) computes
65 ## TODO: them frame-by-frame.
66 function [...] = pwelch(x, ...)
67 ## sort out parameters
68 if nargin < 1,
69 usage("[Pxx, w] = pwelch(x,nfft,Fs,window,overlap,pc,range,units,trend)");
70 endif
71 va_start();
72
73 ## Determine if we are called as pwelch, csd, cohere or tfe
74 if isstr(x)
75 calledby = x;
76 else
77 calledby = "pwelch";
78 endif
79 if !isstr(x)
80 ftype = 1;
81 elseif strcmp(x, 'csd')
82 ftype = 2;
83 elseif strcmp(x, 'cohere')
84 ftype = 3;
85 elseif strcmp(x, 'tfe')
86 ftype = 4;
87 endif
88
89 ## Sort out x and y vectors
90 if ftype!=1
91 x=va_arg(); y=va_arg();
92 first = 4;
93 else
94 y=[];
95 first = 2;
96 endif
97 if (columns(x) != 1 && rows(x) != 1) || ...
98 (!isempty(y) && columns(y) != 1 && rows(y) != 1)
99 error ([calledby, " data must be a vector"]);
100 end
101 if columns(x) != 1, x = x'; end
102 if columns(y) != 1, y = y'; end
103 if !isempty(y) && rows(x)!=rows(y)
104 error ([calledby, " x and y vectors must be the same length"]);
105 endif
106
107 ## interpret remaining arguments
108 trend=nfft=Fs=window=overlap=whole=use_dB=[];
109 ci=-1; ## need to do stupid things with ci
110 pos=0; ## no positional parameters yet interpreted.
111 for i=first:nargin
112 arg = va_arg();
113 if isstr(arg),
114 arg=tolower(arg);
115 if strcmp(arg, 'squared')
116 use_dB = 0;
117 elseif strcmp(arg, 'db')
118 use_dB = 1;
119 elseif strcmp(arg, 'whole')
120 whole = 1;
121 elseif strcmp(arg, 'half')
122 whole = 0;
123 elseif strcmp(arg, 'none')
124 trend = -1;
125 elseif strcmp(arg, 'mean')
126 trend = 0;
127 elseif strcmp(arg, 'linear')
128 trend = 1;
129 else
130 error([calledby, " doesn't understand '", arg, "'"]);
131 endif
132 elseif pos == 0
133 nfft = arg;
134 pos++;
135 elseif pos == 1
136 Fs = arg;
137 pos++;
138 elseif pos == 2
139 window = arg;
140 pos++;
141 elseif pos == 3
142 overlap = arg;
143 pos++;
144 elseif pos == 4
145 ci = arg;
146 pos++;
147 else
148 usage(usagestr);
149 endif
150 endfor
151
152 ## Fill in defaults for arguments that aren't specified
153 if isempty(nfft), nfft = min(256, length(x)); endif
154 if isempty(Fs), Fs = 2; endif
155 if isempty(window), window = hanning(nfft); endif
156 if isempty(overlap), overlap = length(window)/2; endif
157 if isempty(whole), whole = !isreal(x)||(!isempty(y)&&!isreal(y)); endif
158 if isempty(trend), trend=-1; endif
159 if isempty(use_dB),
160 ## don't default to db for cohere, or for returned values
161 use_dB = (ftype!=3 && nargout == 0);
162 endif
163 if isempty(ci), ci=0.95; endif # if ci was not passed in, it would be 0
164
165 ## sort out default confidence intervals
166 if (ci < 0) # ci was not passed in
167 if nargout > 2
168 ci = 0.95;
169 else
170 ci = 0.0;
171 endif
172 endif
173 if (ftype > 2 && ci > 0)
174 error([ calledby, " can't compute confidence intervals" ]);
175 elseif (ci < 0 || ci > 1)
176 error([ calledby, " confidence interval must be between 0 and 1" ]);
177 endif
178
179 ## if only the window length is given, generate hanning window
180 if length(window) == 1, window = hanning(window); endif
181 if rows(window)==1, window = window.'; endif
182
183 ## Normalize the window
184 window = window / norm(window);
185
186 ## compute window offsets
187 win_size = length(window);
188 if (win_size > nfft)
189 nfft = win_size;
190 warning (sprintf("%s fft size adjusted to %d", calledby, n));
191 end
192 step = win_size - overlap;
193
194 ## Determine which correlations to compute
195 Pxx = Pyy = Pxy = [];
196 if ftype!=2, Pxx = zeros(nfft,1); endif # Not needed for csd
197 if ftype==3, Pyy = zeros(nfft,1); endif # Only needed for cohere
198 if ftype!=1, Pxy = zeros(nfft,1); endif # Not needed for psd
199
200 ## Average the slices
201 offset = 1:step:length(x)-win_size+1;
202 N = length(offset);
203 for i=1:N
204 a = x(offset(i):offset(i)+win_size-1);
205 if trend>=0, a=detrend(a,trend); endif
206 a = fft(postpad(a.*window, nfft));
207 if !isempty(Pxx), Pxx = Pxx + a.*conj(a); endif
208 if !isempty(Pxy)
209 b = y(offset(i):offset(i)+win_size-1);
210 if trend>=0, b=detrend(b,trend); endif
211 b = fft(postpad(b.*window, nfft));
212 Pxy = Pxy + a .*conj(b);
213 if !isempty(Pyy), Pyy = Pyy + b.*conj(b); endif
214 endif
215 endfor
216 if (ftype <= 2)
217 ## the factors of N cancel when computing cohere and tfe
218 if !isempty(Pxx), Pxx = Pxx / N; endif
219 if !isempty(Pxy), Pxy = Pxy / N; endif
220 if !isempty(Pyy), Pyy = Pyy / N; endif
221 endif
222
223 ## Compute confidence intervals
224 if ci > 0, Pci = zeros(nfft,1); endif
225 if (ci > 0 && N > 1)
226 if ftype>2
227 error([calledby, ": internal error -- shouldn't compute Pci"]);
228 end
229
230 ## c.i. = mean +/- dev
231 ## dev = z_ci*std/sqrt(n)
232 ## std = sqrt(sumsq(P-mean(P))/(N-1))
233 ## z_ci = normal_inv( 1-(1-ci)/2 ) = normal_inv( (1+ci)/2 );
234 ## normal_inv(x) = sqrt(2) * erfinv(2*x-1)
235 ## => z_ci = sqrt(2)*erfinv(2*(1+ci)/2-1) = sqrt(2)*erfinv(ci)
236 for i=1:N
237 a=x(offset(i):offset(i)+win_size-1);
238 if trend>=0, a=detrend(a,trend); endif
239 a=fft(postpad(a.*window, nfft));
240 if ftype == 1 # psd
241 P = a.*conj(a) - Pxx;
242 Pci = Pci + P.*conj(P);
243 else # csd
244 b=y(offset(i):offset(i)+win_size-1);
245 if trend>=0, b=detrend(b,trend); endif
246 b=fft(postpad(b.*window, nfft));
247 P = a.*conj(b) - Pxy;
248 Pci = Pci + P.*conj(P);
249 endif
250 endfor
251
252 Pci = ( erfinv(ci) * sqrt( 2/N/(N-1) ) ) * sqrt ( Pci );
253 endif
254
255 switch (ftype)
256 case 1, # psd
257 P = Pxx / Fs;
258 if ci > 0, Pci = Pci / Fs; endif
259 case 2, # csd
260 P = Pxy;
261 case 3, # cohere
262 P = Pxy.*conj(Pxy)./Pxx./Pyy;
263 case 4, # tfe
264 P = Pxy./Pxx;
265 endswitch
266
267 ## compute confidence intervals
268 if ci > 0, Pci = [ P - Pci, P + Pci ]; endif
269
270 if use_dB
271 P = 10.0*log10(P);
272 if ci > 0, Pci = 10.0*log10(Pci); endif
273 endif
274
275 ## extract the positive frequency components
276 if whole
277 ret_n = nfft;
278 elseif rem(nfft,2)==1
279 ret_n = (nfft+1)/2;
280 else
281 ret_n = nfft/2 + 1;
282 end
283 P = P(1:ret_n, :);
284 if ci > 0, Pci = Pci(1:ret_n, :); endif
285 f = [0:ret_n-1]*Fs/nfft;
286
287 ## Plot if there is no
288 if nargout==0,
289 unwind_protect
290 if Fs==2
291 xlabel("Frequency (rad/pi)");
292 else
293 xlabel("Frequency (Hz)");
294 endif
295 if ftype==1
296 title ("Welch's Spectral Estimate Pxx/Fs");
297 ytext="Power Spectral Density";
298 elseif ftype==2
299 title ("Cross Spectral Estimate Pxy");
300 ytext="Cross Spectral Density";
301 elseif ftype==3
302 title ("Coherence Function Estimate |Pxy|^2/(PxxPyy)");
303 ytext="Coherence ";
304 else
305 title ("Transfer Function Estimate Pxy/Pxx");
306 ytext="Transfer";
307 endif
308 if use_dB,
309 ylabel(strcat(ytext, " (dB)"));
310 else
311 ylabel(ytext);
312 endif
313 grid("on");
314 if ci>0
315 plot(f, [P, Pci], ";;");
316 else
317 plot(f, P, ";;");
318 endif
319 unwind_protect_cleanup
320 grid("off");
321 title("");
322 xlabel("");
323 ylabel("");
324 end_unwind_protect
325 endif
326
327 if nargout>=1, vr_val(P); endif
328 if nargout>=2 && ci>0, vr_val(Pci); endif
329 if nargout>=2 && ci==0, vr_val(f); endif
330 if nargout>=3 && ci>0, vr_val(f); endif
331
332 endfunction
333
334 %!demo
335 %! Fs=8000;
336 %! [b,a] = cheby1(4,3,2*[500, 1000]/Fs); ## define spectral envelope
337 %! s=0.05*randn(2^11,1); ## define noise
338 %! idx=fix(1:Fs/70:length(s))';
339 %! s(idx)=s(idx)+ones(size(idx)); ## add 70 Hz excitation
340 %! x=filter(b,a,s); ## generate signal
341 %!
342 %! figure(1); subplot(221);
343 %! text(0,0.9,'basic estimate','Units','Normalized');
344 %! pwelch(x',[],Fs); text; # slip in a test for row vs. column vector
345 %! subplot(222);
346 %! text(0,0.9,'nfft=1024 instead of 256','Units','Normalized');
347 %! pwelch(x,1024); text;
348 %! subplot(223);
349 %! text(0,0.9,'boxcar instead of hanning','Units','Normalized');
350 %! pwelch(x,[],[],boxcar(256)); text;
351 %! subplot(224);
352 %! text(0,0.9,'no overlap','Units','Normalized');
353 %! pwelch(x,[],[],[],0); text;
354 %!
355 %! figure(2); subplot(121);
356 %! text(0,0.9,'magnitude units, whole range','Units','Normalized');
357 %! pwelch(x,'whole','squared'); text;
358 %! subplot(122);
359 %! text(0,0.9,'90% confidence intervals','Units','Normalized');
360 %! pwelch(x,[],[],[],[],0.9); text;
361 %! oneplot();
362 %! %----------------------------------------------------------
363 %! % plots should show a chebyshev bandpass filter shape