Monday, November 11, 2013

Debian 7 "Wheezy"
Personal Installation Guide

Jump to:

1. First steps

2. Install some extremelly important stuff

3. JAVA JRE installation

4. Install fluxbox

5. Install LightDM and/or SLIM greeters

6. Set-up Debian Backports repository

1. First steps

As root, execute the following commands.
root$ apt-get install vim vim-scripts vim-doc vim-common vim-gui-common vim-gtk
Modify /etc/sudoers
foo$ vi /etc/sudoers
Add the following lines at the end of the file.
# Add users to list:
foo ALL=(ALL:ALL) ALL

Now, as normal user install the following useful stuff:
foo$ sudo apt-get install dkms
Also, certify that all mouse and touchpad related software is installed.
foo$ sudo apt-get install xserver-xorg-input-mouse xserver-xorg-input-synaptics tpconfig gpointing-device-settings

Bash-completion: install and enable it.
foo$ sudo apt-get install bash-completion
foo$ sudo vim /etc/bash.bashrc
Uncomment or add the following lines:
# enable bash completion in interactive shells
if ! shopt -oq posix; then
   if [ -f /usr/share/bash-completion/bash_completion ]; then
      . /usr/share/bash-completion/bash_completion
   elif [ -f /etc/bash_completion ]; then
      . /etc/bash_completion
   fi
fi

