Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/blwhitenoise.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % BLWHITENOISE return a band limited gaussian distributed white noise
2 %
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 %
5 % DESCRIPTION: BLWHITENOISE return a band limited Gaussian distributed
6 % white noise
7 %
8 % CALL: wn = blwhitenoise()
9 %
10 % INPUTS: npts - number of data points in the series. Note, number of
11 % seconds is npts*fs
12 % fs - sampling frequency
13 % fl - lower bandwidth frequency
14 % fh - higher bandwidth frequency
15 %
16 % OUTPUTS: wn - band limited gaussian white noise
17 %
18 % REFERENCES:
19 %
20 % [1] Kafadar, K., Gaussian white-noise generation for digital signal
21 % synthesis. IEEE Transactions on Instrumentation and Measurement IM-35(4),
22 % 492-495, 1986.
23 %
24 % VERSION: $Id: rand.m,v 1.2 2008/10/24 06:19:23 hewitson Exp $
25 %
26 % HISTORY: 04-01-2010 Luigi Ferraioli
27 % Creation
28 %
29 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
30
31 function Xt = blwhitenoise(npts,fs,fl,fh)
32
33 % generate linear spaced frequency vactor for fft
34 f = (0:floor(npts/2)).*fs./npts;
35 f = f';
36 nf = length(f);
37
38 % amplitude vector
39 amp = zeros(nf,1);
40
41 % phase vactor
42 phs = zeros(nf,1);
43
44 % select frequency range
45 if fh>fs/2
46 fh = fs/2;
47 end
48 if fl<0
49 fl = 0;
50 end
51
52 idxl = f>=fl;
53 idxh = f<=fh;
54 idx = idxl&idxh;
55 n1s = sum(idx);
56
57 phs(idx,1) = -pi + (2*pi).*rand(n1s,1);
58 phs(1,1) = 0;
59 phs(end,1) = 0;
60 amp(idx,1) = sqrt(npts);
61 amp(1,1) = 0;
62 amp(end,1) = 0;
63
64 % get final random phase and constant amplitude
65 phs = [phs;-1.*flipud(phs(2:end-1))];
66 amp = [amp;flipud(amp(2:end-1))];
67
68 % fft signal
69 Xf = amp.*(cos(phs)+1i.*sin(phs));
70
71 % time domain signal
72 Xt = ifft(Xf);
73
74 end
75 % END
76