# pilotpirx/numerik1

### Subversion checkout URL

You can clone with
or
.

Aufgabe 3 fertig

commit d33e49dddcfa6588e809a56b2e66da4e91a434d8 1 parent 3f8efda
aseyboldt authored
Showing with 280 additions and 1 deletion.
1. +79 −0 aufgabe3.py
2. +58 −0 fft.c
3. +137 −0 quickvis.py
4. +6 −1 wscript
79 aufgabe3.py
 @@ -0,0 +1,79 @@ +import ctypes +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +from quickvis import Interact +import guidata.dataset.dataitems as di +from math import floor, log + +_ndpointer_complex = np.ctypeslib.ndpointer(dtype=np.complex, + ndim=1, + flags='C_CONTIGUOUS') + +_libfft = np.ctypeslib.load_library('libfft', 'build') +_libfft.fft.argtypes = [_ndpointer_complex, + ctypes.c_int] + +def fft(array): + n = len(array) + if n == 0: + raise ValueError("length of array is 0") + while n > 1: + if n % 2 != 0: + raise ValueError("length of array is not a power of 2") + n /= 2 + points = np.array(array, dtype=np.complex, copy=True) + _libfft.fft(points, len(points)) + return points + +def ifft(array): + return fft(array.conj()).conj() + +saw = lambda x: x.copy() +rect = lambda x: np.where(x > 0, np.pi, - np.pi) +def triag(x): + y = np.zeros_like(x) + y = np.where((x <= np.pi / 2) & (x >= - np.pi / 2), - 2 * x, y) + y = np.where(x > np.pi / 2, + 2 * x - 2 * np.pi, y) + y = np.where(x < - np.pi / 2, + 2 * x + 2 * np.pi, y) + return y + + +@Interact +def show_fft(fig, + func=di.ChoiceItem("Funktion", + [(saw, "saw"), (rect, "rect"), (triag, "triag")]), + n=di.IntItem("N", default=32)): + + fig.clear() + + #x = np.linspace(- np.pi, np.pi, n) + x = 2 * np.pi * np.arange(n) / n - np.pi + y = func(x) + + gs = gridspec.GridSpec(3, 2) + + ax = fig.add_subplot(gs[0, 0:2]) + ax.plot(x, y.real, label=r'\$\Re(f)\$') + ax.legend(loc='upper left') + + ax = fig.add_subplot(gs[1, 0]) + ax.plot(x, fft(y).real, label=r'\$\Re(F(f))\$') + ax.legend() + + ax = fig.add_subplot(gs[1, 1], sharex=ax) + ax.plot(x, fft(y).imag, label=r'\$\Im(F(f))\$') + ax.legend() + + ax = fig.add_subplot(gs[2, 0]) + ax.plot(x, ifft(fft(y)).real, label=r'\$\Re\tilde FF(f)\$') + ax.legend(loc='upper left') + + ax = fig.add_subplot(gs[2, 1], sharex=ax) + ax.plot(x, ifft(fft(y)).imag, label=r'\$\Im\tilde FF(f)\$') + ax.legend(loc='upper left') + + fig.tight_layout() + +if __name__ == '__main__': + show_fft.show_gui()
58 fft.c
 @@ -0,0 +1,58 @@ +#include +#include +#include +#include +#include + +#define PI 3.1415926535897932384626433832795028841971 + +void fft_recurse(gsl_complex *data, int length, int stride, + gsl_complex *out) +{ + int k; + gsl_complex *inner_out; + + if (length == 1) { + out[0] = data[0]; + return; + } + + inner_out = malloc(length * sizeof(gsl_complex)); + + if (!inner_out) { + exit(1); + } + + fft_recurse(data, length / 2, stride * 2, inner_out); + fft_recurse(data + stride, length / 2, stride * 2, inner_out + length / 2); + + for (k = 0; k < length / 2; k++) { + out[k] = + gsl_complex_add( + inner_out[k], + gsl_complex_mul( + inner_out[k + length / 2], + gsl_complex_polar(1, - 2 * PI * k / length))); + out[k + length / 2] = + gsl_complex_sub( + inner_out[k], + gsl_complex_mul( + inner_out[k + length / 2], + gsl_complex_polar(1, - 2 * PI * k / length))); + } + + free(inner_out); +} + + +void fft(gsl_complex *data, int length) +{ + int k; + double sqrt_N = sqrt(length); + gsl_complex *out = malloc(length * sizeof(gsl_complex)); + fft_recurse(data, length, 1, out); + for (k = 0; k < length; k++) { + data[k] = gsl_complex_div_real(out[k], sqrt_N); + } + free(out); +}
