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
Numerical integration with the Simpson method | |
Difference and approximate derivative. | |
Display message and abort function. | |
Inverse discrete Fourier transform. | |
Trapezoidal numerical integration. |
Estimates the directional wave spectrum from timeseries | |
Evaluates directional spectral characteristics |
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
Comments or corrections to the WAFO group