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)

## CROSS-REFERENCE INFORMATION

This function calls:
 simpson Numerical integration with the Simpson method diff Difference and approximate derivative. error Display message and abort function. ifft Inverse discrete Fourier transform. trapz Trapezoidal numerical integration.
This function is called by:
 dat2dspec Estimates the directional wave spectrum from timeseries dspec2char Evaluates directional spectral characteristics

## 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 %
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