137 quickvis.py
 @@ -0,0 +1,137 @@ +import inspect, sys +import guidata +import guidata.dataset.dataitems as di +import guidata.dataset.datatypes as dt +from guidata.dataset.qtwidgets import DataSetEditGroupBox +from matplotlib.backends.backend_qt4agg import FigureCanvasQTAgg as FigureCanvas +from matplotlib.backends.backend_qt4agg import NavigationToolbar2QTAgg as NavigationToolbar +from matplotlib.pyplot import Figure +from PyQt4 import QtGui, QtCore + +class Interact(QtCore.QObject): + def __init__(self, func): + QtCore.QObject.__init__(self) + self._func = func + + class Form(dt.DataSet): + pass + + self._form = Form + self._process_args() + self._figure = Figure() + self._widget = None + + self._app = guidata.qapplication() + + def _process_args(self): + args, varargs, keywords, defaults = inspect.getargspec(self._func) + if varargs is not None or keywords is not None: + raise ValueError("Illegal arguments") + + single_args = args[:len(args) - len(defaults)] + default_args = args[len(args) - len(defaults):] + + for i, name_default in enumerate(zip(default_args, defaults)): + name, default = name_default + if not isinstance(default, di.DataItem): + raise ValueError("Illegal arguments, not a DataItem: " + name) + self._form._items.append(default) + self._form._items[-1].set_name(name) + self._form._items[-1]._order = i + 1 + + if len(single_args) > 1: + raise ValueError("Illegal arguments") + + def edit_data(self): + self._form.edit() + return self._form + + def _get_qtwidget(self, parent=None): + self._widget = DataSetEditGroupBox(None, self._form) + self._widget.setParent(parent) + return self._widget + + def show_gui(self): + main_window = ApplicationWindow(self) + main_window.show() + self._app.exec_() + + def _apply_func(self): + assert self._widget is not None + args = {} + for widget in self._widget.edit.widgets: + args[widget.item.item._name] = widget.item.get() + try: + self._func(self._figure, **args) + self.emit(QtCore.SIGNAL("figure_changed()")) + except Exception as e: + import traceback + traceback.print_exc() + QtGui.QMessageBox.warning(None, type(e).__name__, e.message) + + + + def __call__(self, *args, **kwargs): + return self._func(*args, **kwargs) + + +class MplCanvas(FigureCanvas): + """Ultimately, this is a QWidget (as well as a FigureCanvasAgg, etc.).""" + def __init__(self, parent=None, figure=None): + + FigureCanvas.__init__(self, figure) + self.setParent(parent) + + FigureCanvas.setSizePolicy(self, + QtGui.QSizePolicy.Expanding, + QtGui.QSizePolicy.Expanding) + FigureCanvas.updateGeometry(self) + + +class ApplicationWindow(QtGui.QMainWindow): + def __init__(self, interact): + QtGui.QMainWindow.__init__(self) + self.setAttribute(QtCore.Qt.WA_DeleteOnClose) + self.setWindowTitle("application main window") + + data_widget = interact._get_qtwidget(parent=self) + mpl_widget = MplCanvas(parent=self, figure=interact._figure) + + l = QtGui.QVBoxLayout() + l.addWidget(data_widget) + l.addStretch() + + data_widget_top = QtGui.QWidget() + data_widget_top.setLayout(l) + + data_widget_dock = QtGui.QDockWidget() + data_widget_dock.setWidget(data_widget_top) + + toolbar = NavigationToolbar(mpl_widget, self) + + self.addDockWidget(QtCore.Qt.LeftDockWidgetArea, data_widget_dock) + self.addToolBar(QtCore.Qt.TopToolBarArea, toolbar) + + self.connect(data_widget, QtCore.SIGNAL("apply_button_clicked()"), + interact._apply_func) + self.connect(interact, QtCore.SIGNAL("figure_changed()"), + mpl_widget.draw) + + self.setCentralWidget(mpl_widget) + + + def fileQuit(self): + self.close() + + def closeEvent(self, ce): + self.fileQuit() + +def beispiel(): + @Interact + def func(fig, a=di.ChoiceItem("Faktor", ["a", "b"])): + import numpy as np + ax = fig.add_subplot(111) + ax.plot(np.linspace(0, 1), a * np.linspace(0, 1)**2) + + func.show_gui() +
7 wscript
 @@ -6,9 +6,14 @@ def options(opt): def configure(ctx): ctx.load('compiler_c python') - ctx.env.append_unique('CFLAGS', ['-O3']) + ctx.env.append_unique('CFLAGS', ['-O3', '-Wall']) ctx.check_python_module('numpy') ctx.check_python_module('matplotlib') + ctx.check_cc(lib='gsl', uselib_store='gsl') + ctx.check_cc(lib='gslcblas', uselib_store='cblas') + ctx.check_cc(lib='m', uselib_store='m') + ctx.check(header_name='gsl/gsl_complex.h') def build(bld): bld.shlib(source='poly.c', target='poly') + bld.shlib(source='fft.c', target='fft', use=['gsl', 'm', 'cblas'])