Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/fftfilt_core.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 % FFTFILT_CORE Simple core method which computes the fft filter.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: Simple core method which computes the fft filter.
5 %
6 % CALL: ao = fftfilt_core(ao, filt, Npad)
7 %
8 % INPUTS: ao: Single input analysis object
9 % Npad: Number of bins for zero padding
10 % filt: The filter to apply to the data
11 % smodel - a model to filter with.
12 % mfir - an FIR filter
13 % miir - an IIR filter
14 % tf - an ltpda_tf object. Including:
15 % - pzmodel
16 % - rational
17 % - parfrac
18 %
19 % VERSION: $Id: fftfilt_core.m,v 1.17 2011/05/30 10:42:32 mauro Exp $
20 %
21 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22
23 function bs = fftfilt_core(varargin)
24
25 bs = varargin{1};
26 filt = varargin{2};
27 Npad = varargin{3};
28
29 if nargin == 4
30 inCondsMdl = varargin{4};
31 else
32 inCondsMdl = [];
33 end
34
35 [m, n] = size(bs.data.y);
36 % FFT time-series data
37 fs = bs.data.fs;
38 % zero padding data before fft
39 if isempty(Npad)
40 Npad = length(bs.data.y) - 1;
41 end
42 if n == 1
43 tdat = ao(tsdata([bs.data.y;zeros(Npad,1)], fs));
44 else
45 tdat = ao(tsdata([bs.data.y zeros(1,Npad)], fs));
46 end
47 % get onesided fft
48 ft = fft_core(tdat, 'one');
49
50 switch class(filt)
51 case 'smodel'
52 % Evaluate model at the given frequencies
53
54 amdly = filt.setXvals(ft.data.x).double;
55 amdly = reshape(amdly, size(ft.data.y));
56
57 amdl = ao(fsdata(ft.data.x, amdly, fs));
58
59 % set units
60 bs.setYunits(simplify(bs.data.yunits .* filt.yunits));
61
62 case {'miir', 'mfir', 'pzmodel', 'parfrac', 'rational'}
63
64 % Check if the frequency of the filter is the same as the frequency
65 % of the AO.
66 if isa (filt, 'ltpda_filter') && fs ~= filt.fs
67 error('### Please use a filter with the same frequency as the AO [%dHz]', fs);
68 end
69
70 % get filter response on given frequencies
71 amdl = resp(filt, plist('f', ft.data.x));
72
73 % set units
74 bs.setYunits(simplify(bs.data.yunits .* filt.ounits ./ filt.iunits));
75
76 case 'ao'
77
78 % check if filter and data have the same shape
79 if size(ft.data.y)~=size(filt.data.y)
80 % reshape
81 amdl = copy(filt,1);
82 amdl.setX(ft.data.x);
83 amdl.setY(reshape(filt.data.y,size(ft.data.x)));
84 amdl.setName(filt.name);
85 else
86 amdl = copy(filt,1);
87 amdl.setName(filt.name);
88 end
89
90 % set units
91 bs.setYunits(simplify(bs.data.yunits .* amdl.data.yunits));
92
93 case 'collection'
94
95 % run over collection elements
96 amdl = ao(fsdata(ft.data.x, ones(size(ft.data.x)), fs));
97 for ii=1:numel(filt.objs)
98 switch class(filt.objs{ii})
99 case 'smodel'
100 % Evaluate model at the given frequencies
101 amdly = filt.objs{ii}.setXvals(ft.data.x).double;
102 amdly = reshape(amdly, size(ft.data.y));
103 amdl_temp = ao(fsdata(ft.data.x, amdly, fs));
104 % set units
105 bs.setYunits(simplify(bs.data.yunits .* filt.objs{ii}.yunits));
106
107 case {'miir'}
108 % get filter response on given frequencies
109 amdly = utils.math.mtxiirresp(filt.objs{ii},ft.data.x,fs,[]);
110 amdl_temp = ao(fsdata(ft.data.x, amdly, fs));
111 % set units
112 bs.setYunits(simplify(bs.data.yunits .* filt.objs{ii}.ounits ./ filt.objs{ii}.iunits));
113
114 case 'ao'
115 % check if filter and data have the same shape
116 if size(ft.data.y)~=size(filt.objs{ii}.data.y)
117 % reshape
118 amdl_temp = copy(filt.objs{ii},1);
119 amdl_temp.setX(ft.data.x);
120 amdl_temp.setY(reshape(filt.objs{ii}.data.y,size(ft.data.x)));
121 amdl_temp.setName(filt.objs{ii}.name);
122 else
123 amdl_temp = copy(filt.objs{ii},1);
124 amdl_temp.setName(filt.objs{ii}.name);
125 end
126 % set units
127 bs.setYunits(simplify(bs.data.yunits .* amdl_temp.data.yunits));
128
129 case 'filterbank'
130 % get filter response on given frequencies
131 amdly = utils.math.mtxiirresp(filt.objs{ii}.filters,ft.data.x,fs,filt.objs{ii}.type);
132 amdl_temp = ao(fsdata(ft.data.x, amdly, fs));
133 % handle units
134 switch lower(filt.objs{ii}.type)
135 case 'parallel'
136 % set units of the output object
137 bs.setYunits(simplify(bs.data.yunits .* filt.objs{ii}.filters(1).ounits ./ filt.objs{ii}.filters(1).iunits));
138 case 'series'
139 % get units from the series
140 sunits = filt.objs{ii}.filters(1).ounits ./ filt.objs{ii}.filters(1).iunits;
141 for jj = 2:numel(filt.objs{ii}.filters)
142 sunits = sunits.*filt.objs{ii}.filters(jj).ounits ./ filt.objs{ii}.filters(jj).iunits;
143 end
144 % set units of the output object
145 bs.setYunits(simplify(bs.data.yunits .* sunits));
146 end
147 end
148 % update response
149 amdl = amdl .* amdl_temp;
150 end
151
152 otherwise
153
154 error('### Unknown filter mode.');
155
156 end
157
158 % Add initial conditions
159 if ~isempty(inCondsMdl) && ~isempty(inCondsMdl.expr.s)
160 inCondsMdl.setXvals(amdl.x);
161 inCondsMdl.setXunits(amdl.xunits);
162 inCondsMdl.setYunits(amdl.yunits*ft.yunits);
163 inCondsEval = inCondsMdl.eval;
164 % Multiply by model and take inverse FFT
165 y = ifft_core(ft.*amdl+inCondsEval, 'symmetric');
166 else
167 % Multiply by model and take inverse FFT
168 y = ifft_core(ft.*amdl, 'symmetric');
169 end
170
171 % split and reshape the data
172 if m == 1
173 bs.setY(y.data.getY(1:n));
174 else
175 bs.setY(y.data.getY(1:m));
176 end
177
178 % clear errors
179 bs.clearErrors;
180
181 end