Creating Bode Plot from Experimental Data

Assuming you have an input signal $ u = A cos(2\pi ft) $ and you measure an output signal
$y = B cos(2\pi ft + \theta) $
The $\frac{B}{A}$ ratio is the gain and $\theta$ is the phase shift for frequency $f$.
You could simply demodulate $y$ by multiplying by $x = cos(2\pi ft) - jsin(2\pi ft)$
$$ z(t) = y*x = \frac{B}{2}(cos(4\pi ft) + cos(\theta) + j(sin(4\pi ft) + sin(\theta))$$
Now simply compute the  average of $z(t)$ over multiple cycles so that the $cos(4\pi ft)$ and $sin(4\pi ft)$ terms will average to 0.
Divide the magnitude of $z_{mean}$ by $A$ to get the gain and compute the angle of $z_{mean}$ to get the phase shift.
The advantage of this method is that by averaging over multiple cycles, you will significantly reduce the noise and improve the accuracy of your phase and gain estimates. Apply this technique for multiple frequencies, you will get your Bode plot.


The (periodic) autocorrelation function of a sinusoid $A\cos(2\pi f_0 t + \theta)$
is

$$\begin{align}
R(\tau) &= \int_0^{f_0^{-1}} A\cos(2\pi f_0 t + \theta)
A\cos(2\pi f_0 (t +\tau)+ \theta)\,\mathrm dt\\
&= \frac{A^2}{2}\int_0^{f_0^{-1}} \cos(2\pi f_0\tau)+\cos(2\pi f_0 (2t +\tau)+ \theta)\,\mathrm dt\\
&= \frac{A^2}{2f_0}\cos(2\pi f_0\tau)
\end{align}$$
which has a peak at $\tau = 0$. Notice in
particular that there is no dependence on $\theta$.   More generally,
if the integration is carried out over an interval $[0,T]$
of arbitrary length, we get
$$\begin{align}
\hat{R}(\tau) &= \int_0^{T} A\cos(2\pi f_0 t + \theta)
A\cos(2\pi f_0 (t +\tau)+ \theta)\,\mathrm dt\\
&= \frac{A^2}{2}\int_0^{T} \cos(2\pi f_0\tau)+\cos(2\pi f_0 (2t +\tau)+ \theta)\,\mathrm dt\\
&= \frac{A^2T}{2}\cos(2\pi f_0\tau) + \left.\frac{A^2}{2}\frac{\sin(2\pi f_0 (2t +\tau)+ \theta)}{4\pi f_0}\right|_0^T\\
&= \frac{A^2T}{2}\cos(2\pi f_0\tau) + \frac{A^2}{2}\frac{\sin(2\pi f_0 (2T +\tau)+ \theta)-\sin(2\pi f_0\tau + \theta)}{4\pi f_0}\\
&= \frac{A^2T}{2} \left[\cos(2\pi f_0\tau) + \frac{\sin(2\pi f_0 (2T +\tau)+ \theta)-\sin(2\pi f_0\tau + \theta)}{4\pi f_0T}\right]
\end{align}$$
which does depend on $\theta$ as well as $T$. Notice that
if $T$ is an integer multiple $kf_0^{-1}$ of the fundamental
period $f_0^{-1}$, then the second term
in brackets disappears entirely and we get
$$\hat{R}(\tau) = kR(\tau).$$
On the other hand, since the numerator in that second term cannot
have value more than $2$, we get that if $T \gg f_0^{-1}$ but
is not an integer multiple of $f_0^{-1}$, the
second term is quite small in magnitude and so the autocorrelation
function that one deduces from $\hat{R}(\tau)$ depends only
very slightly on $\theta$.

With that as prologue, I point out to you that the FFT method of
computing the autocorrelation function of a signal gives the
periodic autocorrelation function of a segment of the
signal, that is, what one gets is effectively an approximation to
$\hat{R}_x(\tau)$ for a periodic signal $x(t)$ that is approximately
$2.5$ periods of a sinusoid repeated over and over again.  In particular,
this periodic repetition is not the same as the sinusoid continuing
smoothly and without break as one would want it to be for the above
results to apply.  Thus, it is not at all surprising that the
results from the
FFT method of computing the autocorrelation function depend more
heavily on the phase of the sinusoid. In approximately
$2.5$ periods, the $sin(\cdot)$ goes from $0$ at $t=0$ to nearly
$0$ again at the end of nearly $2.5$ periods. Thus, the periodic
extension is nearly continuous (though there is an abrupt
change of phase at the transition). On the other hand, with
a $90$-degree phase shift, $\cos(\cdot)$ goes from $1$ at
$t = 0$ to nearly $-1$ at the end of nearly $2.5$ periods,
and so the periodic extension is very discontinuous. This
shows up in the FFT results which don't get a good estimate of
the frequency in the latter case.

If you take the frequency of the sinusoid to be be such
that an integer number of cycles fit into the $N$ that
you are using for your FFT, these artifacts will disappear.


A few conceptual questions about filter, pole, and bilinear

The poles of low pass and high pass Butterworth filters are indeed the same if both filters have the same cut-off frequency. The difference between the two lies in the numerator. A low pass filter has a constant in the numerator, whereas a high pass filter has a constant times $s^N$, where $N$ is the filter order.

The poles of a normalized Butterworth low pass filter (i.e., one with cut-off frequency $\omega_c=1$) lie on a circle with center $s=0$ and radius $1$ in the complex $s$-plane. Of course they are all in the left-half plane:

$$p_k=-e^{j\frac{k\pi}{2N}}\tag{1}\\
k=\pm 1,\pm 3,\ldots,\pm (N-1),\quad N\text{ even}\\
k=0 ,\pm 2,\ldots,\pm (N-1),\quad N\text{ odd}$$

If you have a cut-off frequency $\omega_c\neq 1$, the variable $s$ is transformed as $s\rightarrow \frac{s}{\omega_c}$, so a factor $\frac{1}{s-p_k}$ of the transfer function becomes $\frac{1}{s/\omega_c-p_k}=\frac{1}{\omega_c}\frac{1}{s-\omega_cp_k}$. So the poles are indeed multiplied by the cut-off frequency $\omega_c$, i.e., they lie on a circle with radius $\omega_c$. Their angles of course remain unchanged.

Without pre-warping, the bilinear transformation is usually implemented as

$$s=2f_s\frac{z-1}{z+1}\tag{2}$$

where $f_s$ is the sampling frequency. Expressing $z$ in terms of $s$ gives

$$z=\frac{2f_s+s}{2f_s-s}\tag{3}$$

Equation $(3)$ can be used to determine the pole locations in the complex $z$-plane. If the cut-off frequency of the analog filter is $\omega_c$, its poles are located at $\omega_cp_k$, with $p_k$ given by $(1)$. So the poles in the complex $z$-plane are given by

$$z_{\infty,k}=\frac{2f_s+\omega_cp_k}{2f_s-\omega_cp_k}\tag{4}$$

You can use $(4)$ to check the pole locations you obtained from Matlab (and I don't doubt that they are identical up to numerical inaccuracies).


The point is that the dictionary is redundant or overcomplete, i.e. there is no unique expansion of a given signal, unlike the expansion in an orthonormal basis. This means that one tries to find "the best" expansion of a given signal, e.g. such that the error between the signal and its expansion is small while most expansion coefficients are close to zero, which is advantageous for signal compression. The matching pursuit algorithm is an iterative algorithm trying to find an approximation of the best expansion of a signal given a dictionary. Here you can find another useful overview of the matching pursuit algorithm.


Logarithmic fourier transform(LFT) on audio signal

The simplest and most pragmatic solution is to use a normal FFT of a sufficiently large size that you get the required resolution at the lowest frequency of interest. E.g. if you want 1 Hz resolution at the lowest frequency of interest then you will need a 1 second FFT window, i.e. the FFT size would need to be equal the sample rate, e.g. 44100.

Note that even if you could implement a logarithmic FFT then it would still be bound by the laws of physics (information theory) and you would still need a similar length sample window - all you would gain would be convenience (not having to aggregate output bins) at the expense of performance.


best similarity measure for audio classification

I would like to know which is the best similarity measure that I can adopt to check test signal with training data.

Welcome to engineering. There is no universally "best" approach, as it always depends on your specific models.

Mathematically, as soon as you can define what "similarity" should mean, you might be able to derive something.

Assuming that bird's auditory senses kind of work similar to human hearing, I'd try with the usual audio comparison models, probably a basic MFCC thing.

As a comment: A former advisor of me once said:

  A week in the lab can easily save you up to two hours of reading literature.

Really, there's a whole lot of papers, textbook introductions and wikipedia articles on sound comparison. Start with anything, and read forward. That's going to be faster than asking for a single answer to your acute problem without getting the bigger picture of how things work.


How to make low pass filter using frequency sampling method?

The article is assuming that you already somewhat understand sampling.  As soon as you sample a signal, you need to express all of the frequencies with respect to the sampling rate (or you have to carry the sampling rate inside all of your calculations, which is both tedious and loses some generality).

For discrete-time Fourier transform the definition of the Fourier transform is that you do the integration from $-\pi$ to $\pi$.  The way you integrate a function that's defined as the second one is from calculus class: do the integration only over the non-zero parts of $H_d(\omega)$:

$$h_d(n) = \int_{-\omega_c}^{\omega_c} e^{j n \omega} d\omega$$

You'll find that the result has an infinitely long response, which is a bit awkward for fitting into our real, finite world.  A lot of practical DSP design is doing things like finding an acceptable tradeoff between filter length and how close it comes to what you want.

In your case, you want to cut off at $1000\mathrm{Hz}$ in real-world units, which means you need to scale your cutoff frequency by the sampling rate, so that $\omega_c$ is in radians/sample.  I.e., $\omega_c = \frac{1000\mathrm{Hz}}{F_s}$, where $F_s$ is your sampling rate.


Fourier transform in Matlab and hermitian symmetry

In the 2D (continuous) Fourier domain, the Hermitian symmetry (for functions) writes (conjugate symmetric with respect to the origin):

$$F(-u,-v) = \overline{F}(u,v)\,.$$

For matrix $A_{i,j}$, we should have:

$$ a_{i,j} = \overline{a}_{j,i}$$

which means that the diagonal is real. As you can see, Hermitian matrices and Hermitian functions are slightly different. There is something in the diagonal that needs a center. But diagonals of matrices are real.

Plus, the FFT output matrix is not straightforwardly organized like the (continuous) Fourier transform. Friday night here, I am coming back later on FFT symmetries.


The basic trick is to bound the series above and below. Let us do it on one side, for positive indices.

For any $N> 0$, you can write $N=kN_0+r_N$, with $0\le r_N< N_0$. Then if $a_n$ (here $a_n = |x_n|^2$) is positive, $\sum_{n=0}^{N-1} a_n$ is increasing. Now  $kN_0 \le N< (k+1)N_0$, hence you have:

$$ \sum_{n=0}^{kN_0-1} a_n  \le \sum_{n=0}^{N-1} a_n  \le   \sum_{n=0}^{(k+1)N_0-1} a_n   $$
or

$$ k\sum_{n=0}^{N_0-1} a_n  \le \sum_{n=0}^{N-1} a_n  \le   (k+1)\sum_{n=0}^{N_0-1} a_n   $$

since  $a_n$ is  $N_0$-periodic. Now since

$$\frac{1}{(k+1)N_0} \le \frac{1}{N} \le  \frac{1}{kN_0}$$
by a simple point-wise product of the threefold inequality (valid because positive):

$$ \frac{k}{(k+1)N_0} \sum_{n=0}^{N_0-1} a_n  \le \frac{1}{N}\sum_{n=0}^{N-1} a_n  \le    \frac{(k+1)}{kN_0}\sum_{n=0}^{N_0-1} a_n  \,. $$

Since $k/(k+1)\to 1$ and $(k+1)/k\to 1$ as $N\to\infty$, you have your result, mostly. You can split the original series in $\sum_{-N}^{N}a_n = \sum_{0}^{N} a_n+ \sum_{-N}^{0}a_n - a_0$,  take care of the positive and negative indices separately, and leave $a_0/(2N+1)$ tend to $0$.


The Concepts Behind SVD Based Image Processing

You may think on the SVD as a generalization of the Discrete Fourier Transform.
Namely, it is generates an orthogonal basis to represent the data.
The nice thing about it, it generates the basis according to data (Where the Discrete Fourier Transform basis is the same for any data).
Just like the Fourier Spectrum, you have the "Energy" - The eigenvalue.
The i-th eigen value represent the projection of the data on the i-th basis.
In the image above they show you some Eigen Vectors, the basis functions, of the data. As you can see, it is very similar to Fourier Basis which means the data is well represented by the harmonic discrete signals.

The horizontal axis is the sample number of the Eigen Vector - The basis.
The spectrum is the value of the projection. The vectors are the basis to project upon.
Yes. It is the set of the Singular Values.


How to represent ratio between harmonics in audio spectrum

I think what you are referring to is called the Total Harmonic Distortion (THD). THD is the level of the fundamental divided by the level of the harmonics.  You can express this as either as a ratio (X/Y), a percentage (X/Y * 100) or in dB (20*log10(X/Y)). A percentage is probably what most people are familiar with.


Apply Low Pass (Smoothing) Filter on a Set of Matrices and Reject Outliers

The method used was a parametric representation of the matrix $ H $.

For intuition, thing of rotation matrix in $ \mathbb{R}^{2} $, namely $ H \in \mathbb{R}^{2 \times 2} $.
Since we have series of those (In case of video, a series in time) we can define $ H \left( t \right) $.

As a rotation matrix it can be parameterized by a single parameter $ \theta $ which is the rotation value. Since we are dealing with a time series then we have $ \theta \left( t \right) $.

Now, this is just a simple Time Series and regular methods of smoothing and outlier rejection can be applied.

Another idea would be to find the best base to approximate all matrices (Which with respect to that base the can be diagonalized) and then work on the "Singular / EigenValues".


Matching input buffer size to opus frame size

I ended up solving the problem. Because the input buffer was smaller than 1 frame, I was able to write the input to a circular buffer, wait until I have enough data, then do the processing. The result was the written to a second circular buffer, that would be read every block.


The time span of the zero crossings increases after the final RRC filtering (and the symbol sampling locations converge which is the goal for the benefit of zero ISI but the zero crossing increase in the process is to the detriment of timing recovery!). So if you are using a Gardner TED which is sensitive to this, it is better to have TED prior to RRC filtering, as the timing SNR will be higher.

However a synchronizer such as the Mueller and Mueller which operates on symbol decisions (1 sample per symbol) has better performance after the RRC filter.

Below are details related to the use of the Gardner TED, to show the considerations involved:

See my figures below from when I studied this more closely with the Gardner TED, showing also the trade-off to be made with a lower roll-off (alpha value) for better spectrum efficiency but lower timing SNR due to reduced discriminator slope and higher pattern noise. Shown is the "Mean(TED)" which is the timing discriminator for a Gardner TED (higher slope means higher loop gain/sensitivity), and the pattern noise of the detector versus timing offset. Notice that a higher timing SNR is achieved if you do your timing error detection prior to RRC filtering.  This is because the time span of the zero crossings increase after RRC filtering, while the symbol decision locations converge (see figure of waveform before and after RRC).

For comparison, if desired, I can also later include the M&M synchronizer since I had studied that in similar fashion as well, but for now I bottom lined the comparisons I derived at the end of this post.

Gardner TED for QPSK / QAM:

$TED = I_{2n}(I_{2n+1}-I_{2n-1}+Q_{2n}(Q_{2n+1}-Q_{2n-1}) = real[conj(y_n)(y_{2n+1}-y_{2n-1})]$

This shows the frequency spectrum characteristics of the pattern noise from the Garndner Timing Error Detector:

For comparison (and consideration), below is the Mueller & Mueller (M&M) Synchronizer:

Without elaborating on the Mueller & Mueller synchronizer (M&M) in similar detail, I will at least include below the key take-aways I made when comparing the Gardner TED vs the Mueller & Mueller synchronizer:


In image processing, what is the difference or relationship between noise and artifact?

The difference is that noise may obscure features in an image, while artefacts appear to be features but are not.  If the 'problem' is structured, it is probably an artefact, whereas if it is random, it is probably noise (as a generalisation).

A computed tomography example:  noise will make the image look grainy, and make small differences in contrast difficult (or impossible) to identify.  A streak artefact on the other hand has structure to it, and looks like the patient has a region of low density where they actually don't.  Likewise with ring artefacts - they look like the patient has ring shaped structures within them.  Artefacts are generally caused by assumptions made in the development of the reconstruction algorithm (though not always).

It is similar to the difference between noise and interference in some other fields - interference is structured in some way, having some particular source, while noise is less predictable.

EDIT: Just found something in the glossary of "Computed Tomography" by Willi A. Kalender, 2011 (978-3-89578-317-3):

artifact: part of the contents of an image that does not have a counterpart in the physical object being imaged; in CT images an artifact can appear as false structures (e.g. aliasing artifact, beam hardening artifact, motion artifact, partial volume artifact) or as a corruption of the CT numbers.
noise: contributions to a measured signal due to random processes; in principle these contributions are of statistical nature and do not carry any information about the signal.

If that's of any use to you.


I tried two approaches, naively (using only 3 segments). Surely there would be fancier methods out there.

 RANSAC, supposed to be a robust fitting mechanism. It's easy to stop the algorithm after a number of segments. However it may be difficult to enforce continuity between segments--as seems required in your application-- at least with a simple implementation.

As a proof of concept, I created an image from the data points so that I could use the RANSAC engine available in $ImageLines$, the line detection function of Mathematica.

 Fit a piecewise linear model using a general purpose minimizer. It's easy to enforce segments continuity. Interestingly, testing for residuals and other properties may provide enough information to determine automatically the number of segments--I've not tried it though. That's how it looks in Mathematica:


Simple beat detection algorithm for microcontroller

I did one very simple bpm calculation in matlab and maybe it can help you!

source

Steps:

Decompose the signal in six frequency bands (six level) using DWT
Extrac the Envelope from all sub-bands

To make the Envelope Calculation do you need:

Apply Full Wave Rectification
Apply Low Pass Filter
Apply Downsample
Apply Mean Removal

Autocorrelation


Salt and Paper noise are better dealt and reduced using Median Filter.

The properties of the Salt and Pepper noise make it an outlier in almost any patch of the image.
Hence the best way to remove it is using a method robust to outliers.
Linear filters (Low Pass Filters) are basically weighted mean an hence very sensitive to Outliers which makes them less efficient in this case.

The Median is known and simple robust estimator and hence perform very well in this case.

The way to use it is simple, working by patches / windows.
Iterate through the pixels of the image, for each pixel open a Window of its neighborhood controlled by the parameter which is the Window Radius.
For each window replace the value of the center pixel by the median of the window.

MATLAB Code (From StackExchange GitHub):

% Image Salt and Pepper Noise Filtering
% See http://dsp.stackexchange.com/questions/27298

%% General Parameters and Initialization

clear();
close('all');

% set(0, 'DefaultFigureWindowStyle', 'docked');
defaultLooseInset = get(0, 'DefaultAxesLooseInset');
set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

titleFontSize   = 14;
axisFotnSize    = 12;
stringFontSize  = 12;

thinLineWidth   = 2;
normalLineWidth = 3;
thickLineWidth  = 4;

smallSizeData   = 36;
mediumSizeData  = 48;
bigSizeData     = 60;

randomNumberStream = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');
subStreamNumber = 57162;
set(randomNumberStream, 'Substream', subStreamNumber);
RandStream.setGlobalStream(randomNumberStream);

%% Constants
UNIT8_MIN_VALUE = 0;
UNIT8_MAX_VALUE = 255;

%% Loading Data

mInputImage = imread('../RawData/Image0001.png');
mInputImage = mInputImage(:, :, 1); %<! Assuring Single Channel Image

numRows = size(mInputImage, 1);
numCols = size(mInputImage, 2);
numPixels = numRows * numCols;

%% Parameters
winRadius           = 1; %<! Loacl Window Radius
noiseProbability    = 0.05; %<! Probability of a Pixel to be affected by noise

%% Generating Noisy Image

numNoisyPixels  = round(noiseProbability * numPixels);

vNoisyPixelsIdx         = randperm(numPixels, numNoisyPixels);
vPepperNoisyPixelsIdx   = vNoisyPixelsIdx(1:floor(numNoisyPixels / 2));
vSaltNoisyPixelsIdx     = vNoisyPixelsIdx((floor(numNoisyPixels / 2) + 1):numNoisyPixels);

mNoisyImage                         = mInputImage;
mNoisyImage(vPepperNoisyPixelsIdx)  = UNIT8_MAX_VALUE;
mNoisyImage(vSaltNoisyPixelsIdx)    = UNIT8_MIN_VALUE;

%% Image Median Filtering

mFilteredImage = mNoisyImage;

for iCol = 1:numCols
    for jRow = 1:numRows
        winFirstRowIdx  = max(1, (jRow - winRadius));
        winLastRowIdx   = min(numRows, (jRow + winRadius));
        winFirstColIdx  = max(1, (iCol - winRadius));
        winLastColIdx   = min(numCols, (iCol + winRadius));

        vRowsIdx = [winFirstRowIdx:winLastRowIdx];
        vColsIdx = [winFirstColIdx:winLastColIdx];

        mLocalWin = mFilteredImage(vRowsIdx, vColsIdx);

        mFilteredImage(jRow, iCol) = median(mLocalWin(:));
    end
end

%% Display Results

hFigure = figure('Position', [100, 100, 900, 550], 'Units', 'pixels');
set(hFigure, 'Colormap', gray(256));
hAxes   = axes();
set(hAxes, 'Units', 'pixels')
set(hAxes, 'Position', [50, 50, numCols, numRows]);
hImageObject = image(mInputImage);
set(get(hAxes, 'Title'), 'String', ['Input Image'], ...
    'FontSize', titleFontSize);

hFigure = figure('Position', [100, 100, 900, 550], 'Units', 'pixels');
set(hFigure, 'Colormap', gray(256));
hAxes   = axes();
set(hAxes, 'Units', 'pixels')
set(hAxes, 'Position', [50, 50, numCols, numRows]);
hImageObject = image(mNoisyImage);
set(get(hAxes, 'Title'), 'String', ['Noisy Image'], ...
    'FontSize', titleFontSize);

hFigure = figure('Position', [100, 100, 900, 550], 'Units', 'pixels');
set(hFigure, 'Colormap', gray(256));
hAxes   = axes();
set(hAxes, 'Units', 'pixels')
set(hAxes, 'Position', [50, 50, numCols, numRows]);
hImageObject = image(mFilteredImage);
set(get(hAxes, 'Title'), 'String', ['Filtered Image'], ...
    'FontSize', titleFontSize);

%% Restore Defaults
set(0, 'DefaultFigureWindowStyle', 'normal');
set(0, 'DefaultAxesLooseInset', defaultLooseInset);

Enjoy...


Calculating the network function of a filter

You need to consider the system's frequency response

$$H(j\omega)=\frac{\kappa}{j\omega+a}\tag{1}$$

Now you know that $$H(0)=1\tag{2}$$ and $$|H(j\omega_0)|=\frac12\tag{3}$$ (with $\omega_0=10^4$). Note that $H(0)$ is real-valued, whereas $H(j\omega_0)$ will generally be complex-valued. From (1) and (2) you immediately get $\kappa=a$. Furthermore, since the filter is stable we know that $a>0$ (i.e. the pole must lie in the left half-plane of the complex $s$-plane). Combining (3) with (1) gives

$$\frac{a}{\sqrt{\omega_0^2+a^2}}=\frac12\tag{4}$$

from which it should be easy for you to compute $a$.

Note that you don't need the phase value to specify the system. What should be checked is if the phase value is the actual phase value of the system that you just computed. In order to do this you need to verify that

$$\arg\{H(j\omega_0)\}=-\frac{\pi}{3}\tag{5}$$

Luckily that's the case, otherwise the assumption (1) would be wrong. I leave the proof of (5) up to you.


Digital Lowpass butterworth filter with cut off 500Hz and sampling rate 1.25MSPS

Solution: Take a lower order butterworth filter or sampling frequency

Reason: As discussed here, a high order butterworth filter with a low (relative) cutoff frequency may be numerically unstable due to quantisation noise.

Explanation with figures:

I changed your frequency plot to include the region of interest:

freqz(b,a, logspace(1, 5, 1000), Fs)
ax = findall(gcf, 'Type', 'axes');
set(ax, 'XScale', 'log');

The frequency response is:

which seems very noisy at low frequencies, indicating an unstable filter.

This results in an extremely high output (blue = input, red = output):

If I change the butterworth order to 4, I get:

The output (red) is a shifted version of the input (blue) as can been expected from the frequency response.

Keeping a 6th order butterworth filter, but lowering the sampling frequency also solves your problem.


Is a real Grayscale image always 8 bit per pixel indexed?

I don't think that it is that easy. Grayscale images do have certain benefits when you think of imaging experiments. When the measurement is squeezed into eight bits, even though the measurement is much more precise, then there is a loss of information due to the quantization. See the DICOM format for example that is widely used in medical imaging. It is very common there that the information is encoded in 12 bit, or even more.


Decomposing a DFT into multiple FFT calls

The simplest thing to do if you have access to only a power-of-two FFT is to zero-pad your data going into the $N=2^p$ FFT and interpolate what comes out.  i think that's what MATLAB does for fft() for non-power-of-two DFT.


Laplacian response decays as $\sigma$ increases, and it is the second Gaussian derivative so it is multiplied by $\sigma^2$.

LOG is defined as $\bigtriangledown^2G$ so the scale normalized LOG would be $\sigma^2\bigtriangledown^2G$. You need to get rid of the scaling factor of the Gaussian which is $\sigma^2$.


Multiply a measured quantity of a control-system in the frequency domain during continious operation

I assume that the transfer function that you are using is of a similar form as $V(\eta)$ in the image of your question, with $\eta = \Delta f / f_0$. Often a transfer function is written in terms of the complex Laplace variable $s$, but I assume that your $\eta$ represents some sort of normalized frequency. That is equivalent to using $s = i\,k\,\omega$ with $i$ the unit imaginary number, $k$ some normalization constant and $\omega$ the angular frequency. If you rewrite your transfer function in terms of $s$ instead of $\eta$ you can use one of the many known discretization methods. A simple but good discretization method is the bilinear transform, which requires you to substitute $s$ with the following expression in the transfer function

$$
s = \frac{2}{T} \frac{z - 1}{z + 1}, \tag{1}
$$

with $T$ the sample time of your sensor. After simplifying the resulting (discrete time) transfer function you should have something of the form

$$
V(z) = \frac{b_m\,z^m + b_{m-1}\,z^{m-1} + \cdots + b_0}{z^n + c_{n-1}\,z^{n-1} + \cdots + c_0}, \tag{2}
$$

with $n \geq m$ and all $b_i$ and $c_i$ real scalar constants. When multiplying a discrete signal with the variable $z$ can be interpreted as shifting that signal one sample in time. So given the input signal $a[k]$ (with $k$ an integer used as an index and the time associated with that index could for example be $t = k\,T$) the output $F[k]$ can be related to each other using $(2)$

\begin{align}
(z^n + c_{n-1}\,z^{n-1} + \cdots + c_0)\,F[k] &= (b_m\,z^m + b_{m-1}\,z^{m-1} + \cdots + b_0)\,a[k], \\
F[k+n] + c_{n-1}\,F[k+n-1] + \cdots + c_0\,F[k] &= b_m\,a[k+m] + b_{m-1}\,a[k+m-1] + \cdots + b_0\,a[k].
\end{align}

Solving this for $F[k]$ gives

$$
F[k] = b_m\,a[k+m-n] + b_{m-1}\,a[k+m-1-n] + \cdots + b_0\,a[k-n] - c_{n-1}\,F[k-1] - \cdots - c_0\,F[k-n],
$$

so a weighted sum of previous inputs and outputs. This relation for $F[k]$ allows you to calculate a new output every time you get a new input, so no need to calculate it using chunks. I can be noted that initially you do not know the past values of the input and output. However, usually one set those to zero or some other initial value, which should only change the initial transient behavior of the output (assuming that the filter is stable).


Regarding phase and frequency components in the signal

I find it easier to understand the phase relationships by observing complex tones rather than sines and cosines, which are related using Euler's identity:

$e^{j\theta}= \cos(\theta)+j \sin(\theta)$

$e^{j\theta}$ is simply a phasor with magnitude 1 and angle $\theta$. If that angle is changing with time at a constant rate ($e^{j\omega t})$ the real portion will be a cosine and the imaginary portion will be a sine. If that phasor that is spinning at a constant rate does not start at magnitude = 1 and angle = 0 at time $t = 0$, then we have a phase offset ($e^{j\omega t + \phi})$ where $\phi$ here represents the starting phase at $t = 0$.

Thus a single impulse in the frequency domain is a single spinning phasor in the time domain as shown in the figure below. From that figure the I axis represents the real part of the spinning phasor, which would be a cosine wave, or a sine wave if the starting phase at $t=0$ was $\pi/2$ since $\cos(\theta+\pi/2)= \sin(\theta)$

So to answer your questions with reference to the explanation above:

  Is there any possibility that the phase of the signal remains constant even if the frequency of the signal is changed?

By definition of frequency being phase changing with time, the only signal that has a phase that remains constant is a DC signal. So if you stop the phasor from spinning (at any angle) this would be frequency = 0 (DC).  You can have a constant phase offset for all frequencies, which is the phase you are referring to- this simply means at t=0 every frequency in your scenario has the same phase; but then as time proceeds each frequency would have a phase that would change at a different rate (the slope of the phase vs time IS the frequency!).

  If we observe a real-time signal the signal can be a periodic sine wave signal but not smooth curve as it is accompanied by noise. In such cases how the phase of the signal is varied? Can the noise effect the phase of signal?

Yes of course it does effect the phase, and this is what is called "phase noise". Refer to the plot above of the spinning phasor in the time domain, and now imagine noise being added to it. The noise would appear as a random vector on the end of the phasor (proportional to how strong the noise is), and that vector would change randomly with time. Often this is shown as a round "noise cloud" that is on the end of the vector signifying its statistical distribution with time (but of course at any given instant the noise vector would be at one particular amplitude and phase). This noise cloud (If a complex Gaussian noise) would equally effect the amplitude and phase of the signal, causing both AM (amplitude modulation) and PM (phase modulation), what we would refer to as AM noise and Phase Noise.


Intuition behind the Continuous Wavelet Transform?

I think that the best way to explain the CWT is to start by explaining the Fourier Transform, then move on to explaining the Short-Time Fourier Transform, and then finally explain the CWT as a variation of the STFT.

The Fourier Transform exploits the fact that any decently behaved function can be represented as a sum of sinusoids (i.e. a Fourier series) and that the sinusoid basis posses the property of orthogonality:

$$ \int_{-\infty}^{\infty}  \sin(nx)\,\sin(mx)\,dx = \int_{-\infty}^{\infty}  \cos(nx)\,\cos(mx)\,dx = \begin{cases}
1,  & \text{if $n=m$} \\
0, & \text{if $n\neq m$}
\end{cases}$$

So, since:

$$ e^{iax} = \cos(ax)+i\,\sin(ax)$$

the Fourier Transform is simply doing this integration for all frequencies and keeping track of which outputs are zero (i.e. that frequency is not in the signal) and which are non-zeros (i.e. that frequency is in the data and its output is scaled by how much of it is in there):

$$ S(f) = \int_{-\infty}^{\infty} s(t)\,e^{i2\pi ft}dt $$

In this case you are doing this integration over the entire signal so you can't really tell if the frequency content is changing from the beginning of the signal to the end. One way around this is to compute the Short-Time Fourier Transform: i.e. window the signal, calculate the Fourier transform of the windowed signal, store it, then shift the window down a bit and repeat for all shifts:

$$ S(\tau,f) = \int_{-\infty}^{\infty} w(t- \tau)\,s(t)\,e^{i2\pi ft}dt  $$

where

$$ w(\tau,t)=\begin{cases}
1,  & \text{if $\tau\approx t$} \\
0, & \text{if $\tau \not\approx t$}
\end{cases} $$

The key thing here is that you are calculating the typical Fourier transform but of a new signal that only exists in a localized part of the t-axis. To emphasize this, you can see the new signal whose Fourier transform we are calculating by associating:

$$ S(\tau,f) = \int_{-\infty}^{\infty} [w(t-\tau)\,s(t)]\,e^{i2\pi ft}dt  $$

And here's a graphical example of this showing the new signal for different $\tau$ values and on the right are a few sinusoids to represent what we are using to decompose the signals (i.e. our basis, or kernel).

But, we can also change the association as such without changing the outcome:

$$ S(\tau,f) = \int_{-\infty}^{\infty} s(t)\,[w(t-\tau)\,e^{i2\pi ft}]dt  $$

So this means that instead of windowing our signal, we are windowing our basis functions. But here's the kicker, if we are windowing our basis functions, we don't have to use a constant-size window since we know that a basis function of high frequency will need a shorter window than a basis function of low frequency. This is the whole point of the CWT. It is a decomposition of a signal by "wavelets" (i.e. windowed sinusoids in this case) where the windowing is adaptive to the sinusoid frequency. If the we choose a Gaussian window (as I have chosen in these examples), then our wavelets are called Morlet wavelets (or Gabor wavelets, in some literature).

Finally, you can generalize this for any choice of wavelet that you want. In that generalization, you can describe your wavelet basis functions as being "stretched" and "squeezed" versions of some arbitrary "mother" wavelet, $ \psi $. And so the previous equation can now be written as the final form of the CWT:

$$ S(a,b) = \frac{1}{\sqrt b}\int_{-\infty}^{\infty} s(t)\,\psi (\frac{t-a}{b})dt $$

where $a$ is now what used to be our $\tau$ (i.e. a time shift), and $b$ is called the "scale" which is just a parameter to stretch and squeeze the wavelet (similar to what our parameter $f$ was except now the interpretation is more difficult). And the only reason you have the $\frac{1}{\sqrt b}$ up front is to normalize the wavelets so that they have the same "energy" and you end up comparing apples to apples in your time-frequency representation.

I hope this helps!

-Antonio


This is only a partial answer, but there's a lecture online where Hamming talks about how he came up with his eponymous window. Starting at roughly 15:15 gives the full context.

With a reasonably entertaining story, he credits John Tukey with inventing the theory of windows (for spectrum analysis). However, he introduces the whole subject in the context of using Lanczos sigma factors to reduce Gibbs phenomenon. In addition, in The Art of Doing Science and Engineering (based on the same lectures), he describes how his window is a variation on the Hann window, which he claims was used by von Hann in economics (long before its application in signal processing). That suggests the history goes much further back, depending on how you want to define it.

The book where Tukey first named the Hamming window is The Measurement of Power Spectra from the Point of View of Communications Engineering. Given Hamming's assertion that Tukey invented the theory of windows, it would probably be a good place to start for a deeper understanding of how to design new ones. I think the book is just a reprint of Part I and Part II of his Bell System Technical Journal article, so it's available online.


Maybe this is an appropriate paper for you: Exposing Digital Forgeries in Scientific Images, by Hany Farid (freely available, e.g. at http://www.ists.dartmouth.edu/library/199.pdf).


The PNG format is lossless for RGB24 image data.  However the conversion from YUV to RGB24 is not lossless, as the two formats quantize the color space differently.

To see this, the following animated gif was made by applying your two ffmpeg operations 200 times back-and-forth and collating the resulting 200 images into the gif.

By contrast, the following gif was generated using yuv444p as the YUV format instead of nv12. This format does not do subsampling / as much quantization.


Effect of Sampling frequency on RMS estimation

If the sampling rate meets the nyquist criteria (sampling rate more than twice the highest frequency in the signal) then it won't be a factor in the accuracy of the RMS value.

Since you are measuring velocities the nyquist criteria could be a problem if you don't know the frequencies involved.  The frequencies would be related to the acceleration - higher acceleration means higher frequencies.  I don't know how to tell you to determine the needed sampling rate give the expected acceleration.

If you can't determine the highest frequency (or are limited in your sampling frequency,) then you might run your sampling at the highest available rate, but filter the analog signal with a lowpass with a cutoff frequency lower than half of the sampling rate.  This will avoid the problems that arise when sampling a fast signal with a too slow sampling rate, at the expense of throwing away information if the sampling rate really is too slow.

The number of samples only comes into it when trying to calculate the RMS value of a repeating signal, line a sine wave.  In that case you have to make sure to get a least one full cycle of the lowest frequency.  More could help make up for inaccuracies in the measurement or reduce some of the noise.  In this case, the longer sampling period averages out any measurement errors and noise.

Summary: Make sure that your sampling is fast enough, that the samples are evenly spaced, and that you (obviously) have all the samples for the period of interest.


Is "early vision" chronological or process

I'm pretty sure it means early as in the sorts of things that happen in the first steps of the visual cortex. So things like basic segmentation and shape detection and motion analysis would be early but not higher level reasoning like object recognition.

http://www.cs.duke.edu/~tomasi/papers/tomasi/tomasiEcs00.pdf


This is a saddle point. I would suggest using derivative operators to spot such discontinuities. In fact, Haralick presents a broad overview of these methods in his Topographic Primal Sketch.

The use of approximate Facet model would give you a speed boost. I used them for edge detection here.

There is also one implementation available here.

Even though this question and the topic is rather old, I would like to update that, the method presented in the following work contains no branching and is quite suitable for the shader:

  Birdal, Tolga, Ievgeniia Dobryden, and Slobodan Ilic. X-Tag: A
  Fiducial Tag for Flexible and Accurate Bundle Adjustment. 3D Vision (3DV), 2016 Fourth International Conference on. IEEE, 2016.


What does PSD represent in the spatial domain?

Please see this response. Radiant Flux should take care of the physical units.

The PSD in the spatial domain would represent a "dominant" blob, whose shape is determined by the spectral content of the PSD, underlying the image. (Please note that the PSD is defined over the Fourier Transform of the autocorrelation of a signal).

EDIT:

So, your pixel size is 5mx5m and the actual length of your red line is (let's say) about 500 meters across and contains 164 samples. Let's assume a square pixel size and call it $P_s$, call the length of the line $L$ and the number of samples across the line $L_s$.

You perform an FFT on the profile line ($x_n$) and calculate the amplitude spectrum as the absolute value of the complex spectrum of $x$ (i.e numpy.abs(numpy.fft.fft(myData))). This gives you $|X_k|$.

From that, you can see one sinusoid clearly sticking out of the rest at a frequency bin of 12. What does this mean for the ridge?

The formula that links the bin to the actual spatial frequency is $f = \frac{k}{L_s} \times \frac{1}{P_s}$ where $k$ is the frequency bin.

But, what is $\frac{1}{P_s}$?

The FFT decomposes a waveform into a sum of sinusoids at different frequencies. It will take the waveform from the spatial domain (height  of ridge ($x_n$) at well defined points in space ($n$)) to the frequency domain (amplitude of sinusoids $|X_k|$ that make up the waveform over the same well defined points in space). The variable $k$ ranges from 0 (also known as DC) to $L_s$, corresponding to actual physical frequencies $f$ between DC and $F_s$.

By the way, you only show the real part of the FFT, which means that your $k \in [0 \ldots \frac{L_s}{2}]$ and your $f \in [0 \ldots \frac{F_s}{2}]$.

$F_s$ stands for Sampling Frequency and is usually expressed in Hertz (Hz). 1 Hz means 1 cycle per second. PAUSE.

Imagine a signal in the time domain, it looks like $x(t)$. For every time instance $t$ we collect a sample $x(t)$. If we sample it regularly at some sampling frequency $F_s$ then the distance (the sampling period $T_s$) between successive time instances $t,t+T_s$ is $T_s = \frac{1}{F_s}$. PLAY.

So, you know that you have sampled space once every 5m. That is, $T_s= 5m$. Then, what is $Fs$ when I know $T_s$? ... That's what $\frac{1}{P_s}$ is.

And as you can see from the above simple formulae, the longer the line, the more samples you obtain from the bottom, the more clearly you can tell spatial frequencies apart.

OK, now we know how to translate between what the FFT tells us and actual spatial frequency. For this example $f=\frac{12}{164} \times \frac{1}{5} \approx 0.014Hz$. OK...ok...What does this mean? It means that your ridge, as idealised by a single sinusoid has a period of $T = \frac{1}{0.014} \approx 71.428$ meters. That is, one period of sand raising and then dipping is about 71 meters across.

But, what about height? I read an amplitude of 42 on my single powerful harmonic, does this mean 42 meters?!?

No. Different scientific packages might "implement" slightly different formulas for the FFT. If you are within the Python ecosystem, you are likely to be using numpy. So, for this purpose, take a look at this link. You see how the forward transform is just a sum but the inverse transform has a $\frac{1}{n}$ in front of it? That's because of the way the "averaging" works. So, when you read 42, in reality it is an amplitude of $\frac{42}{L_s} \approx 0.25$ meters or about 25cm. Of course, I am assuming that your height is measured in meters here.

These are the "basics", if you like. There are a few more things that, I suppose, are going to come into play later on.

For example, is a single sinusoid good enough to represent the cross section of a ridge pattern? Personaly, I have no idea, but I am sure that there are ways for you to find out. In general though, you can order the amplitudes of the harmonics of the FFT in decreasing order and keep adding them until you have reached something like 90% of the total amplitude. This might not extend to more than 4-5 coefficients in places where the ridge is very clear. Another way to select coefficients is, of course, by estimating the goodness of fit of a set of $n$ sinusoids to explaining the total variation of the amplitude in the obtained sample, there are many different ways to carry out this step. The more coefficients you include of course, the better your representation but the more difficult the interpretation. However, you can collect these coefficients and correlate them later on with other features from your reef and see if there is something that could explain the ridge pattern better. Maybe stiff corals create strong objects which result in turbulent flow but soft seaweed or flexible flora bends and follows the current and slows it down and contributes to a more laminar flow (?).

Another thing is, how do you judge that what you "see" in front of you is good enough to be characterised as a ridge pattern. The algorithm is unsupervised so it might well go over a patch where the flow is so turbulent that it doesn't leave a well defined pattern. The mathematical pipeline will still do its work and come back with numbers but it is likely that these won't mean anything. You can easily see for example that if we were looking naively for a single powerful sinusoid from [4,4,4,4,5,4,4,4,4,4,4,], we would say 5...But that's not correct. For this reason, I would suggest that as a criterion to whether or not you go ahead with harmonic analysis of a path you first apply a Spectral Flatness metric which will tell you very quickly if the spectrum is "spikey" or close to flat. So, now, you can put a threshold on flatness and say something like "If the spectrum is too flat then don't even bother with FFT and the rest, just characterise this patch of sand as "unsuitable""...or something like that.

Hope this helps.

EDIT:

To clarify points in the ensuing discussion, I am adding here the following code with detailed comments.

import numpy
from matplotlib import pyplot as plt

#Some simple setup
signalDuration = 4; #Seconds (METERS in your case)
samplingFrequency = 100; #Samples Per Second (NUMBER OF SAMPLES ALONG YOUR SONAR PROFILE)
timeVector = numpy.linspace(0,signalDuration,samplingFrequency * signalDuration) #Discrete time (DISCRETE SPACE IN YOUR CASE)
phaseVector = 2.0 * numpy.pi * timeVector #Discrete phase. I am doing this to get rid of mul by 2*pi all the time, don't worry too much.

#Some Signal (the actual shape of it doesn't matter too much here)
Q = 7.0
A = 1.0
B = 1.0
C = 4.0

#With S here, we simulate your readings. The actual signal doesn't matter
S = Q + (A * numpy.sin(8.7*phaseVector) + B * numpy.sin(8*phaseVector)) + C * numpy.sin(2 * phaseVector)

#Let's look at our signal
plt.plot(S);plt.grid(True);plt.title("S in the time domain")

#Fourier Transform (this is the original spectrum)
FS = numpy.fft.fft(S)

#AMPLITUDE SPECTRUM. NOTICE THE SCALING FACTOR!
FS_AMP = 1/float(len(S)) * numpy.abs(FS)

#In a semilogy axis so that we can see all harmonics
plt.figure()
plt.semilogy(FS_AMP);plt.grid(True);plt.title("S in the frequency domain (AMPLITUDE SPECTRUM)")

#In a linear axis so that we can look at the amplitudes of each individual harmonic
plt.figure()
plt.plot(FS_AMP);plt.grid(True);plt.title("S in the frequency domain")

plt.show()

#The first spectral line is your C (because it is at the lowest frequency).
#The second spectral line is your B
#The third spectral line is your A
#
#Look at the spectral line that corresponds to C. We specified 4.0. Fourier Transform tells us 2.0.
#WE SPECIFIED 4.0 AND WE DO GET ABSOLUTE AMP OF 4.0 (remember the symmetry of FFT), similarly for the rest.
#Now, if you ignore Q (which is just a constant that I used to offset the measurements on purpose)
#Harmonic C is approximately 66% of the total amplitude.

Hope this helps.


Think of being at the center of sphere with a very large radius that has a very large number of uniformly distributed noise sources on its interior surface. What you receive at the center ls spherically isotropic noise.  One can envision a similar cylinder for cylindrical isotopic noise.

These are ambient noise models.  The interesting thing is that the noise is spatially correlated between two points and has a sin(x)/x form for spherically isotropic noise.

This is a model typically used in SONAR where the electronic thermal noise of each sensor is lower than the ambient noise.

This model is derived in a number of books or at least referenced in just about every SONAR or Underwater Acoustics book.  I like the presentation in  Underwater Acoustic System Analysis by William S. Burdic (first edition) because he also derives the half sphere and infinite plain models. The only criticism of Burdic's analysis is that it assumes straight line propagation. Van Trees, Detection Estimation and Modulation Theory, Vol 4, Array processing has a section on ambient noise fields with many references.

The really interesting thing of spherically isotropic noise is that for a given temporal frequency, the zeros of the spatial correlation function are separations of $\lambda/2$ which is also the spatial Nyquist criteria.


Your problem can be restated as the following:

"Find an easily computable bijection from the integers $0 .. {n \choose n / 2} - 1$ to the set of binary strings of length $n$ with exactly $n / 2$ zeros and $n / 2$ ones."

Formulated this way, this looks like a combinatorial numbering problem.

Decoding is quite simple and involves accumulating binomial coefficients for each bit set in the codeword. The encoding looks more complicated. Having enough ROM/RAM to store a big binomial coefficients table helps a lot in any case. You might want to repost on cs / math stackexchange, mentioning only the combinatorics aspect (not the digital communications application) if you are in search of more memory-efficient algorithms.

See python code for encoding and decoding here.

High-level explanation.

Let us consider the case $n = 6$ (3 ones and 3 zeros). The goal is to find a mapping between the strings 000111, 001011, 001101 up to 111000... and the numbers 0, 1, 2, 3 up to 19. This is equivalent to coding $\log_2 20 = 4.32$ bits into a 6-bit DC-free code.

There are ${6 \choose 3} = 20$ strings made of 3 zeros and 3 ones. ${5 \choose 3} = 10$ of them start with a 0 and ${5 \choose 2} = 10$ of them start with 1. So we can assign the numbers 0 to 9 to the first group, and 10 to 19 to the second group.

Within the first group, there are ${4 \choose 3} = 4$ which start with 00 and ${4 \choose 2} = 6$ which start with 01. So we can assign the numbers 0 to 3 to the 00 group and 4 to 9 to the second group.

Within the 00* group, there are ${3 \choose 3} = 1$ which start with 000 and ${3 \choose 2} = 3$ which start with 001. So we can assign the number 0 to the first one (000111), and the numbers 1 to 3 to the second group (001*). And so on...

Observe that in all the ${n \choose k}$ we have written, $n$ is always the number of remaining positions, and $k$ the number of unallocated ones. So the partitions of the $0 .. 19$ set we have built are created by adding ${\text{position at which we allocate a one} \choose \text{number of remaining ones}}$ terms.

Let us formalize this procedure. If the bitstring to decode is $b_5b_4b_3b_2b_1b_0$, the number we have associated to it is:

$b_5 {5 \choose b_4 + b_3 + b_2 + b_1 + b_0} + b_4 {4 \choose b_3 + b_2 + b_1 + b_0} + b_3 {3 \choose b_2 + b_1 + b_0} + b_2 {2 \choose b_1 + b_0} + b_1 {1 \choose b_0} + b_0$.

That's all for the decoding (DC-free bitstring to number)!

The encoding (number to DC-free bitstring) works by trying to "fit" the terms of this sum to the integer to encode. For example, given the number 15 to encode:

This number is above ${5 \choose 2}$ so the first digit is 1, and the first term of the sum is $1 {5 \choose 2} = 10$. 2 ones remains to be allocated. 5 remains in the sum.
The remainder is less than ${4 \choose 2} = 6$ so the next digit can't be a 1, so the second term of the sum is $0 {4 \choose 2}$. 2 ones remain to be allocated. 5 remains in the sum.
The remainder is above ${3 \choose 2} = 3$, so the next digit is 1, and the third term of the sum is $1 {3 \choose 2}$. 1 one remains to be allocated. 2 remains in the sum.
The remainder is equal to ${2 \choose 1} = 2$, so the next digit is 1, and the fourth term of the sum is $1 {2 \choose 1}$. All ones are allocated. All terms of the sum are accounted for.
The remaining terms are $0 {1 \choose 0} + 0$.

Thus, the bitstring associated to 15 is 101100.


Relationship between the autocorrelations of X(t) and X(nt)

The answer to the OP's question is more straightforward than rb-j's comments make it out to be.

$\{X(t)\colon -\infty < t < \infty\}$ is a continuous-time WSS random process with autocorrelation function
$$R_X(\tau) = E[X[t)X(t+\tau)], -\infty < \tau < \infty.$$
On the other hand,
$\{Y(n)\colon n \in \mathbb Z\}$ where $Y[n] = X(nT)$ is a discrete-time process whose mean function
\begin{align}
\mu_Y[n] &= E\left[Y[n]\right]\\
&= E[X(nT)]\\
&= \mu_X(nT)\tag{1}\\
&= \mu_X
\end{align}
is a constant since $\{X(t)\}$ is assumed to be WSS and so its mean function $\mu_X(t)$ has fixed value $\mu_X$.  Similarly, the autocorrelation function of $\{Y[n]\}$ is
\begin{align}
R_{Y,Y}[m,n] &= E\left[Y[m]Y[n]\right]\\
&= E[X(mT)X(nT)]\\
&= R_{X,X}(mT, nT),\tag{2}\\
&= R_X((n-m)T)\end{align}
that is, the value of $R_{Y,Y}[m,n]$ depends only on the difference $n-m$ of the arguments $m$ and $n$, and not upon their individual values. We conclude that

  $\{Y(n)\colon n \in \mathbb Z\}$ is a discrete-time WSS random process with mean $\mu_X$ and autocorrelation function
  $$R_Y[m] = E\left[Y[n]Y[n+m]\right] = R_X(mT).$$

Note that just as the random variables $Y[n]$ are regularly spaced samples ($T$ seconds apart) of the $\{X(t)\}$ process, the autocorrelation function $R_Y[m]$ is just regularly spaced samples ($T$ seconds apart) of the autocorrelation function $R_X(\tau)$ of the $\{X(t)\}$ process. Indeed, Eqs. $(1)$ and $(2)$ show that if $\{X(t)\}$ were not WSS, then the mean function $\mu_Y[n]$ would be equal to $E[X(nT)] = \mu_X(nT)$, that is,  the mean function of $\{Y[n]\}$ is just regularly spaced samples of the mean function $\mu_X(t)$ of the $\{X(t)\}$ process while $$R_{Y,Y}[m,n] = E\left[Y[m]Y[n]\right] = E[X(mT)X(nT)] = R_{X,X}\left(mT, nT\right),$$
that is, the autocorrelation function $R_{Y,Y}[m,n]$ is just regularly spaced sample values of $R_{X,X}(t_1, t_2) = E[X(t_1)X(t_2)]$ where the samples are on a square grid.

Since the OP's original question and my answer have resulted in a flood of comments and queries from rb-j that are far too many to be answered in comments, I will add something here to further explain my view of the issues, and then ignore any follow-up queries or comments. If you don't like my answer, down-vote it or write your own answer explaining your view of the matter. Both actions  will be ignored by me.

The OP's question concerns sampling a wide-sense-stationary (WSS) continuous-time process $\{X(t)\colon t \in \mathbb R\}$ to create a discrete-time process $\{Y[n]\colon n \in \mathbb Z\}$ where $Y[n]$ is defined to be $X(nT)$ where $T$ is a known positive constant.  The OP then asks for the relationship between the autocorrelation functions $R_X(t)$ (continuous-time)  and $R_Y[n]$ (discrete-time) of the two processes.

My answer, which I have explained in some detail, to this question,  is that $R_Y[n]$ is just the sampled version of $R_X(t)$:  $$R_Y[n] = R_X(nT)~\text{for all} ~n \in \mathbb Z.$$
This relationship is not a two-way street; given the values of $R_Y[n]$, it is not possible to reconstruct $R_X(t)$. Nor have I claimed anywhere that $R_Y[n]$ and $R_X(t)$ are the same function.  "But, but, but,$\ldots$ what about the sampling theorem? What about ergodicity? What about $\ldots$," you splutter. Please re-read the sampling theorem, especially the hypotheses. If $R_X(t)$ is a strictly band limited function (and nowhere has it been asserted that $R_X(t)$ enjoys this property), then it can be reconstructed from its sample  values provided the sampling rate is adequately high (exactly how high is not relevant). Given only sample values, one can construct infinitely many possible functions that have the specified sample values. If the  sample values are of a function is indeed strictly band limited, then one of these reconstructions will coincide with the function that was sampled. But simply to use the Shannon-Whittaker interpolation formula on $R_Y[n]$ and insist that the result must be $R_X(t)$ is something that boggles my mind, though ymmv. As for ergodicity, the less said the better.


Difference between MATLAB zp2sos and python zpk2sos functions

I am not sure why the first 3 coefficients in python are always different than that of MATLAB.

It's just a scaling. If you call Matlab's zp2sos without the second output argument, you will get the same three values as Python. If you need Python to match Matlab, use the first value of the sos as your "k" and divide the first three values in the sos by k (in your example this yields 1 2 1).

Output ordering for 'p' is different than MATLAB but all values account for. Not sure if it matters?

Depending on your application it actually does matter. For rather complicated reasons, Python's is better and will in most cases generate less numerical noise  and often significantly so.


"Normalizing" my D4 wavelet transformation at each step reduces final image quality

The right answer is you have to perform each of the Update/Predict steps on the input signal completely before doing the next Update/Predict.  What I was doing was walking through the signal, and performing each Update/Predict as I go.

On page 158 of "Ripples", there is a reference implementation.

// s is the signal
#define IEVEN (2*j)
#define IODD (2*j + i)
for( int i = 1 ; i <= n/2 ; i *= 2 )
{
  for( int j = 0 ; j <= n/2 - i ; j += i ) // Must do this Predict step COMPLETLEY
    s[ IEVEN ] += √3 * s[ IODD ] ;

  for( int j = 0 ; j <= n/2 - i ; j += i ) // Then this one..
  {
    int prevEvenIndex = IPREVEVEN ;
    s[ IODD ] -= d4c1*s[ IEVEN ] + d4c2*SAFE_PREV(s,prevEvenIndex) ;
  }

  for( int j = 0 ; j <= n/2 - i ; j += i )
  {
    int nextOddIndex = INEXTODD ;
    s[ IEVEN ] -= SAFE_NEXT(s,nextOddIndex) ;
  }

  for( int j = 0 ; j <= n/2 - i ; j += i )
  {
    s[ IEVEN ] *= d4normEvens ;
    s[ IODD ] *= d4normOdds ;
  }
}

The 98% 0's D4 transform:


Reconstruction of audio signal from Spectrogram

If the spectrogram was computed as the magnitude of short time fourrier transforms from overlapping windows, then the spectrogram contains implicitly some phase information.

The following iterations do the job :

$$x_{n+1} = \text{istft}(S\cdot\exp(i\cdot\text{angle}(\text{stft}(x_n))))$$

$S$ is the spectrogram, $\text{stft}$ is the forward-short time Fourier transform, $\text{isft}$ is the inverse-short time Fourier transform.


What's the unit of spectral intensity in calculation of RGB tristimulus?

You are probably working not with CIE RGB, but CIE XYZ, so these are properly denoted the XYZ tristimulus values, not RGB. Their absolute values do not matter if you want to determine the chromatacities ("color", as we commonly understand it), because you take the ratio of the components to their sum.

Once you have the XYZ tristimulus values or xy chromaticities, you can convert to something like sRGB to yield usable RGB values. You might need to scale the result to 0-255.


Your expected output can be achieved when the input is as follows:

real part = cos(2*pi*[0:63]*3/64)

imag part = 0

So I think you have a couple of issues:

You are missing a factor 1/64 in your description of the input above.
You should use cos() instead of sin() for the real part of the input.

Try fixing issue 1 in your post, and try fixing issue 2 in your code.


inverse discrete fourier transform with plain python

Removing the 2nd term -> array = np.asarray(array, dtype=float) should work


YV12 to RGB - What is wrong with my algorithm?

The output image has lines shifted one against another, which is typically due to incorrect [byte] distance between the lines (stride). That is, the color conversion math is perhaps correct, however you are either taking source pixel from wrong position and/or putting output pixel into wrong place.

For performance reasons YV12 images might be using extended strides in hardware to facilitate performance. For example, 320x240 YV12 image might have a stride of 512 and this means that the data buffer is a 512x240 data space (512x240 bytes for Y, followed by two 256x120 arrays for V and U) where the payload image on the left, and there is invisible padding on the right.


Amplitude of the fft is less than what is expected

The amplitude peak is hidden between FFT result bins for frequencies that are not exactly integer periodic in the FFT's length.  But you can interpolate to find the results for those frequencies.  Try windowed Sinc interpolation of the real and imaginary components, then compute the magnitude of the complex result.  If you don't know the exact frequency, you can use successive approximation to find an interpolated magnitude peak between FFT result bins.


Entropy can be used as a criterion to decide whether a wavelet packet branch should be further decomposed, or to prune a wavelet packet tree. This carries similarities with segmentation techniques like split-and-merge.

Let us take one sign convention for entropy. One hopes that entropy diminishing, locally, across some levels, meas that this part of the image is simplified by
 one wavelet decomposition, and when the entropy starts to increase (or vice-versa), further decomposition levels become useless to simplification. Following the entropy gradient is a way to find some "optimal" wavelet packet splitting.

But first,  optimal for what purpose is unclear from your question so far. Denoising, compression? Entropy can provide a nice guess to predict compression performance with an actual coder, but should be applied to binned/quantized wavelet coefficients, and may be limited unless you use higher-order/cross-entropy to take into account advanced vector coders (as compared to scalar ones).

Second, traditional wavelets packets for image are implemented separably, and notoriously, separable schemes and packets lack in efficient representation for  piecewise regular images, and induce increased aliasing, that may yield mildly efficient decompositions.

Codes implementing such decompositions are presented in Matlab, see for instance: wpdec2: Wavelet packet decomposition 2-D, and the reference paper Coifman, R.R.; M.V. Wickerhauser (1992), “Entropy-based algorithms for best basis selection,” IEEE Trans. on Inf. Theory, vol. 38, 2, pp. 713–718.


What does it take to design a 4th order band-pass filter with -24dB/Octave?

After experimenting with winfilter, generating the following filters mimics the EFX isolators, amazingly :)

2 * Butterworth LP 2nd order @ 300Hz
2 * Butterworth BP 2nd order @ 300Hz/6000Hz
2 * Butterworth HP 2nd order @ 6000Hz

So each signal passes through two filters of the same type, for some reasons I don't know the cause of, a 4th order BP in winfilter simply doesn't work and produces garbage instead.

As a reference, here's what really helped me out in matching the filters properties, white noise recorded for each of the bands then guesses until it matches:

Surprisingly, I'd say the above filters match the EFX about 99%, the only differences are at the edges where the winfilter ones really trim down the frequencies much more.


I can envision advantages (or convenience) in doing this for lower IF signals certainly, as long as the sampling rate is significantly higher than the signal bandwidth to contain the phase dispersion you would get with this approach.  An implementation can be lower cost and simpler in the analog with a single mixer/ADC front-end for example. This becomes less feasible at higher carriers where nonlinear noise contributions can start to dominate (relative SFDR in ADC's for example).

To further comment on cautions with an approach of sampling at "twice the carrier frequency". I assume in this statement that the approach is to approximate the baseband IQ samples by demultiplexing the samples: I1 Q1 I2 Q2 I3 Q3, etc.  It should be pointed out that this is just an approximation, with large errors as the signal bandwidth becomes large relative to the entire digital bandwidth. This is because in this approach the 90° shift is being emulated by a delay, but a delay would only be 90° at ONLY the center frequency of the signal (and varies linearly with frequency everywhere else). To do IQ sampling without delay distortion requires a Hilbert transform which is 90 degrees at all frequencies within the bandwidth of the signal. Note this is of minimal consequence when the bandwidth is significantly smaller than the sampling rate, which would apply to most direct carrier sampling schemes, but is mentioned since the real application to do 2x carrier sampling would not be such cases but lower IF carriers where analog cost and performance is less affected.

So that said, I would likely never sample at twice the carrier frequency, but if the intention is to sample close to that (for some reason), I would choose either 2/.75 = 2.66x or 2/1.25 = 1.6x. Let me explain why below:

So how to do something close to the intention of 2x carrier sampling (which is to sample at a different frequency closest to fs with a single ADC and avoid having to do any IQ in the analog). Obviously we can envision approaches of sampling with a much higher rate to a digital IF and processing IQ in the digital, but to quantify a specific approach that has no additional inaccuracies beyond the analog hardware is to to the following:

Sample the analog signal following all rules of Nyquist (including anti-alias bandpass filtering), meaning ensure the sampling rate is twice the bandwidth of the anti-alias filter, with the additional conditions:

Sample the analog signal such that the resulting digital signal resides at $F_s/4$ where $F_s$ is the sampling rate or any integer fraction of the actual sampling rate.

This means the signal can be centered in the analog prior to sampling at $nF_s/4$ where n = 1,3,5,...    Which n is used is referred to as which "Nyquist Zone". n=1 is the first and primary Nyquist zone, for example.

Once digitized, the signal can be complex down-converted to baseband if complex IQ samples are still desired.  What is great about positioning the signal at $F_s/4$ is that the complex downconverter reduces to multiplying the signal by the following to generate I and Q:

Since the complex rotator to do $Fs/4$ downconversion to baseband is multiplying by $e^{-jn\pi/4}$ = +1, -j, -1,+j,... then the I samples for this multiplier are +1,0,-1,0,... and the Q samples are 0, -1, 0, +1 ... So baseband I and Q is created by doing the following to the sampled sequence:

$$I = +1x_1 + 0x_2 -1x_3 +0x_4 + 1x_5 ...$$

$$Q = 0x_1 - 1x_2 +0x_3 +1x_4 + 0x_5 ...$$

Very convenient: demultiplex and multiply each by +/-1!

Note, if even Nyquist zones are chosen for the analog downconversion, spectral inversion will result. This is easily resolved by rotating the other way: $e^{+jn\pi/4}$.

$$I = 1x_1 + 0x_2 -1x_3 +0x_4 + 1x_5 ...$$

$$Q = 0x_1 + 1x_2 -0x_3 -1x_4 + 0x_5 ...$$

Note how this is done just by changing the sign of the Q.

Now if the complex rotator is used for carrier tracking (which is entirely likely) then we would not use the reduction down to +/-1's in the multiplier, but still benefit by the ideal placement of fs/4 which optimizes the analog anti-alias filter design by centering the signal within the bandwidth from the edges where aliasing would otherwise occur.


Can multi-rate operations of decimation and interpolation can be used to implement a rate change by a factor an irrational number?

You're exactly right, with rational operations, you can't achieve an irrational ratio.
The only way out there is using some interpolative number, i.e. an arbitrary resampler.
However, in signal processing, you basically always deal with system that only run for a finite time and within that have a boundable tolerance against jitter or frequency error - so one common solution to that dilemma is, e.g. in receivers, to do the resampling with a slightly wrong, but rational ratio, and then let the symbol synchronizer, which often does some form of arbitrary resampling internally, anyways, take care of the rest.


discrete-time (slew) rate limiting algorithm

Like the old German saying goes:

If you give a man a filter, everything looks like a signal that needs to be low-pass filtered

or so.

Slew rate and signal bandwidth are closely related. Simply low-pass filter your signal; usually, you'd do that with a linear filter, not an non-linear filter, and especially not a non-linear recursive one, because those are pretty hard to get stable.

By the way, if you look at FIR and IIR filters, you'll notice that we model these in z-Domain as polynomials, but by doing so, we can apply them to the signal in a linear way.


Question based on scaling the FFT for obtaining PSD

The exact location (index) of the Nyquist frequency in your frequency vector depends on the length of the input signal. If the length is odd the Nyquist frequency is not included; if the length is even it is included. Consequently, and assuming the signal is real, either the first $\frac{N}{2}$ or $ \frac{N}{2} + 1$ samples are needed to fully describe your signal in the frequency domain. Note that, like the first sample describing the DC value, the sample describing the Nyquist value only occurs once; it should therefore not be multiplied by two when calculating the PSD of a real signal.

The proper scaling value to be used depends on the exact definition of the DFT. I have not checked MATLAB’s definition but if you want to make sure that the used scaling factor is correct use Parseval’s theorem. Check if the signal power in the time and frequency domain are identical; if not adjust the scaling factor accordingly.


Approximation using a Fourier transform with low pass filter

Depends on what you mean by 'best'.  In terms of least squares, then yes, just knocking off the terms above 1kHz will give the closest approximation.  I guess one way to prove it would be:
$$
     x(t) = \int_{-\infty}^{\infty} X(f)\exp\left( j2\pi ft \right) dt\\
\hat x(t) = \int_{-1k}^{1k} X(f)        \exp\left( j2\pi ft \right) dt
$$
leaving a residual error of
$$
e(t) = \int_{|f| > 1k} X(f)\exp\left( j2\pi ft \right) dt .
$$

Compare this to a 'fiddled' version, using $\tilde X(f)$ for the coefficients:
$$
\hat x(t) = \int_{-1k}^{1k} \tilde X(f) \exp\left( j2\pi ft \right) dt,
$$
which has a residual error of
$$
e(t) = \int_{|f| \le 1k} \left(X(f) - \tilde X(f)\right) \exp\left( j2\pi ft \right) dt +
\int_{|f| > 1k} X(f)\exp\left( j2\pi ft \right) dt ,
$$
which is greater than the error in the first instance unless $\tilde X(f) = X(f) \; \forall f$.

Something like that ought to do it, anyway.

EDIT: It is possible to prove that a projection onto any orthogonal set of functions (doesn't need to be the Fourier set) gives the least squares fit:

$$
\int f_k(t) f_m(t) dt = \cases{
  0, & $k \ne m$ \\
  1/\lambda_k, & $k = m$}
$$
for any orthogonal basis set $\{f_k\}$ (I'm not bothering with weighting functions), where
$$
\lambda_k = \frac{1}{\int f^2(t) dt}
$$
The projection of $g(t)$ onto the set is given by
$$
c_k = \lambda_k \int g(t)f_k(t) dt.
$$
The least squares fit of a set of orthogonal basis function is the minimisation of
$$
\int \left[g(t) - \sum_k d_k f_k(t)\right]^2 dt.
$$
Minimise by setting the derivative (wrt $d_m$) to zero:
$$
2\int\left(g(t) - \sum_kd_kf_k(t)\right)\left(-f_m(t)\right) dt = 0 \\
-\int f_m(t)g(t)dt + \int f_m(t)\sum_kd_kf_k(t)dt = 0 \\
-\int f_m(t)g(t)dt + \sum_kd_k \int f_m(t)f_k(t)dt = 0 \\
-\int f_m(t)g(t)dt + \sum_kd_k \frac{1}{\lambda_m} = 0
$$
which rearranges to give the same value as for $c_k$, thus the least squares fit is exactly the straightforward projection onto the orthogonal set.


Third octave bandpass filter with python

I recently developed a function to easily perform octave and fractional octave filtering, it is available on github: PyOctaveBand
It uses the SOS coefficients and performs downsampling to filter correctly at low frequencies.


Convolution in frequency space with fft and ifft -> num output samples?

What you do is called circular convolution but you want "linear" convolution: Say you have two signals of length $N$ and $M$. Append zeros to both time-domain signals, so that their length is $L=N+M-1$, then calculate the $L$-DFT of both signals, multiply them in frequency-domain and convert the result back to time-domain.

In your case, you should append 44099 zeros to each signal before transforming it to the frequency domain. The result will then be 88199 samples long.

If the available FFT is shorter than 88199 you have to divide the signal into shorter blocks and apply the overlap add method.


sequence to perform filtering? which filter first?

Theoretically (with infinite precision arithmetic) order of LTI filters does not matter. If any of the filters is not LTI then order will be important.

Plus, for finite precision arithmetic one should also take care of numerical error build up and potential instability issues especially with IIR filters where pole locations may become critical


Which Anisotropic Diffusion Method is Best for Denoising Now?

I think "Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE’s" by David Tschumperlé is really good.

I started implementing it using MATLAB once, yet I didn't finish.
If you do implement it, it will be great if you updated us.

Update

Have a look on my implementation at Fast Anisotropic Curvature Preserving Smoothing.


If your noise has independent and identically distributed samples from a zero-mean distribution (for example Gaussian), it is white. Other definitions of white noise also require the distribution to be symmetrical, but that is not required for the spectrum to be flat. Clipping the samples of such white noise will only change the common probability distribution of the samples:

This will affect the mean of the distribution, but if both sides are clipped in a way that preserves the zero mean of the distribution, then there is no "coloring" of the spectrum by a 0 Hz peak and the noise remains white.

Example in Octave (MATLAB clone):

Create random variables from a Gaussian distribution:

N = 65536;
x = normrnd(zeros(N, 1), ones(N, 1));

Clip from above:

y = min(x, 1);

Clip symmetrically from below:

z = max(y, -1);

Plot values of Gaussian random variables (blue), clipped above (green), clipped symmetrically (red). Horizontal axis = index, vertical axis = value:

plot(1:N, x, 1:N, y, 1:N, z)

Plot spectrum of Gaussian random variables (blue), clipped above (green), clipped symmetrically (red). Horizontal axis = frequency, vertical axis = magnitude:

loglog(1:N/2, abs(fft(x)(1:N/2)), 1:N/2, abs(fft(y)(1:N/2)), 1:N/2, abs(fft(z)(1:N/2)));

The sequence that was clipped from above has a spectral peak at 0 Hz. Otherwise, the spectral envelopes are flat. This is not very easy to see because of the increasing density of plot points in the high frequencies in the logarithmic scale plot.


How do I remove all white pixels in a modified image which are not present in the original?

Welcome to Image processing!

(See edit at the bottom)

If I understood you well, thresholding the image with a low threshold will provide a mask for the hand that you can use to remove the background pixels.

That's how it will look in Mathematica (ImagePad here is just used to remove the extra white border on your images):

img = ColorConvert[Import["http://i.stack.imgur.com/krbMs.jpg"], "Grayscale"];
border = BorderDimensions[img];
img = ImagePad[img, -border]

hand = Binarize[img, 0.2]

thresholded = Binarize[ImagePad[Import["http://i.stack.imgur.com/ylo3B.jpg"], -border], 0.1]

You may want to check this other question for an alternative segmentation technique based on ridge filtering.

EDIT: Based on your comment, you would like to only work on the thresholded image and clean it up to leave only the components corresponding to veins. I would look at morphological operations (e.g. Opening, components selection based on shape properties). By the way I feel this is harder than working on a different vein segmentation algorithm (see one proposed above).


Why are the mark and space tones of these FSK signals reversed?

The apparent exchange of preamble symbols between mark and space in separate transmissions indicates that the packet is differentially encoded in some way.  NRZI encoding is used in both Packet Radio and AIS, so you might wish to investigate if that appears to be the case for your unknown packet format as well.


signal analysis and synthesis with different window length and FFT length

MATLAB documentation says for element-wise operators + and .*:

  Numeric or string inputs A and B must either be the same size or have
  sizes that are compatible (for example, A is an M-by-N matrix and B is
  a scalar or 1-by-N row vector).

Further explanation:

  Two inputs have compatible sizes if, for every dimension, the dimension sizes of the inputs are either the same or one of them is 1.

It looks to me like your code violates this if winLen < FFTLen.


Frequency response equivalent to characterise transient response (to sinusoidal inputs)?

As already pointed out in Hilmar's answer, a linear time-invariant (LTI) system is completely described by its impulse response, or - equivalently - by its frequency response, which is the Fourier transform of its impulse response. I think that there is a misunderstanding concerning the Fourier transform, namely that it can only be used to describe steady-state responses to infinitely long sinusoidal input signals, but this is definitely not the case. I think it's illustrative to look at an example:

Let's assume we have a causal real-valued LTI system with frequency response

$$H(\omega)=A(\omega)+jB(\omega)$$

where $A(\omega)$ and $B(\omega)$ are the real and imaginary parts of $H(\omega)$, respectively. Now let's define an input signal to this system:

$$x(t)=u(t)\sin(\omega_0t)$$

Its Fourier transform is

$$X(\omega)=\frac{\pi}{2j}[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)]
+\frac{\omega_0}{\omega_0^2-\omega^2}$$

Using $A(\omega)$ and $B(\omega)$, the imaginary part of the Fourier transform of the output signal $y(t)$ can be written as

$$\Im\{Y(\omega)\}=-\frac{\pi}{2}[A(\omega_0)\delta(\omega-\omega_0)-A(-\omega_0)\delta(\omega+\omega_0)]-\frac{\omega_0B(\omega)}{\omega^2-\omega_0^2}\tag{1}$$

Since $y(t)=0$ for $t<0$ (because this is true for $x(t)$ and because the system is causal), it is fully described by the imaginary part of its Fourier transform:

$$y(t)=-\frac{2}{\pi}\int_{0}^{\infty}\Im\{Y(\omega)\}\sin\omega t\; d\omega,\quad t>0\tag{2}$$

From (1) and (2) we get

$$y(t) = A(\omega_0)\sin\omega_0t+\frac{2\omega_0}{\pi}\int_{0}^{\infty}\frac{B(\omega)}{\omega^2-\omega_0^2}\sin\omega t\;d\omega,\quad t>0\tag{3}$$

Equation (3) fully describes the response of the LTI system to the causal input signal $x(t)=u(t)\sin\omega_0 t$, including transients. This example is supposed to show that the frequency response (or the impulse response) is sufficient for a complete description of the system's response to arbitrary input signals.

EDIT: Derivation of Equation (2):

Let $Y(\omega)=Y_R(\omega)+jY_I(\omega)$ be the Fourier transform of $y(t)$.
We need the following basic properties of the Fourier transform (see here):

$$y(t)\textrm{ real-valued}\Longrightarrow Y(\omega)=Y^*(-\omega)$$
from which follows
$$Y_R(\omega)=Y_R(-\omega)\textrm{ and }Y_I(\omega)=-Y_I(-\omega)$$

If $y_e(t)=0.5[y(t)+y(-t)]$ is the even part of $y(t)$ and $y_o(t)=0.5[y(t)-y(-t)]$ is the odd part of $y(t)$ then

$$y_e(t)\Longleftrightarrow Y_R(\omega)\\
y_o(t)\Longleftrightarrow jY_I(\omega)$$
These relation can be found in the link above.

If $y(t)$ is real-valued the inverse transform can be written as

$$y(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}Y(\omega)e^{j\omega t}d\omega=
\frac{1}{2\pi}\int_{-\infty}^{\infty}\Re\{Y(\omega)e^{j\omega t}\}d\omega=\\
=\frac{1}{2\pi}\int_{-\infty}^{\infty}[Y_R(\omega)\cos(\omega t)-Y_I(\omega)\sin(\omega t)]d\omega=\\
=\frac{1}{\pi}\int_{0}^{\infty}[Y_R(\omega)\cos(\omega t)-Y_I(\omega)\sin(\omega t)]d\omega$$

where the last equality follows from the fact that the integrand is an even function of $\omega$.
Consequently we have

$$y_e(t)=\frac{1}{\pi}\int_{0}^{\infty}Y_R(\omega)\cos(\omega t)d\omega\\
y_o(t)=-\frac{1}{\pi}\int_{0}^{\infty}Y_I(\omega)\sin(\omega t)d\omega$$

With these preliminaries we can finally state the result using the fact that for causal $y(t)$ the following holds for $t>0$ (because $y(-t)=0$ for $t>0$):
$$y(t)=2y_e(t)=2y_o(t),\quad t>0$$

This means

$$y(t)=\frac{2}{\pi}\int_{0}^{\infty}Y_R(\omega)\cos(\omega t)d\omega=\\
=-\frac{2}{\pi}\int_{0}^{\infty}Y_I(\omega)\sin(\omega t)d\omega,\quad t>0$$

The last equality is the same as Equation (2) above.


Is this statement contradicting with the band-limiting theorem?

there is some choice of words i would do differently.

yes truncating the infinite bandlimited function makes it time-limited and then band-unlimited.  but truncating it to an integer number of periods still allows one to fully describe it with discrete samples.

a periodic function can be described with a Fourier series.  the Fourier coefficients will virtually entirely describe the periodic function.

a bandlimited and periodic function can be completely described with a finite set of Fourier coefficients.  also a bandlimited and periodic function can be completely described by sampling a single period with a sufficient number of samples that are spaced closer in time than $\frac{1}{2B}$ where $B$ is the bandlimit.


As Marcus have pointed, the codec in that article refers to a narrower meaning of the term codec specifically applied to digital (computer) hardware systems.

The term codec is acronym for coder and decoder, hence a codec consists of two separate parts (at different ends) working together to perform an information processing task in between.

A coder is an algorithm (implemented in digital hardware, software or analog hardware) that is used to convert a source of information into a form that has some advantages for the purpose of representation, storage or transmission etc. The source after being processed by the coder is said to be as encoded. The most typical application for coders is compressing digial media such as music, speech, images, video so that they can be represented with less number of bits then their raw formats which is advantageous for storage space and transmission bandwidth for a fixed transmission time such as in real time communications or equivalently transmission time for a given fixed bandwidth for file downloading from internet. The algorithm of the coder is specified mathematically, and it can be implemented in software such as those dll files in a Windows operating system or inside specific program folders. It can also be implemented in hardware which is the case for such as the digital cameras where the hardware device is required to convert the captured raw video input stream into a chosen video format for storage purposes as fast and efficient as possible. A software codec running on a general purpose CPU could be weak for such embedded applications. Note that historically analog hardware base encoders were also used for transforming raw media analog inputs into feasible analog formats, the prime example being a TV encoding IC which converts RGB raw analog video stream into an NTSC/ PAL/ SECAM video format which enables a low bandwidth and black and white compatible transmission over RF broadcasting.

The decoder does the inverse; given an encoded data stream, it produces the original raw data stream as accurate as possible, exactly if the codec is a lossess type and approximately if lossy type. The decoder can also be implemented in hardware or software. The hardware implementation is always faster wheres the software is always more flexible.

As a general principle, the coder and decoder pairs are designed such that the encoders are more powerful and expensive, so that decoders can be less powerful and cheaper (form a commercial perspective) which is the case in multimedia broadcasting systems.

Coming to your figure, it's an accepted misusage of the term codec. Where the input audio from the digital microphone (which is another accepted misusage) is converted into a suitable form for representation and transmission in the following digital systems and vice versa. Hence conventionally an ADC chip which converts analog audio into digital PCM stream and a DAC chip shich converts digital PCM stream into analog audio output may be termed as a codec in the digital audio industry.


Scrambler generator polynomial with initial state

This generating (check if it's a primitive polynomial - those are usually used) polynomial representing a simple Linear Feedback Shift Register (LFSR). The circuit is given below, the gate is a XOR gate (an addition over GF(2) ). You can follow the output of this shift register very easily


A zero phase moving average filter is an odd length FIR filter with coefficients

$$h[n]=\begin{cases}\frac{1}{N},&\quad -(N-1)/2\le n\le (N-1)/2\\0,&\quad\text{otherwise}\end{cases}$$

where $N$ is the (odd) filter length. Since $h[n]$ has non-zero values for $n<0$, it is not causal, and consequently, it can only be implemented by adding a delay, i.e. by making it causal.

Note that you can't simply use Matlab's filtfilt function with that filter because even though you would get zero phase (with a delay), the magnitude of the filter's transfer function gets squared, corresponding to a triangular impulse response (i.e., input samples further away from the current sample receive less weight).

This answer explains in more detail what filtfilt does.

See also this related question.


Discard circular/blob like structures from image

You can calculate the area of each connected area in your image. Then, you throw away the smaller areas and only keep the biggest. See this example Python code:

import skimage.io as io
import numpy as np

from skimage.measure import regionprops, label
from skimage.color import label2rgb

I = io.imread("http://i.stack.imgur.com/rtyFJ.png")
bw = I[:,:,0]

bigger = np.zeros_like(bw)
labels = label(bw)

for R in regionprops(labels):
    if R.area > 500:
        # draw the region (I'm sure there's a more efficient way of doing it)
        for c in R.coords:
            bigger[c[0], c[1]] = 1

plt.imshow(bigger)

You can play around with the area threshold, e.g. by calculating a histogram of all areas and choosing the 3 biggest etc. And, you need to define the neighborhood of a pixel to belong to a connected region. In my case, it uses 8-neighborhood (i.e. if two pixels are connected at a corner, they belong to the same region). This is up to your application.

Additionally, you might throw away regions based on other criteria. E.g. if you want to throw away "circular" structures, you might throw regions away, where the area and perimeter are close to pi for example. Or, you look at the extent or orientation or whatever.


Adaptive filtering: Optimum filter length and delay

In order to be able to choose an optimal value for the delay $\Delta$ it's important to understand how the system works. The purpose of the delay is to decorrelate the desired signal $s(n)$ and the signal component $s(n-\Delta)$ at the input of the adaptive filter. This means that $\Delta$ must be chosen such that the autocorrelation $R_{ss}(k)$ of $s(n)$ is (close to) zero for lags greater than $\Delta$:

$$R_{ss}(k)\approx 0,\qquad |k|>\Delta$$

However, we cannot choose $\Delta$ arbitrarily large because the delayed interference at the input of the filter must be correlated with the interference added to the signal, i.e., the autocorrelation $R_{rr}(k)$ of the interference must still be significant at a lag of $\Delta$, otherwise the adaptive filter cannot predict the interference. If we can assume that $r(n)$ is narrow-band compared to $s(n)$, it's always possible to find an appropriate value for $\Delta$.

With an appropriate value for $\Delta$, the adaptive filter will try to predict the interference, i.e., it will try to undo the effect of the delay in the frequency band where the interference has significant frequency components. So the output of the filter will approximate $r(n)$: $y(n)\approx r(n)$. Consequently, the error signal will approximate the desired signal: $e(n)\approx s(n)$.

After having chosen a value for $\Delta$ based on the autocorrelation of $s(n)$, the filter length must be chosen by trial and error. A long filter will give a better suppression at the cost of slower convergence.


How to calculate transfer rate and needed bandwidth

Your calculations are wrong. You have misunderstood the concept. The files may contain an amount of bits as you calculated but transferring that information is a bit different.

First of all, you send information as PACKETS. Here is an example packet type which is being used in Bluetooth.

The only place you send the actual information is called PAYLOAD. As you can see, payload is occupying a smaller area than you expected in a packet. The size of payload can smaller or larger depending which type of packet you are using, or which type of standards you are working with etc. And your bit rate is determined by which modulation type you use, how many bits your one symbol contains etc.

And not to forget, if your conditions are not perfect (which is the general situation) you WILL send some packets which is corrupt and you need to send them again. Which also reduces your actual bit-rate.

So there is different calculations for different combinations of the different techniques.

I hope the answer is clear.

As always,
Have a nice day.


How do you find the phase of the DFT of a discrete signal?

Here is a hint that will help you:

The DFT is cyclical in time and in frequency. For the sequence given by

$$x(n) = [-1,2,-3,2,-1]$$

With x(0) = -3 would be solved using the standard DFT equation that starts at n=0 using

$$x(n) = [-3, 2, -1,-1, 2]$$

From that you can solve for the DFT and then determine easily for each result what it's phase is.


script to check zero crossings to check if microphones phases are aligned

A better approach for use in the presence of noise and distortion is to use the Wiener-Hopf equations which will provide a least-mean-square solution of the effective "channel" between the microphones. The group delay can be determined using scipy.signal.group_delay.

Further details on this approach specific to microphone captures including Matlab/Octave source code (that can be easily ported to Python Scipy) is on this post:

Compensating Loudspeaker frequency response in an audio signal

In that post, the right channel was equalized to the left channel by passing it through the derived compensation for the effective channel between the microphones, but the post details how you can determine the channel itself rather than the compensation.

You would want to do this with sounding signals that fully occupy the usable spectrum of interest (chirps or white noise generators would be good choices) since the channel can only be determined accurately at the frequencies where signal energy exists. You can do this using one received channel compared to the other as I have done in the link when you don't have a copy of the source, or you can do each left and right channel separately compared to the source if you have the exact waveform that was used to create the source sound.

By using this to derive the channel, and then using the group delay on the channel frequency response, you will see the delay for each frequency in your signal (if there is significant variation, this is group delay variation which is a source of distortion (so you can also use this technique to equalize your room when you do have the source signal). To match the distance between the microphones, you would compare the average of the group delay over frequency, and you can weight this average on frequency ranges that are more important to you.


First you need to compute the signal's autocorrelation:

$$R(k)=\sum_nx(n)x(n+k),\quad k=0,1,\ldots,p$$

In your case with $p=2$ you get

$$\begin{align}R(0)&=(-1)(-1)+(-2)(-2)+0\cdot 0+1\cdot 1=6\\
R(1)&=(-1)(-2)+(-2)\cdot 0+0\cdot 1=2\\
R(2)&=(-1)\cdot 0+(-2)\cdot 1=-2\end{align}$$

Then you have to solve a system of linear equations:

$$\begin{bmatrix}R(0)&R(1)\\
R(1)&R(0)\\
\end{bmatrix}
\begin{bmatrix}a_1\\a_2\end{bmatrix}=
-\begin{bmatrix}R(1)\\R(2)\end{bmatrix}$$

which gives the values $a_1=-0.5$ and $a_2=0.5$. Of course $a_0=1$, as you've already noted.

The matrix that occurs in the system of linear equations is a Toeplitz matrix, i.e. all its diagonals have the same elements. As a result, such systems can be solved by specialized and efficient algorithms, such as Levinson's algorithm.


Higher-order moment of output of LTI system

For the case of input process $\{X(t)\}$ being white Gaussian noise with two-sided power spectral density $\frac{N_0}{2}$, the output process $\{Y(t)\}$ is a strictly stationary zero-mean Gaussian process in which all the random variables have the same variance $\frac{N_0}{2}\int_{-\infty}^\infty |H(f)|^2 \,\mathrm df$ almost as you say.  But the key point here is that the $Y$'s are zero-mean Gaussian random variables with known variance, and it is known that if $Y \sim N(0,\sigma_Y^2)$ then
$$E[Y^n] = \begin{cases}\sigma_Y^n
\cdot (n-1)
\times  (n-3)
\times (n-5)\times  \cdots \times 3 \times 1, & n ~ \text{even},\\
0, & n ~ \text{odd}.\end{cases}$$

Now, if the input process is just white noise (not necessarily Gaussian white noise), then the output process is not necessarily Gaussian, and while the mean and variance are as specified above, we cannot say anything about the distribution of the $Y(t)$'s and cannot infer anything about the higher moments of the $Y(t)$'s the way we could for Gaussian output processes.

If the input process is a zero-mean WSS Gaussian process (which is also a strictly stationary process) but not necessarily a _white noise process, than $\{Y(t)\}$ is also a zero-mean WSS Gaussian process (and so strictly stationary too) with variance
$$\sigma_Y^2 = \int_{-\infty}^\infty S_X(f)|H)f)|^2 \,\mathrm df$$ where $S_X(f)$ is the power spectral density of the input process and so once again we are in business and can use $$E[Y^n] = \begin{cases}\sigma_Y^n
\cdot (n-1)
\times  (n-3)
\times (n-5)\times  \cdots \times 3 \times 1, & n ~ \text{even},\\
0, & n ~ \text{odd}.\end{cases}$$


If X is your fft data then you need to take Y=abs(X).^2; to convert it. Your frequency axis is f=4:399;. You can plot the result then using plot(f,Y).


Spectrum of OFDM with raised cosine window

This spectral shaping technique is applied in time domain, after adding the guard interval (GI). But I find it easier to add GI and the cyclic extension for windowing in one step. Let $x(n)$ be an $N$ subcarriers OFDM symbol without guard interval. Then $W + G$ samples are copied to the beginning accounting for guard interval and windowing samples. Additionaly, $W$ samples are copied to the end, also for windowing:
$$
y(n) =
\begin{cases}
x(n+N) & \text{for} & -G-W \leq n \leq -1 \\
x(n) & \text{for} & 0 \leq n \leq N -1 \\
x(n-N) & \text{for} & N \leq n \leq N + W -1
\end{cases}
$$
In the next step, the raised cosine function is applied to the first and last $W$ samples of $y(n)$, respectively. The windowing function $w(n)$ is given by:
$$
w(n)=
\begin{cases}
\cos^2\left( \frac{n+G+1}{W-1} \frac{\pi}{2}\right) & \text{for} & -G-W \leq n \leq -G-1 \\
1 & \text{for} & -G \leq n \leq N - 1 \\
\cos^2\left( \frac{n-N}{W-1} \frac{\pi}{2}\right) & \text{for} & N \leq n \leq N+W-1 \\
0 & \text{otherwise}
\end{cases}
$$
$w(n)$ is similar but not equal to the transfer function of a raised cosine filter often used as impulse shaper in single carrier transmission systems. The two differences are: (1) the raised cosine function is applied in time domain for OFDM systems and in frequency domain for single carrier systems and (2) the "flat top" is usually much longer for OFDM systems, whereas its length is in a fixed relation with the flanks' length, given by the roll-off factor, for single carrier systems.

Finally, the OFDM symbol including GI and spectral shaping is calculated by
$$
z(n) = w(n)y(n)
$$

When transmitting several OFDM symbols $z_i(n)$, two consecutive symbols overlap at $W$ samples. The discrete transmit signal $u(n)$ is therefore given by
$$
u(n)=\sum_{i=-\infty}^\infty z_i(n-i(N+G+W))
$$

The implementation in Matlab should now be straightforward by substituting $n$ with $n' = n + G+ W+ 1$ in the above equations. Leave a comment if not.


Determining whether or not a system is BIBO stable

In your case

$$\int_{-\infty}^{\infty}|h(t)|dt=\int_{-\infty}^{\infty}h(t)dt=
\int_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}\delta(t-2n)dt=
\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}\delta(t-2n)dt$$

And since

$$\int_{-\infty}^{\infty}\delta(t-2n)dt=1$$

(as correctly pointed out by Dilip Sarwate in his comment to Phonon's answer), the above sum is

$$\sum_{n=-\infty}^{\infty}1=1+1+1+\ldots=\infty$$

I believe that this is basically the argument you gave in your question.


The phase of a sinewave that is not exactly periodic relative to the FFT aperture length is discontinuous at the edges.  If you want a sane/visualizable phase (not the phase of a discontinuity) for non-periodic frequencies, try generating your sinewave with a known phase at the center of the FFT (sin(2*pi*foo*(t-N/2))), and use an fftshift to re-reference the reported phase result of the FFT to the center of the FFT aperture as well.  Then the oddness/evenness ratio of your test function, and thus the phase result will remain stable as frequency is swept between FFT bin center frequencies.

The phase of a non-periodic-in-aperture sinewave will contain the phase of the convolution of the signal with the transform of the FFT window (something like a Sinc, which ends up spanning the full width of the FFT), which will also have a phase related to that of the signal.  Even a computed FFT of an exactly periodic sinewave will also produce some values that aren't exactly zero in all the other bins due to quantization noise in the computer arithmetic; thus those essentially zero, but not quite zero, magnitudes will have a phase that may or may not be zero.  And the phase of an actual zero magnitude vector is indeterminate, not necessarily zero.


One possible pre-processing apporach is to use imdilate(I,strel('disk',5)) to thicken the road lines a little bit.

If you are trying houghlines, probably one parameter that you need to tune is the 'FillGap' which will connect the very short segments (you can try to change this parameter yourself to observe the effect, along with the MinLength value):

lines = houghlines(BW,T,R,P,'FillGap',30,'MinLength',15);

However, the road curves are not straight lines, so I won't recommend you to use hough transform. You may try bwmorph(bw,'skel',N); after you thicken the binarized road image. N is the iteration number that you need to tune.