=========================================================================
Find the frequency terms in the data set with the use of filters and FFT
or a "search" algorithm (slower but more accurate). The function returns
the number of oscillations per unit time where the unit time is define by
DT. Eg If DT is in seconds, then freq is the number of oscillations per
second. If DT is the number of turns, then freq is the number of
oscillations per turn.
[freq(s) amplitude(s) eigenvector(s) time_vec_used] = ...
FREQSEARCH(DATA, [FILTER, METHOD, DT, ORDER, RANGE, TOLERANCE, windowfraction])
DATA : input data, can be complex.
FILTER : 'hanning','none' (default)
METHOD : 'fft' (default),'search','spectrum'
DT : timestep between each data point. (default: 1)
The options below are only applicable to the 'search' method.
ORDER (vector) : number of frequency terms to extract. Ordered by relative
strength. So (default: [1])
RANGE : frequency range over which to scan. (default: [0 Inf])
TOLERANCE : Search until freq_(n) - freq_(n-1) < TOLERANCE (default:
1e-10)
windowfraction : How wide should the subsequent search range should be.
(default: 0.03)
Examples:
>> [f a] = freqsearch(data,'hanning','search',1,[1 2 3])
24/01/2006
Eugene
v2.0 - fft and "search" method combines to increase the speed in which
one can analyse the freqency components.
- Use of filters, modified hanning type filter only at the moment.
- Multiple orders for comparisons of multiple resonant frequencies.
The frequency is ordered by amplitude/strength. The first
frequency being the dominant one followed by the second strongest.
19/08/2010 Eugene: added 'spectrum' option to return a spectrogram. When
in this mode, frequency and amplitude will be vectors
that represent the spectrogram. Range will need to be
specified and has to be a vector at the frequencies of
interest.0001 function varargout = freqsearch3(data, varargin) 0002 % ========================================================================= 0003 % Find the frequency terms in the data set with the use of filters and FFT 0004 % or a "search" algorithm (slower but more accurate). The function returns 0005 % the number of oscillations per unit time where the unit time is define by 0006 % DT. Eg If DT is in seconds, then freq is the number of oscillations per 0007 % second. If DT is the number of turns, then freq is the number of 0008 % oscillations per turn. 0009 % 0010 % [freq(s) amplitude(s) eigenvector(s) time_vec_used] = ... 0011 % FREQSEARCH(DATA, [FILTER, METHOD, DT, ORDER, RANGE, TOLERANCE, windowfraction]) 0012 % 0013 % DATA : input data, can be complex. 0014 % FILTER : 'hanning','none' (default) 0015 % METHOD : 'fft' (default),'search','spectrum' 0016 % DT : timestep between each data point. (default: 1) 0017 % 0018 % The options below are only applicable to the 'search' method. 0019 % 0020 % ORDER (vector) : number of frequency terms to extract. Ordered by relative 0021 % strength. So (default: [1]) 0022 % RANGE : frequency range over which to scan. (default: [0 Inf]) 0023 % TOLERANCE : Search until freq_(n) - freq_(n-1) < TOLERANCE (default: 0024 % 1e-10) 0025 % windowfraction : How wide should the subsequent search range should be. 0026 % (default: 0.03) 0027 % 0028 % Examples: 0029 % >> [f a] = freqsearch(data,'hanning','search',1,[1 2 3]) 0030 % 0031 % 24/01/2006 0032 % Eugene 0033 % v2.0 - fft and "search" method combines to increase the speed in which 0034 % one can analyse the freqency components. 0035 % - Use of filters, modified hanning type filter only at the moment. 0036 % - Multiple orders for comparisons of multiple resonant frequencies. 0037 % The frequency is ordered by amplitude/strength. The first 0038 % frequency being the dominant one followed by the second strongest. 0039 % 19/08/2010 Eugene: added 'spectrum' option to return a spectrogram. When 0040 % in this mode, frequency and amplitude will be vectors 0041 % that represent the spectrogram. Range will need to be 0042 % specified and has to be a vector at the frequencies of 0043 % interest. 0044 0045 DEBUG = false; 0046 0047 %====================================================================== 0048 % Parse input 0049 % Set some defaults (mainly for the search method) 0050 % Min and Max number of iterations for search method 0051 MAXIT = 23; 0052 0053 % What filter to use 0054 Nparam = 1; 0055 if nargin >= Nparam + 1 & ischar(varargin{Nparam}) 0056 switch lower(varargin{Nparam}) 0057 case 'hanning' 0058 filter = 'hanning'; 0059 case 'none' 0060 filter = 'none'; 0061 otherwise 0062 error(sprintf('Unknown filter option %s',varargin{Nparam})); 0063 end 0064 else 0065 filter = 'none'; 0066 end 0067 0068 % What method to use 0069 Nparam = Nparam + 1; 0070 if nargin >= Nparam + 1 & ischar(varargin{Nparam}) 0071 method = lower(varargin{Nparam}); 0072 else 0073 method = 'fft'; 0074 end 0075 0076 % Time step between each sample 0077 Nparam = Nparam + 1; 0078 if nargin >= Nparam + 1 0079 dt = varargin{Nparam}; 0080 else 0081 dt = 1; 0082 end 0083 0084 % Number of terms to extract 0085 Nparam = Nparam + 1; 0086 if nargin >= Nparam + 1 0087 order = varargin{Nparam}; 0088 else 0089 order = 1; 0090 end 0091 0092 % Tune range 0093 Nparam = Nparam + 1; 0094 if nargin >= Nparam + 1 0095 range = varargin{Nparam}; 0096 else 0097 range = [0 0.5]; 0098 end 0099 0100 Nparam = Nparam + 1; 0101 if nargin >= Nparam + 1 0102 tolerance = varargin{Nparam}; 0103 else 0104 tolerance = 1e-10; 0105 end 0106 0107 % Determines how much to zoom in when narrowing the search. 0108 % Depending on the tolerance, the optimal value for the windowfraction 0109 % changes. However 4% seems good enough. 0110 Nparam = Nparam + 1; 0111 if nargin >= Nparam + 1 0112 windowfraction = varargin{Nparam}; 0113 else 0114 windowfraction = 0.04; 0115 end 0116 0117 if DEBUG 0118 fprintf('Options selected: filter(%s) method(%s) dt(%11.3e)\n',... 0119 filter, method, dt); 0120 fprintf(' order(%d) range(%f %f) tolerance(%11.3e) windowfraction(%f)\n',... 0121 order(end), range(1), range(2), tolerance, windowfraction); 0122 end 0123 % Finshed parsing input 0124 %====================================================================== 0125 0126 % Define variables 0127 % Define the time or running parameter against which to calculate the 0128 % freqrange. 0129 neval = length(data); 0130 T2 = dt*(neval-1)/2; %-T/2 --> T/2 0131 t = [-T2:dt:T2]'; 0132 0133 eigenvec = zeros(neval,max(order)); 0134 orthvec = zeros(neval,max(order)); 0135 orthvec_ = zeros(neval,max(order)); 0136 a = zeros(1,max(order)); 0137 nu = zeros(1,max(order)); 0138 0139 % Ensure that data and t are column vectors; 0140 data = reshape(data,neval,1); 0141 datareal = isreal(data); 0142 %====================================================================== 0143 0144 % Remove any DC component in the signal 0145 % data = data - 0.5/T2*local_midpointsum(data); 0146 0147 % What filter to apply to the data 0148 if DEBUG; disp('Calculating filter'); end; 0149 usefilter = 0; 0150 switch filter 0151 case 'hanning' 0152 % Window function that increases hight of the fundamental peak to make it 0153 % easier to pickout. 0154 p = 1; % cosine window order 0155 kai = complex( 2^p*(factorial(p))^2*(1+cos(pi*(t/T2))).^p/factorial(2*p) ); 0156 usefilter = 1; 0157 end 0158 % Finished applying filter 0159 %====================================================================== 0160 0161 % What method to use 0162 if DEBUG; disp('Starting calculation'); end; 0163 switch method 0164 case 'fft' 0165 [nu a] = local_calculate_with_fft(data,dt,range); 0166 order = 1; 0167 0168 case 'search' 0169 0170 for k=1:max(order) 0171 % Start the frequency search using a two step approach, first 0172 % use the FFT to get a coarse measurement of the frequency 0173 % followed by the correlation analysis to get a more accurate 0174 % measure of the dominant frequency component. 0175 0176 % FFT 0177 if usefilter 0178 prelim_freq = local_calculate_with_fft(data.*kai,dt,range); 0179 else 0180 prelim_freq = local_calculate_with_fft(data,dt,range); 0181 end 0182 0183 % Will scan this range of frequencies. 0184 freqrange = local_find_new_range(prelim_freq,range(2),range(1),windowfraction); 0185 0186 % Some initial variables. Start with some guess at the 0187 % frequency, mainly for the first difference calculation. 0188 % This is the power spectrum/frequency scan. 0189 psi = zeros(1,length(freqrange)); 0190 freq = zeros(1,MAXIT); 0191 0192 omega_prev = median(freqrange); 0193 difference = 1; 0194 0195 for j=1:MAXIT 0196 % Do the integral that calculates the average <f(t), e^i*freqrange*t>. Not 0197 % including multiplication of some factors like dt, since we only 0198 % need to find where psi is a maximum and extract the corresponding 0199 % freqrange. Vectorising this loop does not help, 0200 % evaluated already. 0201 if usefilter 0202 psi = local_psi_integral(data.*kai,t,freqrange); 0203 else 0204 psi = local_psi_integral(data,t,freqrange); 0205 end 0206 0207 if j >= 1 && j <=1 && DEBUG 0208 figure; plot(freqrange,abs(psi)); 0209 xlabel('freq / Frequency'); ylabel('Arb. Units'); 0210 end 0211 0212 % Calculate the value of freqrange for the maximum psi. 0213 [maxpsi maxind] = max(psi(:)); 0214 freq(j) = freqrange(maxind); 0215 0216 difference = abs(freq(j) - omega_prev); 0217 if difference < tolerance 0218 if DEBUG; fprintf('Difference less than specified tolerance. j=%d\n',j); end; 0219 break; 0220 else 0221 omega_prev = freq(j); 0222 end 0223 0224 % Find new range to seach, zoom in. 0225 freqrange = local_find_new_range(freq(j),freqrange(end),freqrange(1),windowfraction); 0226 0227 psi = zeros(size(freqrange)); 0228 end 0229 if DEBUG; fprintf('FREQ = %20.10e\n',freq(1:j)); end; 0230 % Orthogonal projection to determine the coeffients. Since 0231 % e^i*2pi*freq*t. 0232 0233 eigenvec(:,k) = exp(complex(0,2*pi*freq(j).*(t))); 0234 % Orthogonalize 0235 % sumprojections = zeros(neval,1); 0236 % for ii=1:k-1 0237 % sumprojections = sumprojections + dot(eigenvec(:,k),orthvec(:,ii))/dot(orthvec(:,ii),orthvec(:,ii))*orthvec(:,ii); 0238 % end 0239 % orthvec(:,k) = eigenvec(:,k) - sumprojections; 0240 0241 a(k) = ((0.5/T2)*local_midpointsum(data.*conj(eigenvec(:,k))))*dt; 0242 % a(k) = (0.5/T2)*maxpsi; 0243 nu(k) = freq(j); 0244 0245 % Subtract the component from 'f' function. 0246 data = data - a(k)*eigenvec(:,k); 0247 end 0248 case 'spectrum' 0249 % Return the power spectrum 0250 if usefilter 0251 psi = local_psi_integral(data.*kai,t,range); 0252 else 0253 psi = local_psi_integral(data,t,range); 0254 end 0255 if datareal 0256 nu = psi*dt/T2; 0257 else 0258 nu = 0.5*psi*dt/T2; 0259 end 0260 order = 1:length(psi); 0261 0262 otherwise 0263 error(sprintf('Unknown method option %s',varargin{Nparam})); 0264 end 0265 0266 varargout{1} = nu(order); 0267 if nargout > 1 0268 if datareal 0269 % With only real data the returned amplitudes should also be real. 0270 % And the factor 2 is needed here but not quite sure why just yet. 0271 varargout{2} = 2*abs(a(order)); 0272 else 0273 varargout{2} = a(order); 0274 end 0275 end 0276 if nargout > 2 0277 varargout{3} = eigenvec(:,order); 0278 end 0279 if nargout > 3 0280 varargout{4} = t; 0281 end 0282 0283 % temp = freq(find(freq ~= 0)); 0284 % varargout{1} = temp(end)/(2*pi); 0285 % DEBUG 0286 % fprintf('%i %17.15g\n',j, difference); 0287 0288 0289 function fctnsum = local_midpointsum(fctn) 0290 % Vectorise the midpoint method of integrating a numerical function. 0291 % f_n = [fctn 0]; 0292 % f_nplus1 = [0 fctn]; % shift all the numbers one "space" to the right. 0293 % midpoints = 0.5*(f_n + f_nplus1); 0294 0295 midpoints = 0.5.*(fctn(1:end-1) + fctn(2:end)); 0296 fctnsum = sum(midpoints); 0297 0298 0299 function psi = local_psi_integral(data,t,freqrange) 0300 0301 0302 midpoints = zeros(1,length(t)-1); 0303 omegarange = -freqrange*2*pi; 0304 fctn = zeros(1,length(t)); 0305 psi = zeros(1,length(omegarange)); 0306 0307 for k=1:length(freqrange) 0308 fctn = data.*exp(complex(0,omegarange(k)).*t); 0309 % fctn = data.*complex(cos(freqrange(k).*t),-sin(freqrange(k).*t)); 0310 midpoints = 0.5.*(fctn(1:end-1) + fctn(2:end)); 0311 psi(k) = sum(midpoints); 0312 end 0313 0314 0315 function [freq varargout] = local_calculate_with_fft(data,dt,range) 0316 % Find peak with FFT within "range" of frequencies. 0317 % Auto calculate number of points to calculate fft. Use maximum 0318 nn = [4:15]; 0319 ind = max(find(2.^nn - length(data) < 0)); 0320 Nfft = 2^nn(ind); 0321 0322 % Calculate FFT and the power spectrum 0323 yy = fft(data,Nfft); 0324 Pyy = yy.*conj(yy); 0325 0326 % Corresponding frequency range 0327 f = 1/(dt*Nfft).*(0:Nfft/2); 0328 ii = find(f > range(1) & f < range(2)); 0329 0330 % Find peak 0331 [maxval maxind] = max(Pyy(ii)); 0332 freq = f(ii(maxind)); 0333 if nargout > 1 0334 varargout{1} = abs(yy(ii(maxind)))/(Nfft/2); 0335 end 0336 0337 0338 function freqrange = local_find_new_range(centre,upper,lower,windowfraction) 0339 0340 % Find new range to seach, zoom in. 0341 new_width = (upper - lower)*windowfraction; 0342 % minimum frequency separation 0343 min_freq_sep = new_width/500; 0344 0345 if centre == lower 0346 lowerbound = centre - new_width*2/windowfraction; 0347 upperbound = centre; 0348 0349 freqrange = [lowerbound:(upperbound-lowerbound)/100:upperbound]; 0350 elseif centre == upper 0351 lowerbound = centre; 0352 upperbound = centre + new_width*2/windowfraction; 0353 0354 freqrange = [lowerbound:(upperbound-lowerbound)/100:upperbound]; 0355 else 0356 lowerbound = centre - new_width; 0357 upperbound = centre + new_width; 0358 0359 % num = 15; 0360 % scalefactor = (2*new_width - min_freq_sep*(num+1))/(num+2); 0361 % freqrange = lowerbound + cumsum((1 + cos(0:2*pi/num:2*pi))*scalefactor + min_freq_sep); 0362 0363 scalefactor = (2*new_width - min_freq_sep*16)/17; 0364 freqrange = lowerbound + cumsum((1 + cos(0:2*pi/15:2*pi))*scalefactor + min_freq_sep); 0365 end 0366 0367 % freqrange = freqrange*2*pi; 0368 0369