Updating repository list. First install debian-keyring.
foo$ sudo apt-get install debian-keyring debian-archive-keyring
Most repository paths may be retrieved from https://wiki.debian.org/UnofficialRepositories. Edit /etc/apt/sources.list:
foo$ sudo vim /etc/apt/sources.list
Change/add the following lines (don't forget to add contrib and non-free):
deb http://ftp.pt.debian.org/debian/ wheezy main contrib non-free
deb-src http://ftp.pt.debian.org/debian/ wheezy main contrib non-free

deb http://security.debian.org/ wheezy/updates main
deb-src http://security.debian.org/ wheezy/updates main

# wheezy-updates, previously known as 'volatile'
deb http://ftp.pt.debian.org/debian/ wheezy-updates main contrib non-free
deb-src http://ftp.pt.debian.org/debian/ wheezy-updates main contrib non-free

### 3RD PARTY REPOS
## deb-mutimedia (see http://www.deb-multimedia.org/debian-m)
deb http://ftp.eq.uc.pt/software/unix/Linux/deb-multimedia/ stable main
deb-src http://ftp.eq.uc.pt/software/unix/Linux/deb-multimedia/ stable main

## Google
deb http://dl.google.com/linux/chrome/deb/ stable main
deb http://dl.google.com/linux/talkplugin/deb/ stable main

## Skype
deb http://download.skype.com/linux/repos/debian/ stable non-free

## VirtualBox (see https://www.virtualbox.org/wiki/Linux_Downloads)
deb http://download.virtualbox.org/virtualbox/debian wheezy contrib
Update the repositories cache:
foo$ sudo apt-get update

This may yield several errors related with missing GPG keys, such as:
W: GPG error: http://ftp.eq.uc.pt stable Release: The following signatures couldn't be verified because the public key is not available: NO_PUBKEY 07DC563D1F41B907
Add each of the missing public keys. Exemplifying for the deb-mutimedia repository, we know that the key for this repository is 07DC563D1F41B907, as stated in the error. Thus,
foo$ gpg --keyserver hkp://subkeys.pgp.net --recv-keys 07DC563D1F41B907
foo$ gpg --export --armor 07DC563D1F41B907 | sudo apt-key add -

Finally, update and upgrade the system.
foo$ sudo apt-get update
foo$ sudo apt-get upgrade
foo$ sudo apt-get dist-upgrade

2. Install some extremelly important stuff

Tools I need in my system.
foo$ sudo apt-get install rsync curl lftp wget wput colordiff tkdiff meld
foo$ sudo apt-get install pcmanfm roxterm uzbl hnb
foo$ sudo apt-get install redshift conky xdotool libxdo2 x11-utils

OpenSSH
foo$ sudo apt-get install openssh-server openssh-client

Control version software
foo$ sudo apt-get install git git-gui git-cvs git-svn git-email mercurial subversion subversion-tools bzr bzrtools cvs

Packaging
foo$ sudo apt-get install build-essential debhelper devscripts cdbs dh-make diffutils patch gnupg fakeroot lintian devscripts pbuilder dpatch dput quilt

Printer support
foo$ sudo apt-get install cups cups-pdf system-config-printer printer-driver-hpijs

Apache web server
foo$ sudo apt-get install apache2-mpm-itk apache2-mpm-event apache2-mpm-prefork apache2-mpm-worker mysql-server php-pear php5 php5-gd php5-mysql php5-imagick php5-curl phpmyadmin cronolog

Libre Office
foo$ sudo apt-get install libreoffice libreoffice-gtk

Mplayer
foo$ sudo apt-get install mplayer mencoder mplayer2

VLC
foo$ sudo apt-get install vlc vlc-plugin-sdl vlc-plugin-jack

FFMPEG
foo$ sudo apt-get install ffmpeg x264 libav-tools

Transcode and avidemux
foo$ sudo apt-get install transcode transcode-utils
foo$ sudo apt-get install avidemux avidemux-cli avidemux-common avidemux-plugins

Install dvd stuff (including libdvdcss)
foo$ sudo apt-get install libdvdnav4 libdvdread4 libdvdcss2 regionset

CD and DVD cli burners
foo$ sudo apt-get install wodim genisoimage growisofs

Non-free codecs
foo$ sudo apt-get install gstreamer0.10-plugins-ugly
foo$ wget http://www.deb-multimedia.org/pool/non-free/w/w64codecs/w64codecs_20071007-dmo2_amd64.deb
foo$ sudo dpkg -i w64codecs_20071007-dmo2_amd64.deb
foo$ rm w64codecs_20071007-dmo2_amd64.deb

MS core fonts (with following debian repos: contrib non-free)
foo$ sudo apt-get install ttf-mscorefonts-installer

Flash player (with following debian repos: contrib non-free)
foo$ sudo apt-get install flashplugin-nonfree

3. JAVA JRE installation

First make sure all openjdk java versions are installed.
foo$ sudo apt-get install openjdk-6-jre oracle-java7-installer icedtea6-plugin icedtea-6-plugin

Java package (https://wiki.debian.org/JavaPackage)
foo$ sudo apt-get install java-package

To install JAVA 7 from Oracle, choose one of the following (either 3.a, 3.b or 3.c)

3.a As in https://wiki.debian.org/JavaPackage

Check both wiki.debian.org/JavaPackage and https://d.stavrovski.net/blog/installing-oracle-java7-on-debian-wheezy. Find the latest version number in www.oracle.com/technetwork/java/javase/downloads. In the present case this was Java Platform (JDK) 7u45, thus the file to download was jdk-7u45-linux-x64.tar.gz. Download JAVA 7 from Oracle’s website using wget:
foo$ wget --no-cookies \
   --header "Cookie: gpw_e24=http%3A%2F%2Fwww.oracle.com" \
   "http://download.oracle.com/otn-pub/java/jdk/7/jdk-7u45-linux-x64.tar.gz" \
   -O /tmp/jdk-7u45-linux-x64.tar.gz --no-check-certificate
foo$ make-jpkg /tmp/jdk-7u45-linux-x64.tar.gz
foo$ dpkg -i oracle-j2sdk1.7_1.7.0+update45_amd64.deb

Optionally, for 32 bit systems:
foo$ wget --no-cookies \
   --header "Cookie: gpw_e24=http%3A%2F%2Fwww.oracle.com" \
   "http://download.oracle.com/otn-pub/java/jdk/7u21-b11/jdk-7u45-linux-i586.tar.gz" \
   -O /tmp/jdk-7u45-linux-i586.tar.gz --no-check-certificate
foo$ make-jpkg /tmp/jdk-7u45-linux-i686.tar.gz
foo$ dpkg -i oracle-j2sdk1.7_1.7.0+update45_i686.deb

Debian alternatives system will have it already set-up automatically. Use it to choose what java version to use:
foo$ update-alternatives --display java
foo$ sudo update-alternatives --config java

NOTE: if there are any problem with update-alternatives, configure manually:
foo$ JHome=/opt/java-oracle/jdk1.7.0
foo$ update-alternatives --install /usr/bin/java java ${JHome%*/}/bin/java 20000
foo$ update-alternatives --install /usr/bin/javac javac ${JHome%*/}/bin/javac 20000

3.b Using the one from webupd8team (same as for Ubuntu)

Referring to www.webupd8.org/2012/06/how-to-install-oracle-java-7-in-debian.html:
foo$ echo "deb http://ppa.launchpad.net/webupd8team/java/ubuntu precise main" | sudo tee -a /etc/apt/sources.list.d/java-oracle-webupd8team.list
foo$ echo "deb-src http://ppa.launchpad.net/webupd8team/java/ubuntu precise main" | sudo tee -a /etc/apt/sources.list.d/java-oracle-webupd8team.list
foo$ sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys EEA14886
foo$ sudo apt-get update
foo$ sudo apt-get install oracle-java7-installer

3.c Using the one from duinsoft

Referring to www.fandigital.com/2012/07/install-latest-java-in-debian-ubuntu.html :
foo$ echo "deb http://www.duinsoft.nl/pkg debs all" | sudo tee -a /etc/apt/sources.list.d/java-oracle-duinsoft.list
foo$ sudo apt-key adv --keyserver keys.gnupg.net --recv-keys 5CB26B26
foo$ sudo apt-get update
foo$ sudo apt-get install oracle-java7-installer

4. Install fluxbox

Install fluxbox:
foo$ sudo apt-get install fluxbox

Configure the x-window-manager in update-alternatives:
foo$ sudo update-alternatives --config x-window-manager
  Selection    Path                   Priority   Status
------------------------------------------------------------
* 0            /usr/bin/openbox        90        auto mode
  1            /usr/bin/openbox        90        manual mode
  2            /usr/bin/startfluxbox   50        manual mode

Enter to keep the current selection[*], or type selection number: 

Choose the startfluxbox (number 2). Same for x-session-manager in update-alternatives.
foo$ sudo update-alternatives --config x-session-manager
  Selection    Path                      Priority   Status
------------------------------------------------------------
  0            /usr/bin/startlxde         50        auto mode
  1            /usr/bin/lxsession         49        manual mode
  2            /usr/bin/openbox-session   40        manual mode
* 3            /usr/bin/startfluxbox      30        manual mode
  4            /usr/bin/startlxde         50        manual mode

Enter to keep the current selection[*], or type selection number: 

Choose the startfluxbox (number 3).
NOTE: If any of the startfluxbox pointer is missing, add them:
For x-window-manager:
foo$ sudo update-alternatives --install /usr/bin/x-window-manager \
    x-window-manager /usr/bin/startfluxbox 30 \
    --slave /usr/share/man/man1/x-window-manager.1.gz \
    x-window-manager.1.gz /usr/share/man/man1/fluxbox.1.gz

For x-session-manager:
foo$ sudo update-alternatives --install /usr/bin/x-session-manager \
    x-session-manager /usr/bin/startfluxbox 30 \
    --slave /usr/share/man/man1/x-session-manager.1.gz \
    x-session-manager.1.gz /usr/share/man/man1/fluxbox.1.gz

As normal user add startfluxbox to ~/.xinitrc.bak (in case one later needs to use it, by copying it to ~/.xinitrc).
foo$ vi ~/.xinitrc.bak

#!/bin/bash
###############
#
# ~/.xinitrc
# Executed by startx (run your window manager from here)
#
###############

if [ -d /etc/X11/xinit/xinitrc.d ]; then
        for f in /etc/X11/xinit/xinitrc.d/*; do
                [ -x "$f" ] && . "$f"
        done
        unset f
fi

## SET LOGIN MANAGER
exec startfluxbox

## SET MOUSE
xset m 2 10 &

## RESTORE AUDIO LEVELS
restore_alsa() {
 while [ -z "`pidof pulseaudio`" ]; do
  sleep 0.5
 done
 alsactl -f /var/lib/alsa/asound.state restore
}
restore_alsa &

5. Install LightDM and/or SLIM greeters

5.a LightDM

Install LightDM
foo$ sudo apt-get install lightdm lightdm-gtk-greeter

Change the currend default Display Manager to LightDM
foo$ sudo dpkg-reconfigure lightdm

For changes the the configuration of lightdm:
foo$ sudo vi /etc/lightdm/lightdm.conf

To change the greeter's background (/etc/lightdm/lightdm-gtk-greeter.conf):
foo$ sudo vi /etc/lightdm/lightdm-gtk-greeter.conf

Find the following line:
[greeter]
background=/usr/share/images/desktop-base/login-background.svg

If it points to /usr/share/images/desktop-base/desktop-background use update-alternatives instead.
foo$ sudo update-alternatives --config desktop-background

5.b SLIM

foo$ sudo apt-get install slim

Change the currend default Display Manager to SLIM
foo$ sudo dpkg-reconfigure slim

For changes the the configuration of SLIM:
foo$ sudo vi /etc/slim.conf

Set hidecursor to false and change current theme if needed. Additionally, for fluxbox users, add startfluxbox to the sessions line.
hidecursor      false

current_theme   debian-joy

sessions        /usr/bin/startfluxbox,/usr/bin/openbox-session,/usr/bin/startlxde

To change the SLIM background, change the current_theme background.png file. For example, with current_theme set as debian-joy:
foo$ ls -la /usr/share/slim/themes/debian-joy/
foo$ ln -sf new_background.png /usr/share/slim/themes/debian-joy/background.png

6. Set-up Debian Backports repository

Quoting from https://wiki.debian.org/Backports:
Backports are recompiled packages from testing (mostly) and unstable (in a few cases only, e.g. security updates), so they will run without new libraries (wherever it is possible) on a stable Debian distribution. It is recommended to pick out single backports which fit your needs, and not to use all backports available.
Add Backports repository. Edit /etc/apt/sources.list:
foo$ sudo vim /etc/apt/sources.list
Change/add the following lines:
deb http://ftp.pt.debian.org/debian/ wheezy main contrib non-free
deb-src http://ftp.pt.debian.org/debian/ wheezy main contrib non-free

deb http://security.debian.org/ wheezy/updates main
deb-src http://security.debian.org/ wheezy/updates main

# wheezy-updates, previously known as 'volatile'
deb http://ftp.pt.debian.org/debian/ wheezy-updates main contrib non-free
deb-src http://ftp.pt.debian.org/debian/ wheezy-updates main contrib non-free

## BACKPORTS REPOSITORY deb http://ftp.debian.org/debian wheezy-backports main contrib non-free
### 3RD PARTY REPOS
## deb-mutimedia (see http://www.deb-multimedia.org/debian-m)
deb http://ftp.eq.uc.pt/software/unix/Linux/deb-multimedia/ stable main
deb-src http://ftp.eq.uc.pt/software/unix/Linux/deb-multimedia/ stable main

## Google
deb http://dl.google.com/linux/chrome/deb/ stable main
deb http://dl.google.com/linux/talkplugin/deb/ stable main

## Skype
deb http://download.skype.com/linux/repos/debian/ stable non-free

## VirtualBox (see https://www.virtualbox.org/wiki/Linux_Downloads)
deb http://download.virtualbox.org/virtualbox/debian wheezy contrib
Update the repositories cache:
foo$ sudo apt-get update

As an example, suppose one want's a newer linux kernel. Before adding the Backports repository the list of available kernels was:
foo$ apt-cache search linux-image
linux-headers-3.2.0-4-amd64 - Header files for Linux 3.2.0-4-amd64
linux-headers-3.2.0-4-rt-amd64 - Header files for Linux 3.2.0-4-rt-amd64
linux-image-3.2.0-4-amd64 - Linux 3.2 for 64-bit PCs
linux-image-3.2.0-4-amd64-dbg - Debugging symbols for Linux 3.2.0-4-amd64
linux-image-3.2.0-4-rt-amd64 - Linux 3.2 for 64-bit PCs, PREEMPT_RT
linux-image-3.2.0-4-rt-amd64-dbg - Debugging symbols for Linux 3.2.0-4-rt-amd64
linux-image-2.6-amd64 - Linux for 64-bit PCs (dummy package)
linux-image-amd64 - Linux for 64-bit PCs (meta-package)
linux-image-rt-amd64 - Linux for 64-bit PCs (meta-package), PREEMPT_RT
nvidia-kernel-3.2.0-4-amd64 - NVIDIA binary kernel module for Linux 3.2.0-4-amd64
After adding the Backports repository:
foo$ apt-cache search linux-image
linux-headers-3.2.0-4-amd64 - Header files for Linux 3.2.0-4-amd64
linux-headers-3.2.0-4-rt-amd64 - Header files for Linux 3.2.0-4-rt-amd64
linux-image-3.2.0-4-amd64 - Linux 3.2 for 64-bit PCs
linux-image-3.2.0-4-amd64-dbg - Debugging symbols for Linux 3.2.0-4-amd64
linux-image-3.2.0-4-rt-amd64 - Linux 3.2 for 64-bit PCs, PREEMPT_RT
linux-image-3.2.0-4-rt-amd64-dbg - Debugging symbols for Linux 3.2.0-4-rt-amd64
linux-image-2.6-amd64 - Linux for 64-bit PCs (dummy package)
linux-image-amd64 - Linux for 64-bit PCs (meta-package)
linux-image-rt-amd64 - Linux for 64-bit PCs (meta-package), PREEMPT_RT
nvidia-kernel-3.2.0-4-amd64 - NVIDIA binary kernel module for Linux 3.2.0-4-amd64
linux-headers-3.12-0.bpo.1-amd64 - Header files for Linux 3.12-0.bpo.1-amd64
linux-headers-3.12-0.bpo.1-rt-amd64 - Header files for Linux 3.12-0.bpo.1-rt-amd64
linux-headers-3.13-0.bpo.1-amd64 - Header files for Linux 3.13-0.bpo.1-amd64
linux-image-3.12-0.bpo.1-amd64 - Linux 3.12 for 64-bit PCs
linux-image-3.12-0.bpo.1-amd64-dbg - Debugging symbols for Linux 3.12-0.bpo.1-amd64
linux-image-3.12-0.bpo.1-rt-amd64 - Linux 3.12 for 64-bit PCs, PREEMPT_RT
linux-image-3.12-0.bpo.1-rt-amd64-dbg - Debugging symbols for Linux 3.12-0.bpo.1-rt-amd64
linux-image-3.13-0.bpo.1-amd64 - Linux 3.13 for 64-bit PCs
linux-image-3.13-0.bpo.1-amd64-dbg - Debugging symbols for Linux 3.13-0.bpo.1-amd64
linux-image-amd64-dbg - Debugging symbols for Linux amd64 configuration (meta-package)
linux-image-rt-amd64-dbg - Debugging symbols for Linux rt-amd64 configuration (meta-package)
nvidia-kernel-3.12-0.bpo.1-amd64 - NVIDIA binary kernel module for Linux 3.12-0.bpo.1-amd64
To install the 3.12 linux kernel:
foo$ sudo apt-get install linux-image-3.12

Thursday, February 2, 2012

On CFT, DFT and their normalization

For me to do a FFT is a bit uncommon. In fact I rather prefer to do a DFT, to avoid having spurious lobes and to try to skip the process of windowing and all that. It is always with care that I approach it, as although one has a basic knowledge of its workings, it is kind of an old lore with dust when one does not use it very often. Well the thing is, although one knows a bit of Fourier analysis, to handle an FFT algorithm is sometimes handling a black box and the worse happens when we think we know what we're doing, but in truth we don't.

I'm writing this because I wanted to compare a the DFT and the CFT of a specific function, to which the first apparent difficulty was to normalize the DFT output. It is trivial to find how one should normalize the power spectral densities (PSD) after doing an FFT with MATLAB, OCTAVE or Python. But to normalize already on the DFT output, that is kind of difficult to find. But alas, I knew my function and my CFT, so this turned to be quite simple. The function is an Agnesi shape (or a version of the Cauchy distribution):
$$ f(x)=h\,\left[1+\left(\frac{x}{a}\right)^2\right]^{-1} $$
to which the continuous Fourier transform yields:
$$ F(\xi)=\int\limits_{-\infty}^{\infty} f(x)\,e^{-2\pi\,i\,x\,\xi}\,dt = \pi\,a\,h\,e^{-2\pi\,|\xi|\,a} $$
So it is easy to find that the amplitude at the zero frequency is equivalent to the area below the spatial (or time) domain curve:
$$ F(\xi=0)=\pi\,a\,h\;\equiv\;\int\limits_{-\infty}^{\infty} f(x)\,dx $$
Thus I want the output of my DFT to have a zero coefficient with the same properties, which gives the normalization I must apply.

Advancing to the DFT, one should first define what is exactly the discrete signal. Defining a time array with N equaly spaced bins of size $\Delta t$:
$$ \mathsf{t}=\{0,\,1,\,\ldots,\,N-1\} \times \Delta t \quad\Rightarrow\quad \mathsf{t}_n=n\,\Delta t\,, \;n=0,\,1,\,\ldots,\,N-1 $$
and the respective sampled signal:
$$ \mathsf{f}=\{ f_0,\,f_1,\,\ldots,\,f_{N-1}\} $$
The sampling frequency is the inverse of the sampling period, $f_s={\Delta t}^{-1}$ .

The DFT is thus defined as:
$$ \operatorname{DFT}(\mathsf{f}) = \mathsf{\widehat{F}}_k = \sum\limits_{n=0}^{N-1} \mathsf{f}_n \, e^{-2\pi\,i\,\frac{n}{N}\,k} \;,\quad k = 0,\,1,\,\ldots,\,Nk-1 $$
As I prefer to have the same size for the array returned by the DFT, $Nk$, as the number of samples I have available in my signal, $N$, I sometimes perform a DFT instead of the more faster FFT (when this number is not a power of two).

This is the same definition as the implementation in Numpy, as well as MATLAB and OCTAVE. As it is stated in the Numpy function reference list:
The values in the result follow so-called "standard" order: If A = fft(a, n), then A[0] contains the zero-frequency term (the mean of the signal), which is always purely real for real inputs. Then A[1:n/2] contains the positive-frequency terms, and A[n/2+1:] contains the negative-frequency terms, in order of decreasingly negative frequency.

By making some simple tests, one finds that the DFT coefficients should be multiplied by the sampling period $\Delta t$ (or divided by the sampling frequency $fs$). This is because the zero frequency coefficient of the DFT is equal to the sum of the elements in the array $\mathsf{f}$. If one wants to be consistent with the results of the CFT (recall that $F(0)$ is equal to the area below function $f(t)$), multiplying the DFT by $\Delta t$ yields a zero coefficient that is also the area beneath the discrete array $\mathsf{f}$, using the most basic numeric integration method:
$$ \mathsf{\widehat{F}}_0 = \sum\limits_{n=0}^{N-1} \mathsf{f}_n \quad\Rightarrow\quad \mathsf{F}_0 = \mathsf{\widehat{F}}_0 \, \Delta t \quad\Rightarrow\quad \mathsf{F}_n = \mathsf{\widehat{F}}_n \, \Delta t $$
Now for the tricky part: when I did my DFT, I did it on an array where the x coordinate was defined as follows:
$$ \mathsf{x} = \{0,\,1,\,\ldots,\,N-1\} \times \Delta x + \mathrm{xmin} \\ \ldots = \{\mathrm{xmin}, \, \mathrm{xmin}-\Delta x,\,\ldots, \, -\Delta x,\, 0,\, \Delta x,\, \ldots,\, \mathrm{xmax}-\Delta x,\, \mathrm{xmax}\} $$
and my f array as:
$$ \mathsf{f}_n = f(\mathsf{x}_n) $$
This yielded the following results for the DFT (both real and imaginary parts):



where,the wavenumber is simply $k = 2 \pi \xi$.
At first, not realizing the correctness of the DFT output, I used python functions fftshift and ifftshift, such that:


F = fft(ifftshift(f), n=NDFT)


which is the same has having the following signal:



Although this would give a DFT that agreed with the CFT, conceptually it was not right.

After further analysis my mistake was evident: although the CFT is done on both directions, i.e., it is done using an improper integral spanning from minus infinity to infinity, the DFT is done on a finite signal, where the first index is 0 and the last is $N-1$. Thus when performing the DFT, is over a signal shifted by a distance $\mathrm{xmin}$. Analytically this is the same as:
$$ g(x) = f(x-\gamma) \quad\Rightarrow\quad G(\xi) = e^{-2\pi\,i\,\gamma\,\xi} \, F(\xi) $$
When viewing from the discrete perspective, the shifting is the same as if an index $l$ was applied, thus:
$$ \mathsf{g}_n = \mathsf{f}_{n-l} \quad\Rightarrow\quad \mathsf{G}_k = e^{-2 \pi\,i\,k\,\frac{l}{N}} \, \mathsf{F}_{k} $$
where:
$$ l = \begin{cases} \frac{N-1}{2} \,, & \text{if} \; \mod(N, 2) \ne 0 \\ \frac{N}{2}-0.5 \,, & \text{if} \; \mod(N, 2) = 0 \end{cases} $$
As I want to compare the DFT with the CFT performed on the Agnesi shape without any shifting, I will obtain $\mathsf{G}$ and from it compute $\mathsf{F}$, such that:
$$ \mathsf{F}_k = e^{+2\pi\,i\,k\,\frac{l}{N}} \, \mathsf{G}_{k} $$
which in python code, is equivalent to:

F = np.fft.fft(f, n=NDFT)

if 0==(NDFT % 2):
    l = N/2-0.5
else:
    l = (N-1)/2

F = F * np.exp(2*np.pi * 1j * np.arange(0, NDFT) * l/float(N))

With this correction the DFT yields:





Power Spectral Density (PSD)


Regarding the power spectral density, I obtained:


When searching for ways to normalize the PSD, a common way is (same as in the MATLAB User's Guide):


PSD = np.abs(F)**2/f.size/FS


or using the period (spatial resolution in this case):


PSD = np.abs(F)**2 * DX/f.size



Resorting to the definition of the energy spectral density:
$$ \Phi(\xi) = F(\xi) \cdot \text{conj}\left(F(\xi)\right) \equiv |F(\xi)|^2 $$
and the PSD for an analytic continuous signal:
$$ S(\xi) = \lim\limits_{T\rightarrow\infty} \frac{1}{T} \Phi(\xi) = \lim\limits_{T\rightarrow\infty} \frac{ F(\xi) \cdot \text{conj}\left(F(\xi)\right) }{T} $$
Thus, when looking back to the way the DFT was normalized, in order to resemble the CFT output, generically one can write:
$$ \text{DFT}(\mathsf{f}) = \widehat{\mathsf{F}} \quad \Rightarrow\quad |\widehat{\mathsf{F}}|^2 \, \frac{\Delta t}{N} \;\equiv\; \frac{ \; |\widehat{\mathsf{F}} \, \Delta t|^2 }{\Delta t \, N} \;=\; \frac{|\mathsf{F}|^2}{T} $$
In the above notation, the hat was used to represent the unprocessed output from the DFT, while $\mathsf{F}$ without the hat is the DFT normalized by $\Delta t$. If one obtains the DFT and use it directly to get the PSD, as it is done in the MATLAB reference page, it becomes difficult to distinguish why should the PSD be normalized in such a way, i.e., dividing by $f_s\,N$. Normalizing the DFT output it becomes clear that the PSD is obtained by computing the energy spectra and dividing it by the integral time span of the discrete finite signal. It becomes easier to understand how the PSD relates to the energy spectra and the DFT.



Amplitude or Spectral Density (ASD)


When plotting the frequency content of a signal it is common to resort to the amplitude instead of the PSD. This is useful when dealing with a signal containing several sinusoids at different frequencies, such as:
$$ f(t) = a_1 \sin(2 \pi \, \omega_1 \, t) + a_2 \sin(2 \pi \, \omega_2 \, t) \ldots + a_m \sin(2 \pi \, \omega_m \, t) \ldots + a_M \sin(2 \pi \, \omega_M \, t) $$
Then it would be advantageous if, when viewing the spectrum of $f(t)$, besides identifying the frequencies $\omega_m$, one would also see the correct amplitudes $a_m$ on the vertical axis. The amplitude spectra is given by the magnitude of the DFT, being also related with the energy spectra:
$$ \mathsf{A}_k = \frac{ \text{abs}(\, \widehat{\mathsf{F}}_k ) }{N} \equiv \frac{ \text{abs}(\, \mathsf{F}_k ) }{N \, \Delta t} $$
where the division by $N$ is present to normalize the amplitude,
such that the zero coefficient becomes equal to the mean of the array $\mathsf{f}_n$. Using the DFT output normalized by multiplying with $\Delta t$, again one observes that the normalizing quantity of the ASD becomes the integral time span of the signal. Analytically, the ASD is related to both the energy and power spectra by:
$$ A(\xi) = \frac{ |F(\xi)| }{N \, \Delta t} = \frac{\sqrt{\Phi(\xi)}}{T} \equiv \sqrt{ \frac{ \mathrm{S}_f(\xi) }{T} } $$
Now, the DFT is composed by the zero frequency, a set of floor(NDFT/2) coefficients for positive frequencies an equal set of floor(NDFT/2) coefficients for negative frequencies. The amplitude $a_m$ is thus returned in the ASD, but half at the positive frequency $\omega_m$, while the other half is in the negative frequency $-\omega_m$. AS one is only interested in plotting the positive frequencies, to have the correct value of $a_m$ when plotting the ASD:
$$ \widetilde{\mathsf{A}}_k = 2\,\mathsf{A}_k \;,\quad \text{for}\quad k>0 \\ \widetilde{\mathsf{A}}_0 = \mathsf{A}_0 $$
To finish this post, an example for a signal composed of two sinusoids:



This figure was generated with the following python code:



Thursday, December 15, 2011

2D root finding algorithm using a recursive bisection method for convex quadrilaterals wich allows to find more than 1 zero

Lately I've been occupied in finding critical points in a 2D vector field. These are locations on a X, Y domain where both components of a vector field go to zero, i.e., (U, V) = (0,0). This is the same has having two functions, f(x,y) and g(x,y), and I want to find the set positions (xp,yp) where simultaneously f(xp,yp)=0 and g(xp,yp)=0.
This is not a topic I'm studying and I end up finding these points by mapping the domain with streamlines and verifying for patterns which indicate that a critical point is present (either saddles, nodes or focus). By verifying I meant visual inspection, or creating streamlines and looking at those, which is very boring.

Fortunately I've found the work of A. Globus, C. Levit, and T. Lasinski, where for the 2D case they do this process analytically by fitting bilinear interpolations and verifying the several particular cases (no solution inside the quadrilateral, 1 solution, 2 solutions, imaginary roots or infinite solutions).

Throughout this process I thought of a method of refining a root inside a quadrilateral, akin to the bisection method, but different. Instead of dividing to cut my quadrilateral in order to refine my solution, I create a tree where I further divide my element, letting me search for more roots in other locations. The best way to implement this was with a recursive function, as the quadrilateral divisions and root finding tests are the same to the children of the first element and subsequently. The test to see if a quadrilateral should be refined further is very basic: if the U values at the vertexes are all either >0 or <0, and the same happens to the V values at the vertexes, then this quadrilateral does not have any zero. Of course, there are functions where the topology would destroy this basic technique, such as one that has 4 roots inside the quadrilateral, located in a way that U and V is positive at all vertexes.

However for some functions this would allow to find more than one root. The vector field functions, U and V, may be discrete (leading to the bilinear interpolation case) or have an analytic representation. This would be achieved by changing the ufun and vfun.
Unfortunately it is not robust enough to sort out between repeated zeros (for example, if 2 or 4 children of a quadrilateral have the same vertex which coincides with a root, neither to work well under infinite solutions.

Either way it may be a useful piece of code for the future.

Friday, October 21, 2011

Colorbar titles in matplotlib

In matplotlib sometimes you want to put some title in a colorbar, when producing a filled contour plot or similar. One way is to simply:

Another way is to use the set_title or the set_xlabel of your colorbar axes. For the first:

However in this way I don't know how to set the title at a distance from the axes. For the second, this can be done using the labelpad:

Tuesday, October 11, 2011

Avoid possible division by zero by limiting the divisor

The friendly "division by zero" warning has appeared to almost anyone that had a pocket calculator. Although this can be a nice greeting from our little digital friends that is bidding us to give it some more though on our problem and the required computations, it is generally a pain if one has a large code and has to wait several hours, if not days, to finish a computation. It is sometimes natural that in the progress of the numerical steps needed to solve a PDE where the initial and some boundary conditions are dubious, that the solution needs to wander a bit before giving up and returning those anxiously awaited results.

So, for an operation "a/b" where "b" is may be 0 at some point, a limiter can be applied such that a minimum threshold, dictated by parameter "SMALL", forces the value to be higher than zero and also to avoid underflows.

Welcome

I am human, thus I forget things. Well, sometimes these are quite useful things that in the instances when these are really needed, coincidently it is when these are harder to remember. I have used notebooks for a long time but when these are also needed, they are never around, which is a beautiful fact adding to the arguments of Murphy's law. As such, tired of searching for these multiple pieces of code when a difficulty comes by, I've decided to create this blog to put the stuff that I learn, from other people or by trial and error, so that it is available through some simple mouse clicks.

My daily activity concerns programming to solve problems related with numerical things.
The main programming I work with is FORTRAN 90 and the interpreted languages I use are Python (mainly numpy, scipy, matplotlib), Matlab, Maxima. My preferred shell is Bash and my text editor is good old Vim. Regarding text processing, Latex is a must.