Mo Logo [Home] [Lexikon] [Aufgaben] [Tests] [Kurse] [Begleitmaterial] [Hinweise] [Mitwirkende] [Publikationen]
Mathematik-Online-Aufgabensammlung: Lösung zu

Aufgabe 1545: Programm zur Approximation der Fourier-Koeffizienten

A B C D E F G H I J K L M N O P Q R S T U V W X Y Z
Schreiben Sie ein MATLAB-Programm c=fourier_coefficients(fct,k,tol), das die Fourier-Koeffizienten $ c_0,\dots,c_{k-1}$ einer reellen Funktion fct mit Hilfe des FFT-Algorithmus berechnet $ (c_{-k}=\bar{c_k})\,.$ Beginnen Sie mit $ n=2^\ell \geq 2k$ Funktionsauswertungen und verdoppeln Sie die Punktzahl $ n$, bis alle Approximationen sich um höchstens tol unterscheiden.
Die Matlab Funktion c=fft(x) liefert die Koeffizienten

$\displaystyle c_k=\sum\limits_{\ell=0}^{n-1} x_\ell \exp(-2\pi\mathrm{i}k\ell/n)\,. $

Die Riemann-Summe für $ 2n$ Intervalle für das Integral der Fourier-Koeffizienten ist

$\displaystyle c_k \approx \dfrac{1}{2n} \sum\limits_{\ell=0}^{2n-1} f(2\pi \ell/(2n)) \exp(-2\pi\mathrm{i}k \ell/(2n))\,. $

Diese kann in zwei Teilsummen aufgeteilt werden

$\displaystyle c_k$ $\displaystyle \approx$ $\displaystyle \dfrac{1}{2}\dfrac{1}{n} \sum\limits_{\ell=0}^{n-1} f(2\pi 2\ell/... ...its_{\ell=0}^{n-1} f(2\pi (2\ell+1)/(2n)) \exp(-2\pi\mathrm{i}k (2\ell+1)/(2n))$  
  $\displaystyle =$ $\displaystyle \dfrac{1}{2}\dfrac{1}{n} \sum\limits_{\ell=0}^{n-1} f(2\pi \ell/n... ...pi \ell/n +\pi\ell/n)) \exp(-2\pi\mathrm{i}k \ell/n) \exp(-\pi\mathrm{i}k/n)\,.$  

Der erste Summand ist 1/2 mal die Riemansumme für die halbe Intervallzahl, der zweite ist 1/2 mal die Riemann-Summe für die Intervall-Mittelpunkte bei der halben Intervallzahl mit dem Faktor $ \exp(-\pi\mathrm{i}k/n)$ korrigiert.

Unter Vermeidung von doppelten Funktionsauswertung ergibt sich also die Funktion 


function [c, steps] = fourier_coefficients(fct,k,tol)

% function c = f_coef(fct,k,tol)
% Fourier coefficients c_0, ..., c_k-1 of fct 
% error <= tol, steps: number of refinements 

MAX_STEPS = 8; 

n = nextpow2(k); % for oiptimal fft
x = [0:n-1]*2*pi/n;

% initial Riemann sum
c = fft(feval(fct,x))/n; 
c = c(1:k); % n >= k; only first k coefficients are relevant

for steps = 1:MAX_STEPS
   % Riemann sum at interval mid-points
   d = fft(feval(fct,x+pi/n))/n;
   % combine Riemann sums as in FFT-recursion
   d = d(1:k) .* exp(-pi*i*[0:k-1]/n);
   c = (c+d)/2;
   % check tolerance: c-c_old = c - (2c -d) = d-c
   if max(abs(d-c))<tol; 
      return; 
   end;
   % refinement, double points
   n = 2*n; 
   x = [0:n-1]*2*pi/n;
end;
[Aufgabe]
  automatisch erstellt am 27. 11. 2007