FSCREATE Create Fourier Series Function.
[FSfun,Kn] = FSCREATE(...) creates a function handle FSfun containing the
complex exponential Fourier Series created from its input arguments. The
optional output Kn is a row vector containing the complex exponential
Fourier Series coefficients: Kn = [K(-N) ... K(0) ... K(N)]
FSfun = FSCREATE(t,f,N,TYPE) computes the Fourier Series of a signal
tabulated in the real-valued data vectors t and f.
t is a monotonically increasing time vector with the signal period being
T = t(end)-t(1). f(1)=f(end) is required.
N is the number of harmonics. If N is empty or not given, N = 64.
TYPE is a string describing the signal type described by the data t and f.
If TYPE is not given or if TYPE = 'foh' (First Order Hold), the Fourier
Series coefficients are computed by assuming that the input data is
piecewise linear between the data points.
If TYPE = 'zoh' (Zero Order Hold), the Fourier Series coefficients are
computed by assuming that the input data is piecewise constant between
the data points.
The FFT is not used, so no aliasing occurs.
FSfun = FSCREATE(WAVE,T,N,P) returns the Fourier Series of known
waveforms. WAVE is a character string identifying the desired waveform,
T is the period, and N is the number of harmonics computed. If N is empty
or not given, N = 64. P is a vector of optional parameters required for
some waveforms. See below for waveforms and parameters.
FSfun = FSCREATE(Kn,T) creates a Fourier Series function using the Fourier
series vector Kn and period T. Kn is a row vector of the complex
exponential Fourier series coefficients Kn = [K(-N) ... K(0) ... K(N)]
where K(n) is the n-th term, e.g., K(0) is the DC term, and N is the
number of harmonics.
FSfun2 = FSCREATE(FSfun,'Op',P) performs the operation specified by the
string 'Op' on the Fourier Series function handle FSfun, returning a new
Fourier series function handle. P is a parameter needed for some 'Op'.
Operation DESCRIPTION
diff differentiate the Fourier Series.
int integrate the Fourier Series, ignoring the DC or average term
and returning a Fourier Series having zero average.
mirror time mirror (Fourier Series of f(-t)).
smooth apply Blackman window to the Fourier Series (minimize Gibb's)
trim trim negligible real and imaginary elements.
even extract the even time part (real parts).
odd extract the odd time part (imaginary parts).
halfwave extract the halfwave part, (odd harmonics).
nodc remove DC or average value term (K(0) = 0).
dc P = desired DC or average value (K(0) = P).
add P = amount to add to DC or average value (K(0)=K(0)+P).
delay P = normalized delay to apply to FSfun, P = 1 is one period.
delay P = 'odd' or P = 'even' apply delay to achieve given symmetry.
scale P = amplitude scale factor to apply to FSfun (P*Kn)
period P = desired period to associate with Fourier Series.
resize P = number of harmonics desired in result. If P<N terms are
deleted; if P>N zeros are padded onto Kn.
tag P = tag character string to store with function handle
FSfun = FSCREATE(FSfunA,'Op',FSfunB) performs the mathematical operation
specified by the string 'Op' on the underlying time signals in the input
Fourier Series functions. 'Op' is one of the following '+' (addition),
'-' (subtraction), '*' (multiplication). The period of FSfunA and FSfunB
are assumed to be equal.
FSCREATE(...,Tag) where character string Tag is the last input argument,
stores the string Tag in the function handle for possible retrieval later.
Tag can be any character string.
FSCREATE('WAVE',T,N,P) Waveforms
WAVE DESCRIPTION
'square' square wave, odd symmetry, zero mean,
'saw' sawtooth, positive slope, odd symmetry, zero mean,
P = fractional fall time, 0 < P < 0.5. If not given, P = 0.
'rsaw' reversed sawtooth, same as 'sawtooth' but negative slope
P = fractional fall time, 0 < P < 0.5. If not given, P = 0.
'triangle' triangle, even symmetry, zero mean.
'pulse' pulse train, positive valued, even symmetry
P(1) = fractional duty cycle or fractional time High,
0 < P(1) < 1. If not given, P(1) = 1/2.
P(2) = fractional rise and fall time, 0 < P(2) < (1-P(1))/2.
If not given, P(2) = 0.
'bipolar' bipolar pulse train, odd symmetry, zero mean,
P(1) = fractional duty cycle, 0 < P(1) < 1. If not given,
P(1) = 2/3.
P(2) = fractional rise and fall time, 0 < P(2) < (1-P(1))/4.
If not given, P(2) = 0.
'trap' trapezoidal wave, odd symmetry, zero mean,
P = fractional duty cycle, 0 < P < 1. If not given, P = 2/3.
'full' full wave rectified sine wave with the period equal to that
of the sine wave being rectified.
'half' half wave rectified sine wave.
's2s' sine-to-square wave, odd symmetry, zero mean,
P = normalized morph between a sine wave and a square wave,
P = 0 is a sine wave, P = 1 is a square wave,
0 < P < 1 is a flat-topped wave with sinusoidal transitions.
If not given, P = 1/2.
'sine' sine wave, P = harmonic number. If not given, P = 1.
'cosine' cosine wave, P = harmonic number. If not given, P = 1.
'dc' constant or dc value.
--------------------------------------------------------------------------
The created Fourier Series function handle FSfun provides the following:
SYNTAX DESCRIPTION
FSfun(t) evaluate function at the points in numerical array t.
FSfun('tag') return the Tag string stored in FSfun.
FSfun('coef') return Fourier Series coefficient row vector, Kn.
FSfun('period') return function period, T.
FSfun('size') return number of harmonics, N.
FSfun('avg') return DC or average value.
FSfun('msv') return mean square value (Parseval's Theorem).
FSfun('max') return maximum function value.
FSfun('min') return minimum function value.
FSfun('thd') return total harmonic distortion relative to fundamental.
FSfun('one') return one-sided line spectra amplitude vector, i.e.,
[ |K(0)| 2*|K(1)| ... 2*|K(N)| ]
FSfun('phase') return one-sided line spectra phase vector in Degrees,
i.e., [ angle(K(0)) angle(K(1)) ... angle(K(N)) ]
FSfun('sine') return a vector of the one-sided Fourier Series sine
coefficients, sine portion of trigonometric FS form
FSfun('cosine') return a vector of the one-sided Fourier Series cosine
coefficients, cosine portion of trigonometric FS form
FSfun('all') return a structure having field names equal to the above
string arguments with contents equal to associated data
as described above, e.g., out.coef is the coefficients.
FSfun('plot') create a TIME plot of one period of the Fourier Series
FSfun('spectra') create a STEM plot of the one-sided spectra
FSfun( {'area', [Tmin Tmax]} ) returns the area under the Fourier Series
over the time range Tmin <= t <= Tmax.
FSfun( {'plot', [Tmin Tmax]} ) creates the TIME plot over the time range
Tmin <= t <= Tmax.
FSfun( {'spectra', [Nmin Nmax]} ) creates the STEM plot over the positive
harmonic index range Nmin to Nmax.
--------------------------------------------------------------------------0001 function [FSfun,Kn]=fscreate(varargin) 0002 %FSCREATE Create Fourier Series Function. 0003 % [FSfun,Kn] = FSCREATE(...) creates a function handle FSfun containing the 0004 % complex exponential Fourier Series created from its input arguments. The 0005 % optional output Kn is a row vector containing the complex exponential 0006 % Fourier Series coefficients: Kn = [K(-N) ... K(0) ... K(N)] 0007 % 0008 % FSfun = FSCREATE(t,f,N,TYPE) computes the Fourier Series of a signal 0009 % tabulated in the real-valued data vectors t and f. 0010 % t is a monotonically increasing time vector with the signal period being 0011 % T = t(end)-t(1). f(1)=f(end) is required. 0012 % N is the number of harmonics. If N is empty or not given, N = 64. 0013 % TYPE is a string describing the signal type described by the data t and f. 0014 % If TYPE is not given or if TYPE = 'foh' (First Order Hold), the Fourier 0015 % Series coefficients are computed by assuming that the input data is 0016 % piecewise linear between the data points. 0017 % If TYPE = 'zoh' (Zero Order Hold), the Fourier Series coefficients are 0018 % computed by assuming that the input data is piecewise constant between 0019 % the data points. 0020 % The FFT is not used, so no aliasing occurs. 0021 % 0022 % FSfun = FSCREATE(WAVE,T,N,P) returns the Fourier Series of known 0023 % waveforms. WAVE is a character string identifying the desired waveform, 0024 % T is the period, and N is the number of harmonics computed. If N is empty 0025 % or not given, N = 64. P is a vector of optional parameters required for 0026 % some waveforms. See below for waveforms and parameters. 0027 % 0028 % FSfun = FSCREATE(Kn,T) creates a Fourier Series function using the Fourier 0029 % series vector Kn and period T. Kn is a row vector of the complex 0030 % exponential Fourier series coefficients Kn = [K(-N) ... K(0) ... K(N)] 0031 % where K(n) is the n-th term, e.g., K(0) is the DC term, and N is the 0032 % number of harmonics. 0033 % 0034 % FSfun2 = FSCREATE(FSfun,'Op',P) performs the operation specified by the 0035 % string 'Op' on the Fourier Series function handle FSfun, returning a new 0036 % Fourier series function handle. P is a parameter needed for some 'Op'. 0037 % 0038 % Operation DESCRIPTION 0039 % diff differentiate the Fourier Series. 0040 % int integrate the Fourier Series, ignoring the DC or average term 0041 % and returning a Fourier Series having zero average. 0042 % mirror time mirror (Fourier Series of f(-t)). 0043 % smooth apply Blackman window to the Fourier Series (minimize Gibb's) 0044 % trim trim negligible real and imaginary elements. 0045 % even extract the even time part (real parts). 0046 % odd extract the odd time part (imaginary parts). 0047 % halfwave extract the halfwave part, (odd harmonics). 0048 % nodc remove DC or average value term (K(0) = 0). 0049 % dc P = desired DC or average value (K(0) = P). 0050 % add P = amount to add to DC or average value (K(0)=K(0)+P). 0051 % delay P = normalized delay to apply to FSfun, P = 1 is one period. 0052 % delay P = 'odd' or P = 'even' apply delay to achieve given symmetry. 0053 % scale P = amplitude scale factor to apply to FSfun (P*Kn) 0054 % period P = desired period to associate with Fourier Series. 0055 % resize P = number of harmonics desired in result. If P<N terms are 0056 % deleted; if P>N zeros are padded onto Kn. 0057 % tag P = tag character string to store with function handle 0058 % 0059 % FSfun = FSCREATE(FSfunA,'Op',FSfunB) performs the mathematical operation 0060 % specified by the string 'Op' on the underlying time signals in the input 0061 % Fourier Series functions. 'Op' is one of the following '+' (addition), 0062 % '-' (subtraction), '*' (multiplication). The period of FSfunA and FSfunB 0063 % are assumed to be equal. 0064 % 0065 % FSCREATE(...,Tag) where character string Tag is the last input argument, 0066 % stores the string Tag in the function handle for possible retrieval later. 0067 % Tag can be any character string. 0068 % 0069 % FSCREATE('WAVE',T,N,P) Waveforms 0070 % WAVE DESCRIPTION 0071 % 'square' square wave, odd symmetry, zero mean, 0072 % 'saw' sawtooth, positive slope, odd symmetry, zero mean, 0073 % P = fractional fall time, 0 < P < 0.5. If not given, P = 0. 0074 % 'rsaw' reversed sawtooth, same as 'sawtooth' but negative slope 0075 % P = fractional fall time, 0 < P < 0.5. If not given, P = 0. 0076 % 'triangle' triangle, even symmetry, zero mean. 0077 % 'pulse' pulse train, positive valued, even symmetry 0078 % P(1) = fractional duty cycle or fractional time High, 0079 % 0 < P(1) < 1. If not given, P(1) = 1/2. 0080 % P(2) = fractional rise and fall time, 0 < P(2) < (1-P(1))/2. 0081 % If not given, P(2) = 0. 0082 % 'bipolar' bipolar pulse train, odd symmetry, zero mean, 0083 % P(1) = fractional duty cycle, 0 < P(1) < 1. If not given, 0084 % P(1) = 2/3. 0085 % P(2) = fractional rise and fall time, 0 < P(2) < (1-P(1))/4. 0086 % If not given, P(2) = 0. 0087 % 'trap' trapezoidal wave, odd symmetry, zero mean, 0088 % P = fractional duty cycle, 0 < P < 1. If not given, P = 2/3. 0089 % 'full' full wave rectified sine wave with the period equal to that 0090 % of the sine wave being rectified. 0091 % 'half' half wave rectified sine wave. 0092 % 's2s' sine-to-square wave, odd symmetry, zero mean, 0093 % P = normalized morph between a sine wave and a square wave, 0094 % P = 0 is a sine wave, P = 1 is a square wave, 0095 % 0 < P < 1 is a flat-topped wave with sinusoidal transitions. 0096 % If not given, P = 1/2. 0097 % 'sine' sine wave, P = harmonic number. If not given, P = 1. 0098 % 'cosine' cosine wave, P = harmonic number. If not given, P = 1. 0099 % 'dc' constant or dc value. 0100 % 0101 %-------------------------------------------------------------------------- 0102 % The created Fourier Series function handle FSfun provides the following: 0103 % 0104 % SYNTAX DESCRIPTION 0105 % 0106 % FSfun(t) evaluate function at the points in numerical array t. 0107 % 0108 % FSfun('tag') return the Tag string stored in FSfun. 0109 % FSfun('coef') return Fourier Series coefficient row vector, Kn. 0110 % FSfun('period') return function period, T. 0111 % FSfun('size') return number of harmonics, N. 0112 % FSfun('avg') return DC or average value. 0113 % FSfun('msv') return mean square value (Parseval's Theorem). 0114 % FSfun('max') return maximum function value. 0115 % FSfun('min') return minimum function value. 0116 % FSfun('thd') return total harmonic distortion relative to fundamental. 0117 % FSfun('one') return one-sided line spectra amplitude vector, i.e., 0118 % [ |K(0)| 2*|K(1)| ... 2*|K(N)| ] 0119 % FSfun('phase') return one-sided line spectra phase vector in Degrees, 0120 % i.e., [ angle(K(0)) angle(K(1)) ... angle(K(N)) ] 0121 % FSfun('sine') return a vector of the one-sided Fourier Series sine 0122 % coefficients, sine portion of trigonometric FS form 0123 % FSfun('cosine') return a vector of the one-sided Fourier Series cosine 0124 % coefficients, cosine portion of trigonometric FS form 0125 % FSfun('all') return a structure having field names equal to the above 0126 % string arguments with contents equal to associated data 0127 % as described above, e.g., out.coef is the coefficients. 0128 % 0129 % FSfun('plot') create a TIME plot of one period of the Fourier Series 0130 % FSfun('spectra') create a STEM plot of the one-sided spectra 0131 % 0132 % FSfun( {'area', [Tmin Tmax]} ) returns the area under the Fourier Series 0133 % over the time range Tmin <= t <= Tmax. 0134 % FSfun( {'plot', [Tmin Tmax]} ) creates the TIME plot over the time range 0135 % Tmin <= t <= Tmax. 0136 % FSfun( {'spectra', [Nmin Nmax]} ) creates the STEM plot over the positive 0137 % harmonic index range Nmin to Nmax. 0138 %-------------------------------------------------------------------------- 0139 0140 % D.C. Hanselman, University of Maine, Orono, ME 04469 0141 % MasteringMatlab@yahoo.com 0142 % Mastering MATLAB 7 0143 % 2006-01-01, revised 2006-01-03, 2006-03-08, 2006-05-18, 2006-07-03 0144 0145 % parse inputs 0146 N=64; 0147 FSData.tag=''; 0148 narg=nargin; 0149 if narg<2 0150 error('At Least Two Input Arguments Required.') 0151 end 0152 if narg>2 && ischar(varargin{end}) && ... 0153 ~any(strcmpi(varargin{end},{'odd','even'})) && ... 0154 ~any(strcmpi(varargin{end},{'foh','zoh'})) && ... 0155 (ischar(varargin{end-1}) && ~isequal(lower(varargin{end-1}),'tag')) 0156 FSData.tag=varargin{end}; % grab tag provided 0157 varargin(end)=[]; % strip tag from input 0158 narg=narg-1; 0159 end 0160 switch class(varargin{1}) % find class of first input argument 0161 case {'double' 'single'} 0162 if isequal(numel(varargin{1}),numel(varargin{2})) % fscreate(t,f,N,TYPE) 0163 t=varargin{1}(:); 0164 f=varargin{2}(:); 0165 if narg>=3 && isnumeric(varargin{3}) 0166 N=varargin{3}; 0167 if fix(N)~=N || abs(N)~=N || N<1 0168 error('FSCREATE:InputError','N Must be a Positive Integer.') 0169 end 0170 end 0171 if ischar(varargin{end}) % Type is specified 0172 Tref={'foh' 'zoh'}; 0173 idx=strncmpi(varargin{end},Tref,1); 0174 Type=char(Tref(idx)); 0175 if isempty(Type) 0176 error('FSCREATE:InputError','Unknown TYPE Argument.') 0177 end 0178 else 0179 Type='foh'; 0180 end 0181 Kn=local_getFS(t,f,N,Type); 0182 FSData.Kn=Kn; 0183 FSData.T=t(end)-t(1); 0184 FSfun=@(t) FourierSeries(FSData,t); 0185 else % fscreate(Kn,T) 0186 Kn=varargin{1}(:).'; 0187 N=(length(Kn)-1)/2; 0188 if N~=fix(N) || max(abs(conj(Kn(N:-1:1))-Kn(N+2:end)))>sqrt(eps) 0189 error('FSCREATE:InputError', ... 0190 ['First Input Argument must be a valid Fourier Series\n', ... 0191 'vector containing the complex exponential series\n', ... 0192 'coefficients in increasing harmonic order with\n', ... 0193 'K(-n)=conj(K(n)) where n is the harmonic number.']); 0194 end 0195 T=varargin{2}; 0196 if numel(T)~=1 || T<=0 || ~isreal(T) 0197 error('Second Input Argument Must be the Positive Scalar Period.') 0198 end 0199 FSData.Kn=Kn; 0200 FSData.T=T; 0201 FSfun=@(t) FourierSeries(FSData,t); 0202 end 0203 case 'char' % fscreate('WAVE',T,N,P) 0204 wave=varargin{1}; 0205 T=varargin{2}; 0206 if narg==4 0207 P=varargin{4}; 0208 else 0209 P=[]; 0210 end 0211 if narg>=3 && ~isempty(varargin{3}) 0212 N=varargin{3}; 0213 if fix(N)~=N || abs(N)~=N || N<1 0214 error('FSCREATE:InputError','N Must be a Positive Integer.') 0215 end 0216 end 0217 Kn=local_getwave(wave,T,N,P); 0218 FSData.Kn=Kn; 0219 FSData.T=T; 0220 FSfun=@(t) FourierSeries(FSData,t); 0221 case 'function_handle' 0222 if ~ischar(varargin{2}) 0223 error('FSCREATE:InputError','Second Input Argument Must be a String.') 0224 end 0225 if narg==3 && ... 0226 isa(varargin{3},'function_handle') % fscreate(FSfunA,'Op',FSfunB) 0227 0228 FSfunA=varargin{1}; 0229 NA=FSfunA('size'); 0230 KnA=FSfunA('coef'); 0231 FSData.T=FSfunA('period'); 0232 0233 FSfunB=varargin{3}; 0234 NB=FSfunB('size'); 0235 KnB=FSfunB('coef'); 0236 0237 zAB=zeros(1,NA-NB); 0238 zBA=zeros(1,NB-NA); 0239 switch lower(varargin{2}) 0240 case '+' 0241 Kn=[zBA KnA zBA]+[zAB KnB zAB]; 0242 case '-' 0243 Kn=[zBA KnA zBA]-[zAB KnB zAB]; 0244 case '*' 0245 Kn=conv(KnA,KnB); 0246 iDC=(length(Kn)+1)/2; 0247 N=max(NA,NB); 0248 Kn=Kn(iDC-N:iDC+N); 0249 otherwise 0250 error('FSCREATE:InputError','Unknown Mathematical Operator.') 0251 end 0252 FSData.Kn=Kn; 0253 FSfun=@(t) FourierSeries(FSData,t); 0254 0255 else % fscreate(FSfun,'Op') or fscreate(FSfun,'Op',P) 0256 FSfun=varargin{1}; 0257 Kn=FSfun('coef'); 0258 N=(length(Kn)-1)/2; 0259 iDC=N+1; 0260 T=FSfun('period'); 0261 wo=2*pi/T; 0262 if narg==3 0263 P=varargin{3}; 0264 end 0265 switch lower(varargin{2}(1:min(2,length(varargin{2})))) 0266 case 'ta' % tag 0267 if ischar(P) 0268 FSData.tag=P; 0269 else 0270 error('FSCREATE:InputError','Character String Tag Required.') 0271 end 0272 case 'di' % differentiate 0273 Kn=1j*wo*(-N:N).*Kn; 0274 case 'in' % integrate 0275 nn=-N:N; 0276 idx=nn~=0; 0277 iKn=zeros(size(Kn)); 0278 iKn(idx)=Kn(idx)./(nn(idx)*wo*1i); 0279 iKn(iDC)=0; 0280 Kn=iKn; 0281 case 'mi' % time mirror 0282 Kn=Kn(end:-1:1); 0283 case 'sm' % Exact Blackman smoothing 0284 m=linspace(-pi,pi,length(Kn)); 0285 Kn=Kn.*(7938 + 9240*cos(m) + 1430*cos(2*m))/18608; 0286 case 'tr' % trim negigible elements 0287 tol=sqrt(eps); 0288 mKn=abs(Kn); 0289 rKn=abs(real(Kn)); 0290 iKn=abs(imag(Kn)); 0291 b=mKn<tol*max(mKn); 0292 Kn(b)=0; % elements with small magnitude 0293 b=rKn<tol*max(rKn) | rKn<tol*iKn; 0294 Kn(b)=1i*imag(Kn(b)); % elements with small real part 0295 b=iKn<tol*max(iKn) | iKn<tol*rKn; 0296 Kn(b)=real(Kn(b)); % elements with small imaginary part 0297 case 'ev' % Even part 0298 Kn=real(Kn); 0299 case 'od' % Odd part 0300 Kn=1i*imag(Kn); 0301 case 'ha' % Halfwave part 0302 idx=rem(N,2)+1:2:length(Kn); 0303 Kn(idx)=0; 0304 case 'no' % no DC part 0305 Kn(iDC)=0; 0306 case 'dc' % set DC value 0307 if narg==2 || numel(varargin{3})>1 0308 error('FSCREATE:InputError','Scalar Third Argument P Required.') 0309 end 0310 Kn(iDC)=real(P); 0311 case 'ad' % set DC value 0312 if narg==2 || numel(varargin{3})>1 0313 error('FSCREATE:InputError','Scalar Third Argument P Required.') 0314 end 0315 Kn(iDC)=Kn(iDC)+real(P); 0316 case 'de' % delay 0317 if narg==2 0318 error('FSCREATE:InputError','Third Argument P Required.') 0319 end 0320 if ischar(P) 0321 a1=angle(Kn(iDC+1))/(2*pi); 0322 switch lower(P(1)) 0323 case 'o' 0324 d=a1+1/4; 0325 case 'e' 0326 d=a1; 0327 otherwise 0328 error('FSCREATE:InputError','Unknown Delay Requested.') 0329 end 0330 else 0331 d=P(1); 0332 end 0333 Kn=exp(-2j*pi*(-N:N)*d).*Kn; 0334 case 'sc' % scale amplitude 0335 if narg==2 || numel(varargin{3})>1 0336 error('FSCREATE:InputError','Scalar Third Argument P Required.') 0337 end 0338 Kn=P*Kn; 0339 case 'pe' % set period 0340 if narg==2 || numel(varargin{3})>1 0341 error('FSCREATE:InputError','Scalar Third Argument P Required.') 0342 end 0343 T=varargin{3}; 0344 case 're' % resize 0345 if narg==2 || numel(varargin{3})>1 0346 error('FSCREATE:InputError','Scalar Third Argument P Required.') 0347 elseif fix(P)~=P || abs(P)~=P || P<1 0348 error('FSCREATE:InputError','P Must be a Positive Integer.') 0349 end 0350 if P>N % pad with zeros 0351 zpadd=zeros(1,P-N); 0352 Kn=[zpadd Kn zpadd]; 0353 elseif P<N % delete excess terms 0354 Kn=Kn(iDC-P:iDC+P); 0355 else % P=N no work 0356 return 0357 end 0358 otherwise 0359 error('FSCREATE:InputError','Unknown Operation Requested.') 0360 end 0361 FSData.Kn=Kn; 0362 FSData.T=T; 0363 FSfun=@(t) FourierSeries(FSData,t); 0364 end 0365 otherwise 0366 error('Unknown First Input Argument.') 0367 end % END of Primary Function 0368 %-------------------------------------------------------------------------- 0369 %-------------------------------------------------------------------------- 0370 function Kn=local_getFS(t,f,N,type) % fscreate(t,f,N,type) Subfunction 0371 T=t(end)-t(1); 0372 if T<=0 0373 error('Period Must be Positive.') 0374 end 0375 if ~isreal(f) 0376 error('Input Must be Real-Valued.') 0377 end 0378 if abs(f(1)-f(end))> eps*(1+abs(max(f)-min(f))) 0379 error('Input Must Describe Exactly One Period of the Function.') 0380 end 0381 t=t/T; 0382 wo=2*pi; 0383 df=diff(f); 0384 dt=diff(t); 0385 if any(dt<0) 0386 error('Time Points Must be Nondecreasing.') 0387 end 0388 switch type 0389 case 'foh' 0390 Nl=length(f)-1; 0391 idx=dt~=0; 0392 m=zeros(Nl,1); 0393 m(idx)=df(idx)./dt(idx); 0394 b=f(1:end-1)-m.*t(1:end-1); 0395 n=1:N; 0396 mm=repmat(m,1,N); 0397 bb=repmat(b,1,N); 0398 t1=repmat(t(1:end-1),1,N); 0399 t2=repmat(t(2:end),1,N); 0400 nn=repmat(n,Nl,1); 0401 Fp=sum((1j*(mm.*t2+bb)./(nn*wo)+mm./(nn*wo).^2).*exp(-1j*nn*wo.*t2) - ... 0402 (1j*(mm.*t1+bb)./(nn*wo)+mm./(nn*wo).^2).*exp(-1j*nn*wo.*t1)); 0403 Kn=[conj(Fp(end:-1:1)) trapz(t,f) Fp]; 0404 case 'zoh' 0405 n=1:N; 0406 Ni=length(f)-1; 0407 ff=repmat(f(1:end-1),1,N); 0408 t1=repmat(t(1:end-1),1,N); 0409 t2=repmat(t(2:end),1,N); 0410 nn=repmat(n,Ni,1); 0411 0412 Fp=sum(1j*ff./(nn*wo).*(exp(-1j*nn*wo.*t2) - exp(-1j*nn*wo.*t1))); 0413 Kn=[conj(Fp(end:-1:1)) trapz(t,f) Fp]; 0414 end 0415 %-------------------------------------------------------------------------- 0416 function Kn=local_getwave(wave,T,N,P) % fscreate('WAVE',T,N,P) Subfunction 0417 lenP=length(P); 0418 lkn=2*N+1; 0419 Kn=zeros(1,lkn); 0420 0421 switch lower(wave(1:min(3,length(wave)))) 0422 case 'squ' % squarewave 0423 Kn(N+2:2:lkn)=-2j*(((1:2:N)*pi).^(-1)); 0424 Kn(1:N)=conj(Kn(lkn:-1:N+2)); 0425 case 'tri' % triangle 0426 Kn(N+2:2:lkn)=4*(((1:2:N)*pi).^(-2)); 0427 Kn(1:N)=conj(Kn(lkn:-1:N+2)); 0428 case 'ful' % fullwave rectified sine 0429 Kn(N+1)=2/pi; 0430 n=2:2:N; 0431 Kn(N+3:2:lkn)=2./(pi*(1-n.^2)); 0432 Kn(1:N)=conj(Kn(lkn:-1:N+2)); 0433 case 'hal' % halfwave rectified sine 0434 Kn(N+1)=1/pi; 0435 Kn(N+2)=-j/4; 0436 n=2:2:N; 0437 Kn(N+3:2:lkn)=1./(pi*(1-n.^2)); 0438 Kn(1:N)=conj(Kn(lkn:-1:N+2)); 0439 case 'dc' % dc or constant 0440 Kn(N+1)=1; 0441 case 'sin' % sine wave 0442 if lenP==1 0443 n=P; 0444 else 0445 n=1; 0446 end 0447 if fix(n)~=n || n<1 || n>N 0448 error('Harmonic Number Must be an Integer Between 1 and N.') 0449 end 0450 Kn=zeros(1,N); 0451 Kn(n)=-1i/2; 0452 Kn=[conj(Kn(end:-1:1)) 0 Kn]; 0453 case 'cos' % cosine wave 0454 if lenP==1 0455 n=P; 0456 else 0457 n=1; 0458 end 0459 if fix(n)~=n || n<1 || n>N 0460 error('Harmonic Number Must be an Integer Between 1 and N.') 0461 end 0462 Kn=zeros(1,N); 0463 Kn(n)=1/2; 0464 Kn=[Kn(end:-1:1) 0 Kn]; 0465 case 'saw' % sawtooth 0466 if lenP==0 % ideal sawtooth 0467 Kn(N+2:lkn)=2j*(((1:N)*2*pi).^(-1)); 0468 Kn(1:N)=conj(Kn(lkn:-1:N+2)); 0469 else % finite fall time sawtooth 0470 a=min(.4999,P(1)); 0471 if a~=P(1) 0472 warning('FSCREATE:InputError',... 0473 'Requested Fall Time Must be Less Than 1/2.') 0474 end 0475 t=[0 a .5 1-a 1]; 0476 f=[0 -1 0 1 0]; 0477 [dum,Kn]=fscreate(t,f,N); 0478 end 0479 case 'rsa' % reverse sawtooth 0480 [dum,Kn]=fscreate('saw',T,N,P); 0481 Kn=Kn(end:-1:1); 0482 case 'tra' % trapezoid 0483 if lenP==0 0484 p=2/3; 0485 else 0486 p=min(abs(P(1)),.999); 0487 end 0488 if lenP>0 && p~=P(1) 0489 warning('FSCREATE:InputError',... 0490 'Requested Duty Cycle Must Be Less Than 1.') 0491 end 0492 a=(1-p)/2; 0493 b=1-a; 0494 npi=(1:2:N)*pi; 0495 jnpia=1j*a*npi; 0496 jnpib=1j*b*npi; 0497 Kn(N+2:2:lkn)=(( (1+1j*a*npi).*exp(-jnpia)... 0498 - (1+jnpib).*exp(jnpia)... 0499 + 1j*npi).*(npi.^(-2))./a... 0500 + ( (1/a-1)*exp(jnpia)... 0501 - exp(-jnpia) -1/a)*1j./npi); 0502 Kn(1:N)=conj(Kn(lkn:-1:N+2)); 0503 case 'pul' % Pulsetrain 0504 if lenP==0 0505 p=1/2; 0506 else 0507 p=min(abs(P(1)),.999); 0508 end 0509 if lenP>0 && p~=P(1) 0510 warning('FSCREATE:InputError',... 0511 'Requested Duty Cycle Must Be Less Than 1.') 0512 end 0513 if lenP<2 || P(2)==0 % ideal zero rise time pulse 0514 Kn(N+1)=p; 0515 arg=pi*p*(1:N); 0516 Kn(N+2:lkn)=p*sin(arg)./(arg+eps); 0517 Kn(1:N)=conj(Kn(lkn:-1:N+2)); 0518 else % finite rise time pulse 0519 r=min((1-p)/2-10*eps,abs(P(2))); 0520 if r~=P(2) 0521 warning('FSCREATE:InputError',... 0522 'Requested Rise Time Too Long.') 0523 end 0524 a=p/2; 0525 t=[0 a a+r 0.5 1-a-r 1-a 1]; 0526 f=[1 1 0 0 0 1 1]; 0527 [dum,Kn]=fscreate(t,f,N); 0528 end 0529 case 'bip' % Bipolar Pulsetrain 0530 if lenP==0 0531 p=2/3; 0532 else 0533 p=min(abs(P(1)),.999); 0534 end 0535 if lenP>0 && p~=P(1) 0536 warning('FSCREATE:InputError',... 0537 'Requested Duty Cycle Must Be Less Than 1.') 0538 end 0539 if lenP<2 || P(2)==0 % ideal zero rise time pulsetrain 0540 a=(1-p)/2; b=1-a; 0541 jnpi=(1:2:N)*pi*1j; 0542 Kn(N+2:2:lkn)=(exp(-jnpi*b)-exp(-jnpi*a))./(-jnpi); 0543 Kn(1:N)=conj(Kn(lkn:-1:N+2)); 0544 else % finite rise time pulse 0545 r=min((1-p)/4-10*eps,abs(P(2))); 0546 if r~=P(2) 0547 warning('FSCREATE:InputError',... 0548 'Requested Rise Time Too Long.') 0549 end 0550 a=p/4; 0551 b=.25; 0552 c=.75; 0553 t=[0 b-a-r b-a b+a b+a+r 0.5 c-a-r c-a c+a c+a+r 1]; 0554 f=[0 0 1 1 0 0 0 -1 -1 0 0]; 0555 [dum,Kn]=fscreate(t,f,N); 0556 end 0557 case 's2s' % sine to square morph 0558 if lenP==0 0559 p=1/2; 0560 else 0561 p=P(1); 0562 end 0563 if p<0 || p>1 0564 error('FSCREATE:InputError',... 0565 'Morphing Parameter Must be Between 0 and 1.') 0566 end 0567 if p==0 % sine wave 0568 Kn=zeros(1,N); 0569 Kn(1)=-1i/2; 0570 Kn=[conj(Kn(end:-1:1)) 0 Kn]; 0571 elseif p==1 % square wave 0572 Kn(N+2:2:lkn)=-2j*(((1:2:N)*pi).^(-1)); 0573 Kn(1:N)=conj(Kn(lkn:-1:N+2)); 0574 else % morph between sine and square 0575 kp=zeros(1,N); 0576 b=(1-p)/2; 0577 m=1:2:N; 0578 bm=b*m; 0579 kp(m)=(2*b/pi)*( 2i*bm.*exp(-1i*bm*pi)-1 )./(4*bm.^2-1)... 0580 - 2i*cos(bm*pi)./(m*pi)... 0581 + (2*b/pi)*( 1+ 2i*bm.*exp(1i*bm*pi))./(4*bm.^2-1); 0582 Kn=[conj(kp(end:-1:1)) 0 kp]; 0583 end 0584 otherwise 0585 error('FSCREATE:InputError','Unknown Waveform Input.') 0586 end 0587 %-------------------------------------------------------------------------- 0588 % Fourier Series Function 0589 %-------------------------------------------------------------------------- 0590 function y=FourierSeries(FSData,t) 0591 % Fourier Series Evaluation and Manipulation using function handles 0592 % 0593 % This function gets called when FSfun(...) is issued where FSfun is a 0594 % function handle returned by FSCREATE. 0595 0596 Kn=FSData.Kn; % FS coefficient vector 0597 wo=2*pi/FSData.T; % fundamental frequency 0598 N=(length(Kn)-1)/2; % highest harmonic 0599 iDC=N+1; % index of dc component 0600 0601 if isnumeric(t) && isreal(t) % Evaluate Fourier Series 0602 nwo=wo*(1:N); % positive harmonic indices 0603 ko=real(Kn(iDC)); % average value 0604 tsize=size(t); 0605 Knp=Kn(iDC+1:end).'; % positive frequency coefs 0606 y=ko+2*(real(exp(1j*t(:)*nwo)*Knp))'; 0607 y=reshape(y,tsize); 0608 elseif iscell(t) && ischar(t{1}) % FSfun( {'command', [range]} ) 0609 if length(t)~=2 || ~isnumeric(t{2}) || numel(t{2})~=2 || ~isreal(t{2}) 0610 msg=['Two Element Cell Required. '... 0611 'Second Cell Must be a Two Element Numeric Vector.']; 0612 error(msg) 0613 end 0614 switch lower(t{1}(1:min(3,length(t{1})))) 0615 case 'are' 0616 n=1:N; 0617 nwo=n*wo; 0618 Ikn=Kn(iDC+n)./(1i*nwo); % Kn/(jn*wo) is integral 0619 tmm=t{2}; 0620 y=diff(2*(real(exp(1j*tmm(:)*nwo)*Ikn.')))+Kn(iDC)*diff(tmm); 0621 case 'ste' 0622 0623 if any(fix(t{2})~=t{2}) || any(t{2}<0) 0624 error('Harmonic Indices Must be Nonnegative Integers.') 0625 end 0626 xdata=(min(t{2}):min(max(t{2}),N))+1; 0627 subplot(2,1,1) 0628 ydata=abs([Kn(iDC) 2*Kn(iDC+1:end)]); 0629 ydata=ydata(xdata); 0630 stem(xdata,ydata) 0631 ylabel('Amplitude') 0632 xlabel('Harmonic Index') 0633 title('Fourier Series Line Spectra Plot') 0634 ydata=angle(Kn(iDC:end))*180/pi; 0635 ydata=ydata(xdata); 0636 subplot(2,1,2) 0637 stem(xdata,ydata) 0638 ylabel('Phase') 0639 xlabel('Harmonic Index') 0640 case 'plo' 0641 tt=linspace(t{2}(1),t{2}(2),min(max(5*N,100),500)); 0642 nwo=wo*(1:N); 0643 ko=real(Kn(iDC)); 0644 Knp=Kn(iDC+1:end).'; 0645 ydata=ko+2*(real(exp(1j*tt(:)*nwo)*Knp))'; 0646 plot(tt,ydata) 0647 title('Fourier Series Plot.') 0648 otherwise 0649 error('Uknown Input Argument.') 0650 end 0651 0652 elseif ~ischar(t) 0653 error('Unknown Input Argument.') 0654 else % FSfun('string') 0655 switch lower(t(1:min(3,length(t)))) 0656 case 'tag' 0657 y=FSData.tag; 0658 case {'coe' 'kn'} % coefficient vector 0659 y=Kn; 0660 case {'per' 't'} % period 0661 y=FSData.T; 0662 case 'siz' % size or highest harmonic 0663 y=N; 0664 case {'dc' 'ave' 'avg'} % dc or average value 0665 y=Kn(iDC); 0666 case 'msv' % mean square value 0667 y=real(Kn*Kn'); 0668 case 'max' % maximum or peak value 0669 nwo=wo*(1:N); 0670 Knp=Kn(iDC+1:end).'; 0671 t=linspace(-0.1,1.2,11*N)*FSData.T; 0672 [dum,idx]=max(real(Kn(iDC))+2*(real(exp(1j*t(:)*nwo)*Knp))); 0673 t=linspace(t(idx-1),t(idx+1),101); 0674 y=max(real(Kn(iDC))+2*(real(exp(1j*t(:)*nwo)*Knp))); 0675 case 'min' % minimum value 0676 nwo=wo*(1:N); 0677 Knp=Kn(iDC+1:end).'; 0678 t=linspace(-0.1,1.2,11*N)*FSData.T; 0679 [dum,idx]=min(real(Kn(iDC))+2*(real(exp(1j*t(:)*nwo)*Knp))); 0680 t=linspace(t(idx-1),t(idx+1),101); 0681 y=min(real(Kn(iDC))+2*(real(exp(1j*t(:)*nwo)*Knp))); 0682 case 'thd' % total harmonic distortion 0683 idx=iDC+[-1 1]; 0684 fn=Kn(idx); 0685 Kn(idx)=[]; 0686 y=sqrt(real(Kn*Kn')/real(fn*fn')); 0687 case 'one' % one-sided line amplitude spectra 0688 y=abs([Kn(iDC) 2*Kn(iDC+1:end)]); 0689 case 'pha' % one-sided line phase spectra 0690 y=angle(Kn(iDC:end))*180/pi; 0691 case 'sin' % trig sine coefficients 0692 y=-2*imag(Kn(iDC+1:end)); 0693 case 'cos' % trig cosine coefficients 0694 y=2*real(Kn(iDC+1:end)); 0695 case 'all' % return all data in a structure 0696 y.tag=FSData.tag; 0697 y.coef=Kn; 0698 y.period=FSData.T; 0699 y.size=N; 0700 y.avg=Kn(iDC); 0701 y.msv=real(Kn*Kn'); 0702 0703 nwo=wo*(1:N); 0704 Knp=Kn(iDC+1:end).'; 0705 t1=linspace(-0.1,1.2,11*N)*FSData.T; 0706 [dum,idx]=max(real(Kn(iDC))+2*(real(exp(1j*t1(:)*nwo)*Knp))); 0707 t2=linspace(t1(idx-1),t1(idx+1),101); 0708 y.max=max(real(Kn(iDC))+2*(real(exp(1j*t2(:)*nwo)*Knp))); 0709 0710 [dum,idx]=min(real(Kn(iDC))+2*(real(exp(1j*t1(:)*nwo)*Knp))); 0711 t2=linspace(t1(idx-1),t1(idx+1),101); 0712 y.min=min(real(Kn(iDC))+2*(real(exp(1j*t2(:)*nwo)*Knp))); 0713 0714 idx=iDC+[-1 1]; 0715 fn=Kn(idx); 0716 Kn(idx)=[]; 0717 y.thd=sqrt(real(Kn*Kn')/real(fn*fn')); 0718 0719 y.one=abs([Kn(iDC) 2*Kn(iDC+1:end)]); 0720 y.phase=angle(Kn(iDC:end))*180/pi; 0721 y.sine=-2*imag(Kn(iDC+1:end)); 0722 y.cosine=2*real(Kn(iDC+1:end)); 0723 0724 case {'spe' 'ste'} % create stem plot 0725 xdata=0:N; 0726 subplot(2,1,1) 0727 ydata=abs([Kn(iDC) 2*Kn(iDC+1:end)]); 0728 stem(xdata,ydata) 0729 ylabel('Amplitude') 0730 xlabel('Harmonic Index') 0731 title('Fourier Series Line Spectra Plot') 0732 ydata=angle(Kn(iDC:end))*180/pi; 0733 subplot(2,1,2) 0734 stem(xdata,ydata) 0735 ylabel('Phase') 0736 xlabel('Harmonic Index') 0737 case 'plo' % create time plot 0738 t=linspace(0,FSData.T,min(max(5*N,100),500)); 0739 nwo=wo*(1:N); 0740 ko=real(Kn(iDC)); 0741 Knp=Kn(iDC+1:end).'; 0742 ydata=ko+2*(real(exp(1j*t(:)*nwo)*Knp))'; 0743 plot(t,ydata) 0744 title('Fourier Series Plot.') 0745 otherwise 0746 error('Uknown Input Argument.') 0747 end 0748 end 0749 %--------------------------------------------------------------------------