Home > wafo > multidim > fourier.m

fourier

PURPOSE ^

Returns Fourier coefficients.

SYNOPSIS ^

[a,b]=fourier(t,x,T,M,N);

DESCRIPTION ^

 FOURIER Returns Fourier coefficients. 
  
   CALL:  [a,b] = fourier(t,x,T,M); 
  
     a,b  = Fourier coefficients size M x P 
      t   = vector with N values indexed from 1 to N. 
      x   = vector or matrix of column vectors with data points size N x P. 
      T   = primitive period of signal, i.e., smallest period.  
            (default T = diff(t([1,end])) 
      M-1 = no of harmonics desired (default M = N) 
      N   = no of data points (default length(t)) 
  
  FOURIER finds the coefficients for a Fourier series representation 
  of the signal x(t) (given in digital form).  It is assumed the signal 
  is periodic over T.  N is the number of data points, and M-1 is the 
  number of coefficients. 
  
  The signal can be estimated by using M-1 harmonics by: 
                     M 
  x(i) = 0.5*a(1) + sum [a(n)*c(n,i) + b(n)*s(n,i)] 
                    n=2 
  where 
    c(n,i) = cos(2*pi*(n-1)*t(i)/T) 
    s(n,i) = sin(2*pi*(n-1)*t(i)/T) 
  
  Note that a(1) is the "dc value". 
  Remaining values are a(2), a(3), ... , a(M). 
   
  Example 
   T = 2*pi;M=5; 
   t = linspace(0,4*T).'; x = sin(t); 
   [a,b] = fourier(t,x,T,M) 
  
  See also  fft

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [a,b]=fourier(t,x,T,M,N); 
002 %FOURIER Returns Fourier coefficients. 
003 % 
004 %  CALL:  [a,b] = fourier(t,x,T,M); 
005 % 
006 %    a,b  = Fourier coefficients size M x P 
007 %     t   = vector with N values indexed from 1 to N. 
008 %     x   = vector or matrix of column vectors with data points size N x P. 
009 %     T   = primitive period of signal, i.e., smallest period.  
010 %           (default T = diff(t([1,end])) 
011 %     M-1 = no of harmonics desired (default M = N) 
012 %     N   = no of data points (default length(t)) 
013 % 
014 % FOURIER finds the coefficients for a Fourier series representation 
015 % of the signal x(t) (given in digital form).  It is assumed the signal 
016 % is periodic over T.  N is the number of data points, and M-1 is the 
017 % number of coefficients. 
018 % 
019 % The signal can be estimated by using M-1 harmonics by: 
020 %                    M 
021 % x(i) = 0.5*a(1) + sum [a(n)*c(n,i) + b(n)*s(n,i)] 
022 %                   n=2 
023 % where 
024 %   c(n,i) = cos(2*pi*(n-1)*t(i)/T) 
025 %   s(n,i) = sin(2*pi*(n-1)*t(i)/T) 
026 % 
027 % Note that a(1) is the "dc value". 
028 % Remaining values are a(2), a(3), ... , a(M). 
029 %  
030 % Example 
031 %  T = 2*pi;M=5; 
032 %  t = linspace(0,4*T).'; x = sin(t); 
033 %  [a,b] = fourier(t,x,T,M) 
034 % 
035 % See also  fft 
036  
037  
038 %History: 
039 % Revised pab 22.06.2001 
040 %  -updated help header to wafo style. 
041 %  -vectorized for loop. 
042 %  -added default values for N,M, and T 
043 % 
044 % ME 244L Dynamic Systems and Controls Laboratory 
045 % R.G. Longoria 9/1998 
046 % 
047  
048 t = t(:); 
049 if nargin<5|isempty(N), N = length(t);end 
050 if nargin<4|isempty(M), M = N;        end 
051 if nargin<3|isempty(T), T = diff(t([1,end]));end 
052 szx = size(x); 
053 P=1; 
054 if prod(szx)==N, 
055   x = x(:); 
056 elseif szx(1)==N 
057   P = prod(szx(2:end)); 
058 else 
059   error('Wrong size of x') 
060 end 
061  
062  
063 switch 0 
064   case -1, 
065        % Define the vectors for computing the Fourier coefficients 
066   % 
067   a = zeros(M,P); 
068   b = zeros(M,P); 
069   a(1,:) = simpson(x); 
070  
071   % 
072   % Compute M-1 more coefficients 
073   tmp  = 2*pi*t(:,ones(1,P))/T; 
074   % tmp =  2*pi*(0:N-1).'/(N-1);  
075   for n1 = 1:M-1, 
076     n = n1+1; 
077     a(n,:) = simpson(x.*cos(n1*tmp)); 
078     b(n,:) = simpson(x.*sin(n1*tmp)); 
079   end 
080    
081   a = 2*a/N; 
082   b = 2*b/N; 
083    
084   case 0, 
085      % 
086   a = zeros(M,P); 
087   b = zeros(M,P); 
088   a(1,:) = trapz(t,x); 
089  
090   % 
091   % Compute M-1 more coefficients 
092   tmp  = 2*pi*t(:,ones(1,P))/T; 
093   % tmp =  2*pi*(0:N-1).'/(N-1);  
094   for n1 = 1:M-1, 
095     n = n1+1; 
096     a(n,:) = trapz(t,x.*cos(n1*tmp)); 
097     b(n,:) = trapz(t,x.*sin(n1*tmp)); 
098   end 
099   a = a/pi; 
100   b = b/pi; 
101    
102   case 1, 
103   % Define the vectors for computing the Fourier coefficients 
104   % 
105   a = zeros(M,P); 
106   b = zeros(M,P); 
107    
108   % 
109   % Compute the dc-level (the a(0) component). 
110   % 
111   % Note: the index has to begin with "1". 
112   % 
113    
114   a(1,:) = sum(x); 
115  
116   % 
117   % Compute M-1 more coefficients 
118   tmp  = 2*pi*t(:,ones(1,P))/T; 
119   % tmp =  2*pi*(0:N-1).'/(N-1);  
120   for n1 = 1:M-1, 
121     n = n1+1; 
122     a(n,:) = sum(x.*cos(n1*tmp)); 
123     b(n,:) = sum(x.*sin(n1*tmp)); 
124   end 
125   a = 2*a/N; 
126   b = 2*b/N; 
127 case 2, 
128    % Define the vectors for computing the Fourier coefficients 
129   % 
130   a = zeros(M,P); 
131   b = zeros(M,P); 
132   a(1,:) = trapz(x); 
133  
134   % 
135   % Compute M-1 more coefficients 
136   tmp  = 2*pi*t(:,ones(1,P))/T; 
137   % tmp =  2*pi*(0:N-1).'/(N-1);  
138   for n1 = 1:M-1, 
139     n = n1+1; 
140     a(n,:) = trapz(x.*cos(n1*tmp)); 
141     b(n,:) = trapz(x.*sin(n1*tmp)); 
142   end 
143    
144   a = 2*a/N; 
145   b = 2*b/N; 
146 case 3 
147   % Alternative:  faster for large M, but gives different results than above. 
148   nper = diff(t([1 end]))/T; %No of periods given 
149   if nper == round(nper),  
150     N1 = N/nper; 
151   else 
152     N1 = N; 
153   end 
154  
155   % Fourier coefficients by fft 
156   Fcof1 = 2*ifft(x(1:N1,:),[],1); 
157   Pcor = [1; exp(sqrt(-1)*[1:M-1].'*t(1))]; % correction term to get 
158                                               % the correct integration limits 
159   Fcof = Fcof1(1:M,:).*Pcor(:,ones(1,P)); 
160   a = real(Fcof(1:M,:)); 
161   b = imag(Fcof(1:M,:)); 
162 end 
163 return 
164  
165 %Old call: kept just in case 
166  
167 % 
168 % Compute the dc-level (the a(0) component). 
169 % 
170 % Note: the index has to begin with "1". 
171 %8 
172 a(1) = 2*sum(x)/N; 
173  
174 % 
175 % Compute M-1 more coefficients 
176 for n = 2:M, 
177   sumcos=0.0; 
178   sumsin=0.0; 
179   for i=1:N, 
180     sumcos = sumcos + x(i)*cos(2*(n-1)*pi*t(i)/T); 
181     sumsin = sumsin + x(i)*sin(2*(n-1)*pi*t(i)/T); 
182   end 
183   a(n) = 2*sumcos/N; 
184   b(n) = 2*sumsin/N; 
185 end 
186  
187  
188 return 
189  
190  
191

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003