SPL values from FFT of microphone signal

Let's assume that signal you are analysing is sinusoid with amplitude $A$:

$x=a\sin{2\pi f_{0} t}$

Its RMS value of amplitude is then: $\dfrac{a}{\sqrt{2}}$ as you noticed in your code. Before performing DFT it is good practice to window signal as it decreases leakage, etc.

After transforming this signal into the frequency domain, you will obtain a two-sided spectrum. Then we have to take the magnitude of these complex values. Probably you would notice that all of them are very high. That's because energy is not normalised. Normalisation can be done via dividing values of magnitude by energy of your window. As people tend to take signal "as it is" (which is equal to applying rectangular window) the spectrum is then divided by a number of samples $N$. For different types of windows this factor that is equal to the sum of all window samples.

Additionally, we are only interested in one half of a spectrum, therefore the amplitude of all samples must be multiplied by $2$ to compensate the loss of energy, except DC component (which appears only once) and Nyquist bin (if it's present).

Lastly, as you want to analyse RMS of the signal, you must understand the following. All values of single-sided spectrum are of height $\dfrac{A_{i}^{2}}{2}$, which is equal to $ \left( \dfrac{A_{i}}{\sqrt{2}}\right)^{2}$, where $\dfrac{A_{i}}{\sqrt{2}}$ is the RMS magnitude amplitude for the i-th frequency component. Translating this into Python code you obtain:

import numpy as np
import matplotlib.pyplot as plt

fs = 5e5       # Sampling frequency [Hz]
duration = 10  # Duration of the test signal [s]

N = int(fs * duration)
t = np.arange(N, dtype=float) / fs

f = 1000  # Frequency of the test signal [Hz]
sensitivity = 0.004  # Sensitivity of the microphone [mV/Pa]
p_ref = 2e-5 # Reference acoustic pressure [Pa]
amp = sensitivity * np.sqrt(2)  # Amplitude of sinusoidal signal with RMS of 4 mV (94 dB SPL)
signal = amp * np.sin(2 * np.pi * f * duration * t)  # Signal [V]

# Calculate the level from time domain signal
rms_time = np.sqrt(np.mean(signal**2))
db_time = 20 * np.log10(rms_time / sensitivity / p_ref)

# Apply window to the signal
win = np.hamming(N)
signal = signal * win

# Get the spectrum and shift it so that DC is in the middle
spectrum = np.fft.fftshift( np.fft.fft(signal) )
freq = np.fft.fftshift( np.fft.fftfreq(N, 1 / fs) )

# Take only the positive frequencies
spectrum = spectrum[N//2:]
freq = freq[N//2:]

# Since we just removed the energy in negative frequencies, account for that
spectrum *= 2
# If there is even number of samples, do not normalize the Nyquist bin
if N % 2 == 0:
    spectrum[-1] /= 2

# Scale the magnitude of FFT by window energy
spectrum_mag = np.abs(spectrum) / np.sum(win)

# To obtain RMS values, divide by sqrt(2)
spectrum_rms = spectrum_mag / np.sqrt(2)
# Do not scale the DC component
spectrum_rms[0] *= np.sqrt(2) / 2

# Convert to decibel scale
spectrum_db = 20 * np.log10(spectrum_rms / sensitivity / p_ref)

# Compare the outputs
print("Difference in levels: {} dB".format(db_time - spectrum_db.max()))

plt.semilogx(freq, spectrum_db)
plt.xlim( (1, fs/2)  )
plt.ylim( (0, 120) )
plt.grid('on')
plt.xlabel("Frequency [Hz]")
plt.ylabel("SPL [dB]")

plt.show()

  Difference in levels: -1.4210854715202004e-14 dB

As you can see, this is producing consistent results for both time and frequency domain. Good luck!


The shortest answer to that question is to divide your PSD by the squared sum of window values $w_i$:

$$PSD(n) = \dfrac{2|X(n)|^2}{S^2} $$

where:

$n=0\ldots N/2 $

$S = \sum_{i=0}^{N-1}w_i$

You can read more on that topic in the following article:

  Heinzel G., et al - Spectrum and spectral density estimation by the
  DFT, including a comprehensive list of window functions and some new
  flat-top windows.


The variance is a good inhomogenity measure.
For example, take 20*20 patches and calculate the empirical variance of the patch using these 400 samples.
For very simple purpose, this could work.

Edit : After viewing your example images, try to just play around with simple tools like manual thresolding and edge detection.

Edit 2: Added Example with different patch sizes (40, 10, 2)


How to make sinewives with array input, for-loop functions and operate array with others [matlab]?

I really don't like doing someone's homework so I'll drop the loops and make a few other modifications and leave it up to you figure it out.

I don't really understand why you filtered the input but I kept it in:

clc
clear all
close all
%%
Ts = .001;
N = 1000;
t = (0:999)*Ts;
%f = [2 10 15 20 30 40 50 60 80 90 100 150 200 300 400];
% f = make_an_array;
f=logspace(log10(2),log10(400),50);
%%
input=exp(1j*2*pi*f'*t);
[b,a] = butter(4,.02,'High');
[b1,a1] = butter(2,0.12);
input = filter(b1,a1,input);
output = filter(b,a,input);
G=output./input;
figure(1)
semilogx(f,20*log10(max(abs(G'))))
xlabel('frequency in Hz')
ylabel('Gain dB')
%A = max(output);
grid on
%disp(A)


Matched filter for variable time scaling?

One first quick and dirty method could be to first store a sort of scale-space library of your original waveform $w(t)$, computed  at different scales (discretized), and then perform a "filter bank" of matched filters in parallel on the same observed signal. This proposition is a kind of  dual of the non-feasible solution you evoked. This way, you can save time from having the waveform scaled offline, and do not have to it on the signal. Then, you would probably have to detect the best match at a given time. Ambiguities may appear by superposition: at a given time $t_0$, the signal could look like $a w(t-t_0) + b w\left(\frac{t-t_0}{s}\right)$, where $s$ is a scale parameter.

To see how much better could an offline solution be, you may have a look at a scale-invariant matched filter. Thanks to you question, I just discovered the recent paper:  The Mellin Matched Filter, A. Monakov, 2015:

  In propagation channels signals are undergone translation and dilation
  changes. The Fourier and Mellin transforms are natural foundations of
  the analysis of wideband signals. Three interconnected problems, where
  the Mellin transform plays a key role, are considered in the paper: i)
  synthesis of a linear matched filter that is invariant to signal
  scale; ii) estimation of the signal scale; and iii) the wideband
  ambiguity function and its properties. The filter that is invariant to
  the signal scale and maximizes the output SNR is synthesized in the
  paper using the Mellin transform. The filter is named the Mellin
  matched filter. Introduction of the Mellin matched filter allows to
  give a new interpretation of the signal scale estimation and the
  wideband ambiguity function. Concept of the Mellin matched filter
  facilitates to produce a physical interpretation of the resolution of
  wideband signals and the time-frequency-scale uncertainties.

I have no idea about how it performs in practice.


transient response of simple digitalized RC low-pass filter

If you assume a zero initial condition (i.e., if the system is switched on at $n=0$, then $y[-1]=0$), then the system is fully described by its impulse response, i.e. its response to the unit impulse $x[n]=\delta[n]$ (which is zero everywhere, except at $n=0$, where it is $1$).

Write the system equation as

$$y[n]=b_0x[n]+b_1x[n-1]-a_1y[n-1]\tag{1}$$

In your example $b_0=b_1$. If you apply $x[n]=\delta[n]$, you get the following response:

$$\begin{align}y[0]&=b_0\\y[1]&=b_1-a_1b_0\\
y[2]&=-a_1(b_1-a_1b_0)\\
y[3]&=a_1^2(b_1-a_1b_0)\\\vdots\end{align}$$

So we get for the impulse response

$$h[n]=\begin{cases}0,&n<0\\b_0,&n=0\\
(-a_1)^{n-1}(b_1-a_1b_0),&n>0\end{cases}\tag{2}$$

The general input-output relation is described by the convolution of the input signal and the impulse response:

$$y[n]=\sum_{k=0}^{\infty}h[k]x[n-k]\tag{3}$$

E.g., the step response can be easily computed from the impulse response:

$$a[n]=\sum_{k=0}^{n}h[k]\tag{4}$$


How to ensure anti-aliasing in image resampling at lower resolution?

A "1D image" is just a discrete signal. All the classic signal theory applies:

  Will such transformation introduce aliasing (and if so, what is the correct way to resample with averaging over pixels?)

So, what happens when you downsample is decimation. The aliases are the "folding" of existing spectrum over the conserved part of the spectrum; imagine you have a signal whose spectrum (for simplicity, assume that's its discrete Fourier transform) is

$[A\,|\,B\,|\,C\,|\,D\,|\,E\,|\,F\,|\,G\,|\,H]$

and you simply throw away half of the samples, the resulting spectrum is

$\frac 12 [A+E| B+F| C+G| D+H]$.

  Is there some math procedure, using which I can answer questions like #1 myself?

So what you'd need to do is suppress $[E\,|\,F\,|\,G\,|\,H]$ before throwing away the samples. Something that does that is usually called a low pass filter (LPF).

LPFs are often implemented as Finite Impulse Response Filter (FIR), where you take the current input sample, multiply it with a factor, add a previous input sample (also multiplied with a, possibly different, factor).

As it turns out, your averaging is indeed such an LPF FIR. It has the taps (coefficients) of

$[0.5, 0.5]$. To see this, consider your $2i+1$th sample to be your "current" one, and the $2i$th to be "one in the past".

Now, such a FIR can be analyzed; it's z-domain frequency response is this:

Notice how at the vertical bar, i.e. the point below which you want to preserve and above you want to suppress the spectrum, this only has a value of 0.5; in fact, about 22% of the total pixel "energy" that passes through your filter is in the "upper" spectral half. That's too much, and you'll get aliasing.

You simply might want to build a better filter. This implies you'll add more taps, i.e. more samples right and left of your current sample

There's plenty of tools out there to design filters. Scipy's filter design tools will absolutely do; just design a low pass filter.

For example, I've used GNU Radio's gr_filter_design tool, and specified the passband to end just below the half-band boundary (at frequency $0.45f_{sample}$), and the stopband to start just above that boundary ($0.55f_{sample}$)¹. Suppression should be 30dB, i.e. the magnitude in the $\le\frac1{1000}$ of the passband in the stopband. The resulting filter:

It has 13 taps, i.e. one coefficient for the "center" pixel and one for each of the six pixels left and right of it:

-0.00405147997662425, 0.009050002321600914, -0.02354922890663147, 0.04652633145451546, -0.07230162620544434, 0.09264583140611649, 0.9033603072166443, 0.09264583140611649, -0.07230162620544434, 0.04652633145451546, -0.02354922890663147, 0.009050002321600914, -0.00405147997662425

Notice the symmetry (typical for the class of linear-phase FIRs).

Original | Filtered with your averaging filter, then downsampled | filtered with the 13-tap filter, then downsampled (please view in full size; it doesn't make sense to compare things scaled by your browser)².

Visual comparison shows that the readability of the text suffers a bit more from the long filter response, while the reproduction of the fine structures of the dark cacao powder atop the tasty, juicy FFTiramisu is is better with the 13tap filter.
This is most probably an effect of the filter design (it's a design processed based on the Hamming window). Pre-scaling low pass filters employed by picture processing software often aren't even linear; that makes modeling significantly more complicated, but since sensory input in humans typically has logarithmic charactistics, I presume that hides "border effects" better.

  How, in general, discrete signals are anti-aliased? From DSP books I know that for continuous signals you just apply low-pass filter to your input signal and then sample whatever you got after the low-pass transform. But here is the problem: for images from the very beginning I have only discrete signal and I cannot figure out how to use continuous coordinate low-pass filters for discrete signals. May be I should treat an image as a piece-wise continuous function, which is constant within each interval [i,i+1)[i,i+1) ?

Well, in DSP downsampling is a very common thing to do, so you do with digital filters, as explained above. Don't consider your signal to be continuous; it's not.

¹ FIRs can't have infinitely steep responses. As a general rule, the closer the passband edge is to the stopband edge, the more taps you will need to realize the filter.

² Processing was done using Python, numpy, scipy.signal.lfilter, scipy.misc.im{read,save}.


The frequency response function (FRF) fails to detect the antiresonance of a system

I agree with A_A that the chirp is too fast.

I adjusted a few things in the code

increase time to $600$ (seconds)
start chirp at $0.5\textrm{ Hz}$
end chirp at $2\textrm{ Hz}$
removed zero padding (not required here)

And the result is...

The full code is

% Identification of the Frequency Response of a Transfer Function with
% Resonance and Antiresonance peaks
%% Cleaning the house
clc;
clear;
close all;
%% Definition of the System
wn_z1 = 10*2*pi;           % Anti-Resonance Natural frequency 10 Hz
wn_p1 = 80*2*pi;           % Resonance Natural frequency 80 Hz
chi_z = 0;                 % Damping factor of the zeros
chi_p = 0;                % Damping factor of the poles
s     = tf('s');
sys_1 = (s^2+2*chi_z*wn_z1*s+ wn_z1^2)/(s^2+2*chi_p*wn_p1*s+wn_p1^2);
zpk(sys_1)
%% Definition of the Chirp Input
h                    = 0.001;
time                 = 0:0.001:600;
freq_ini             = 0.5*2*pi; % Initial frequency considered (0.1 Hz)
freq_final           = 2*2*pi; % Final frequency considered (240 Hz)
chirp_input          = 1.5*chirp(time,freq_ini,time(end),freq_final,'logarithmic');

frequency_evolution  = logspace(log10(freq_ini),log10(freq_final),length(time));
figure(1);
subplot(2,1,1);plot(time,chirp_input); xlabel('time (s)'); ylabel('Chirp    amplitude')
subplot(2,1,2);plot(time,frequency_evolution/(2*pi)); xlabel('time (s)'); ylabel('Frequency (Hz)');
%% Execution of the simulation
[sys_output,~]       = lsim(sys_1,chirp_input,time);
%% Computation of the FRF
fft_input   = fft(chirp_input);                                           %% Considering zero padding
fft_output  = fft(sys_output );                                           %% Considering zero padding
fft_input   = fft_input(1:length(fft_input)/2);
fft_output  = fft_output(1:length(fft_output)/2);
fft_output  = fft_output(:);
fft_input   = fft_input(:);
FRF              = abs(fft_output)./abs(fft_input);
frequency_vector = linspace(0,1/(2*h),length(FRF));

%% Comparing the Bode Plot (Matlab) with the FRF (own implementation)
opts = bodeoptions('cstprefs');
opts.FreqUnits = 'Hz';
close all;
bodemag(sys_1,opts,'r'); hold on;
semilogx(frequency_vector, 20*log10(FRF ),'b')
xlim([0.1,frequency_vector(end)])
legend('Matlab Bode','FRF Identification')


Is $A\cos(\omega t+\theta)$ a Gaussian random process?

I'm pretty sure that even algebraically, $Z$ will not be Gaussian.  You're multiplying a Gaussian by a random variable that has a distribution that is effectively the histogram of a (co)sine function.

Let's do it empirically, though, following the example here.

Looking first at the densities, you can see that the density estimate for $Z$ (red) is much more central than that of $A$ (green).

Plotting the qqnorm and qqline plots for both, are below. As you can see, the one for $Z$ does not coincide much with the red line.

First, for $A$:

Then, for $Z$:

Regarding stationarity: the signal $Z$ is not stationary at all because it is cyclostationary.

R Code Below

# 28338
T <- 10000
sigma  <- 1

A <- rnorm(T,0,sigma)
theta <- runif(1,0, 2 * pi)

fs <- 100
t <- (1:T)/fs
omega <- 2* pi * 0.0892987
Z <- A * cos(omega*t + theta)

# First diagram
hist(sin(omega*t))

# Second diagram
dA <- density(A)
dZ <- density(Z)
plot(dZ, col="red", lwd=5)
lines(dA, col="green", lwd=5)

# Third diagram
qqnorm(A)
qqline(A, col = 2)

# Fourth diagram
qqnorm(Z)
qqline(Z, col = 2)


Output of LTI (in time and frequency $\omega$ domain) : when input goes through LPF

Since this appears to be homework, I can give you a few generic pointers

LTI systems maintain periodicity: if the input is periodic so is the output and vice versa.
If a signal is finite in length it can't be periodic
If a signal has finite bandwidth, it must be infinite in time.
You can model the input as $x(t)$ as cosine multiplied with a rectangular function.
Multiplication in the time domain is convolution in the frequency domain.


Is there a difference in bicubic interpolation of images vs patches?

The short answer is no. The order of operations doesn't matter and the final result would be similar within reasonable approximation errors.

The long answer is:

A given image $I$ has a pixel resolution of $M \times N$ pixels and a pixel size of $P \times Q$mm. Therefore, the image physical size is $M \cdot P \times N \cdot Q$mm.

If you downsample that (with any interpolation method) by some factor $d \in \left(0 \ldots 1\right)$, then its new physical size will be $d \cdot M \cdot P \times d \cdot N \cdot Q$. Obviously, here $d \cdot M \le 32$ and $d \cdot N \le 32$.

By the way, since at this point we are talking about reduction in physical size, this means that some of the original pixels, some of the original actual information content, had to be thrown away

If you then select the $32 \times 32$ patch, you are selecting some "space" out of the "shrunk" image, some of which information content has now been thrown away. ($32 \cdot d \cdot M \cdot P$)

If you were to select a $32 \times 32$ patch out of the original $M \cdot P \times N \cdot Q$ image, then you select a patch at full resolution.

Once you have selected the patch, the rest of the image is "thrown away". Any downsampling of that patch will produce a new patch whose dimensions (and actual information content) will be lower than $32 \times 32$. If you are thinking at this point "Oh, its alright, I will cut out a $32 \times 32$ patch, downsample it to (for instance) $12 \times 12$ and then upsample that patch back to $32 \times 32$" then, in terms of actual information content you are not gaining anything and you are back to the first case.

In fact, practically, because of the double interpolation (one to downsample and then one to upsample), this second case will have more opportunity for numerical errors but in terms of actual information content, actual spatial resolution, you don't gain anything.

The "dangerous" bit here is the downsampling. Once that takes place, that is where actual information is thrown away. This is why, when you are trying to resample a signal by a non-integer factor $\frac{u}{v}$, you first upsample by $u$ (where no new information content is added but the same information is distributed over larger space) and then downsample by $v$ (where actual information get's limited).

In your case, you have two different routes that both include downsampling by the same factor at some point. Their end result will be similar within reasonable approximations afforded by the interpolation steps.

Hope this helps.


Relation between signal processing and control systems engineering?

There is a lot of overlap but some differences in emphasis. Control Engineering is also older than DSP.  If you have a traditional EE education, you don’t really make much of a distinction.

State variables are the more typical perspective in Controls. The first edition of Oppenheim and Schafer 1975, had a chapter on state variables, but they dropped it over the years.  You need to understand state variables to do Kalman Filtering which is an area of overlap. Linear Estimation and Linear Controls are duals of each-other.

I would also say that hybrid continuous/discrete time systems are more common in Controls but there are many examples for DSP as well.

DSP is almost always done on uniform sampling. State Variables can work with nonuniform sampling as well.

I’ve never heard of anti causal Control System but forward backward filtering in time is common in DSP. Controls are inherently causal. The one sided Laplace transform is more common in controls.

Stability in feed back loops is important in both areas.  An advanced control systems class will cover topics like Lyaponov stability.  You typically don’t see that covered in DSP but there are DSP papers that use that technique.

Control Theory shows up in mechanical engineering. DSP shows up in finance. There is a-lot of both in robotics which also uses computer vision.

In RADAR, waveforms and filtering are more DSP at the front end, but the tracking systems at the back end are more Controls like.

If I had to use a single word to describe each.

Controls: feed back

Signal Processing: sensing

or maybe using a phrase

Controls: in-the-present

DSP: in-the-groove


FM demodulation amplifier with saturation

The receiver's output does not depend on the input amplitude.  From the point of view of the receiver, the output only depends on the frequency/phase of the input.

An fm signal in an ideal environment has a constant amplitude.  Another way of putting it is that fm modulation has a constant envelope.  Often this is advantageous because it allows the transmit amplifier to run at saturated output and not impact the modulation since amplitude is not varied as opposed to QAM and other non-constant envelope modulation types.  The same should hold true for the receiving side of the equation.  In practice I imagine you'd still need to watch out for overloading the receiver and causing harmonics and other nasties.

A write-up on my pyFmRadio blog post illustrates FM's amplitude independence from the modulating audio source.


What is the delta parameter in $\tt scal2freq$?

The delta argument in the scal2freq MATLAB function is the sampling period that relates to the data sampling frequency (i.e. delta = 1/Fs).  A good description of this function can be found here: http://au.mathworks.com/help/wavelet/ref/scal2frq.html it also includes an example including how to calculate and define the delta argument.  I hope this helps!


I think all it's saying is that we can write $H(\Omega)$ as:

$$
H(\Omega) = \left| H(\Omega) \right | e^{j \angle H(\Omega) }
$$

It also seems to be assuming that $h(t)$ is real-valued so that $H(-\Omega) = H^*(\Omega)$.

The function $H(\Omega)$ is just a map from $\mathbf{R} \mapsto \mathbf{C}$ (from the reals to the complex numbers).  For any specific value of $\Omega$ this generates a single complex number $H(\Omega)$. And, as you say, any complex number resolves to its exponential notation.


$l_2/l_2$ guarantee on sparse Fourier transform

A signal $s$ approximation from $x$ is often quantified by   the closeness of the approximant $s'$ and the desired behavior. Traditionally, one uses norms ($\ell_p$, $p\ge 1$),  pseudo-norms ($0<p < 1$) or the so-called $\ell_0$ count index. A very traditional formula uses a linear transformation $\phi$, and the approximation operator $A$, $s'=A(\phi x))$:

$$ \ell_q(s -A(\phi x))\,.$$

Minimization is often looked for in specific signal classes.
Looking for a $k$-sparse approximation is notoriously complicated: satisfying that $s'$ has at most $k$ nonzero components is generally NP-hard (B. K. Natarajan, Sparse Approximate Solutions to Linear Systems, 1995). So one may want a faster  approximant to be bounded by a measure based on the error with the best $k$-sparse  solution:

$$ \ell_q(s -s')\ \le C\min_{s_k} \ell_p(s -s_k)\,,$$

where $C$ is a constant as close to $1$ as possible. Sometimes, this error is guaranteed with some probability $p$. One often finds the terms "mixed norm error guarantee".

Many references can be found in the compressed sensing literature. On of the early source is A. Cohen, W. Dahmen, and R. DeVore. Compressed Sensing and best k-term approximation, 2006. Other references include $\ell_2/\ell_2$-for each sparse recovery with low risk, A. C. Gilbert et al.  You can find an $\ell_2/\ell_1$ in
  One sketch for all: Fast algorithms for Compressed Sensing, A. C. Gilbert, M. J. Strauss, J. A. Tropp.

A nice survey could be Sparse Recovery Using Sparse Matrices, Anna Gilbert, Piotr Indyk.


Orthogonality of filter impulse response to its even shift

Since this is clearly homework, let's help you with hints:

The set of  tricks for such problems is typically pretty limited; try (and that's really what solves the issue)

Multiplying the equation with elements from your right-hand term from left and/or right, so that identity matrices form.
Writing your matrix/operator property down as structure of your matrix
applying transposes, and hermitians until things happen.


filtering difference b/w freq and time domain

Multiplication in the frequency domain is the equivalent of convolution (digital filtering) in the time domain - one of the fundamental properties of the Fourier transform.
Using an FFT in practice to filter requires stitching together overlapping chunks and a large delay if in real-time.  However in the frequency domain you can immplement brick-wall filters as easy as any, which is much more costly in time-domain.
The algorithms for fft convolution include "overlap-add" and "overlap-save", these are useful search terms on the matter.
The Python scipy.signal library has an implementation of convolve that automatically selects direct filtering or fft convolution based on an estimate of which will be faster.


First Order Hold discrete-time approximation to first order continuous-time linear system

The solution comes from assuming that
$$
x(t) = x[k] + \frac{x[k+1]-x[k]}{T}(t-kT), \ \ t\in[kT,(k+1)T)
$$
where $x[k]:=x(kT)$, this is, a first order hold or a linear interpolation of $x[k]$ and $x[k+1]$ in between $kT$ and $(k+1)T$. Now, we simply compute the EXPLICIT solution of the differential equation (at instants $t=kT$) by integrating it from $t=kT$ towards $t=(k+1)T$. This should be enough to derive the result but let me get into more detail. First, multiply the ODE by $\frac{1}{\tau}e^{t/\tau}$ to obtain:
$$
\frac{dy}{dt}e^{t/\tau} + \frac{1}{\tau}e^{t/\tau}y(t) = \frac{1}{\tau}e^{t/\tau}x(t)
$$
and notice that $\frac{d}{dt}\left(e^{t/\tau}y(t)\right) = \frac{dy}{dt}e^{t/\tau} + \frac{1}{\tau}e^{t/\tau}y(t)$ and thus:
$$
\frac{d}{dt}\left(e^{t/\tau}y(t)\right) = \frac{1}{\tau}e^{t/\tau}x(t)
$$
Now integrate both sides:
$$
\int_{kT}^{(k+1)T}\frac{d}{dt}\left(e^{t/\tau}y(t)\right)dt = \int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt
$$
the left hand side results (by fundamental theorem of calculus) into:
$$
\int_{kT}^{(k+1)T}\frac{d}{dt}\left(e^{t/\tau}y(t)\right)dt = \left.e^{t/\tau}y(t)\right|_{kT}^{(k+1)T}  =e^{(k+1)T/\tau}y[k+1] - e^{kT/\tau}y[k]
$$
Thus,
$$
e^{(k+1)T/\tau}y[k+1] = e^{kT/\tau}y[k] + \int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt
$$
Now, simply divide by $e^{(k+1)T/\tau}$ to obtain
$$
y[k+1] = e^{-T/\tau}y[k] + e^{-(k+1)T/\tau}\int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt
$$
As you can see we are almost done. We just have to compute:
$$
\begin{aligned}
&e^{-(k+1)T/\tau}\int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}x(t)dt =\\ &e^{-(k+1)T/\tau}\int_{kT}^{(k+1)T}\frac{1}{\tau}e^{t/\tau}\left(x[k] + \frac{x[k+1]-x[k]}{T}(t-kT)\right)dt
\end{aligned}
$$
which we can do explicitly since $x[k]$ and $x[k+1]$ are just constants. After some algebra you can factor terms with $x[k]$ and terms with $x[k+1]$ and you will obtain precisely what you wanted. If you want to check, I plugged this into wolfram alpha here ("x1" is $x[k]$ and "x2" is $x[k+1]$).
This kind of discretizations are nice since, if $x(t)$ does behave as a piece-wise linear continuous function, the discretization is NOT an approximation, but an equivalent (and exact) representation of the ODE, this is: since it was obtained by explicit computation of the solution of the ODE, no information is lost about the output samples $y[k]$. (in contrast to approximating $dy/dt \approx (y[k+1]-y[k])/T$ or some other approximation as the bilinear one)
NEW REMARK: I just want to point out that this same procedure can be performed for other high order systems, or even for other types of holders. Recall that any LTI system with scalar input $x$ and scalar output $y$ can be represented by its transfer function, or equivalently by a state space representation
$$
\frac{dz}{dt} = Az + Bx,\ \ \  y = Cz
$$
with $z(t)\in\mathbb{R}^n$ with $n$ the order of the system, $A$ a square matrix, $B$ a column vector and $C$ a row vector. Recall that the poles of the transfer function are the eigenvalues of $A$ and zeros of the transfer function depend only on the coefficients on contained in $C$. Thus, we can solve for the system on $y[k]$ explicitly using a similar procedure as before:
$$
z[k+1] = e^{AT}z[k] + e^{A(k+1)T}\int_{kT}^{(k+1)T}e^{-A\tau}Bx(\tau)d\tau, \ \ y[k]=Cz[k]
$$
where $e^{A\tau}$ is the exponential matrix. Finally, plugging in the particular form of $x(t)$ results in the final form of the discretization (in discrete time state space representation). This can be used to derive (exact or "consistent") discretizations for ZOH, FOH, high order holders or more exotic holders.


Why not to simply use either:

imagesc(z);

Or:

h=pcolor(z)
set(h, 'LineStyle','None')

Although if you really want to use imshow then provide a set of extra parameters to scale the plotting range. It's more computational and memory efficient.

imshow(z, [min(min(z)) max(max(z))])


Common methods of transmitting digital data through an extremely noisy audio signal?

The scheme you are using is called On/Off Keying.  It is not terribly efficient, but it is simple and gets the job done.

When you say that the signal is 10 dB below the noise floor I suspect what you mean is that if you add up all the signal energy and all of the noise from 0.3 - 14 kHz the signal is 10 dB weaker, but that the signal uses a much narrower frequency range, making the signal stronger at very specific frequencies.  If it wasn't, you would not be getting reliable data transmission without doing long synchronous averaging, and it sounds like you aren't.

There are a number of different ways you could do this.  The state of the art would be to use an OFDM signal in the 17-21 kHz band.  A much simpler but not as efficient way to go (but still much more efficient than your OOK) would be to use a QPSK signal in the 17-21 kHz band.  If you have enough transmission power you could modify it to 8-PSK or 16-QAM for higher throughput.  Since it sounds like the transmission power is more of a limitation than the bandwidth (you can get a lot more data through 4 kHz of bandwidth than 80 bits/s), I would add a forward error correction (FEC) code to the data stream.

With QPSK you can get (ignoring the potential reduction in throughput from the error correction code) 1 bit/s/Hz, so you would have an upper limit of 4 kbits/s.  With 8-PSK it would be 6 kbits/s, and with 16-QAM it is 8 kbits/s.  The price of going with the higher data rates is that you are more likely to have errors.

The data rate for OFDM is more complicated.  Given your level of DSP knowledge I would go with QPSK, and then modify from there if needed.  It sounds like you would be thrilled with anything even close to 4 kbits/s, so that should be plenty good enough.


Determining the period of a discontinuous function

This will give you a plot of the autocovariance up to lag 100 samples:

plot(autocov(a,100))

There you can clearly get the period of your signal.

Another approach is to explicitly get the time of each pulse:

pulses = a > 0.1;
leading_edges = diff(pulses) > 0;
times = find(leading_edges > 0);
periods = diff(times)

ans =
    55   56   56   56   56

Explanation: pulses is 1 wherever there is a pulse. leading_edges gets a 1 only on the first sample of a pulse. times gives the sample number of each leading edge.

With this result you can perform your own statistics or estimation. This second method is computationally more efficient than using transforms. Computing time is linear with number of samples.


What is the meaning of a Nyquist diagram in a vibration signal

The "Nyquist plot" you have is really just the Fourier transform of the vibration data.  That the plot is "circular" doesn't mean much: it just means that the phase of the signal is continually increasing and that the magnitude is (roughly) constant.

I wouldn't call such a plot a Nyquist plot.  Nyquist plots are plots of the transfer function of a system. You may be able to use your measured signal and the reference to attempt to generate a Nyquist plot, but it's probably not going to be just a matter of multiplying / dividing Fourier transforms of the signals.

Generally, plotting the Fourier transform of a signal on the complex plane doesn't show much useful information.  More useful is to plot the magnitude against frequency, as this shows the dominant frequencies in the signal.

The reason the phase plots are "chaotic" (sic) is that the signal is noisy.  Taking the argument of a noisy complex signal will nearly always generate discontinuous phase plots.

I believe an Argand diagram is the same as your "Nyquist plot".


Ways to compute the n-the derivative of a discrete signal

The estimation of derivative is straightforward:

$$x'(n)~=\frac{x(n+1)-x(n-1)}{2}$$

$$x''(n)~={x(n+1)-2*x({n})+x({n-1})}$$

or if you have a signal sampled at $t_i=i\Delta t$,
it is

$$x'(t_{i})~=\frac{x(t_{i+1})-x(t_{i-1})}{2*\Delta t}$$

$$x''(t_{i})~=\frac{x(t_{i+1})-2*x(t_{i})+x(t_{i-1})}{(\Delta t)^2}$$

What you are interested in may be how to smooth the estimation. And yes you can use some recursive filters such as  $y(n) = a \cdot x(n) + (1-a) \cdot y(n-1)$, or just implement your estimation with some simple windows (Hann window for example).

To achieve the high SNR without distorting the signal very much, Savitzky–Golay filter also smooths your data by fitting successive sub-sets of adjacent data points with a low-degree polynomial with linear least squares.

EDIT

Matlab code for an N-th derivative of signal row vector x

dx = x; %'Zeroth' derivative

for n = 1:N % Apply iteratively

    dif = diff(dx,1); % First derivative
    first = [dif(1) dif];
    last = [dif dif(end)];

    dx = (first+last)/2;
end


Difference(relation?) between filter banks and wavelet decomposition

If we stick to the linear version and discrete versions of filter banks and wavelets, filter banks represent  the generic tool, and wavelets can be implemented as a specific instance of iterated $2$-band filter banks satisfying some additional properties, namely that low-pass spaces are embedded dyadically.

In other words: get a single-level $2$-band perfect reconstruction filter-bank, add some regularity so that iterations converge toward something legit, and you have a concrete Discrete Wavelet Transformation.

Kernel types have almost nothing to do with it, from a theoretical perspective. Kernels are chosen with respect to, mainly, implementation (fast or not) or signal properties (regularity).

That was the simple answer. If you want to dig deeper into it, you can look at continuous wavelets (that remain filter banks, as the admissibility condition impose a wavelet to be, somehow, a bandpass), nonlinear wavelets, morphological or hybrid filter banks, etc.

When I was a student, I was teached filter banks by Maurice Bellanger, one (the?) father of polyphase decomposition (Digital filtering by polyphase network: Application to sample-rate alteration and filter banks, 1976). As a young person, I just heard about wavelets, and asked him about it. His answer was: "they are just some special filter banks".


Given a filter (coefficients) how can I make it minimum phase?

I've seen Julius' MATLAB code and I know what it does.

Essentially, given an LTI filter with impulse response, $h[n]$, and frequency response:

$$\begin{align}
 H(e^{j \omega}) &\triangleq \Big| H(e^{j \omega}) \Big| \, e^{j \phi(\omega) } \\
&= \sum\limits_{n=-\infty}^{\infty} h[n] \, e^{-j \omega n}
\end{align}$$

Then $\Big| H(e^{j \omega}) \Big| > 0$ is the magnitude response and $\phi(\omega)$ is the phase response.

To be a minimum-phase filter, the phase response (expressed in radians) is the negative of the Hilbert Transform of the natural log of the magnitude response.

The natural log of the magnitude response, $\log \big| H(e^{j \omega}) \big|$, is just like decibels (dB) but is expressed in a different dimensionless unit, the neper.  8.685889638 dB is equal to 1 neper.  Essentially the nepers is the real part of the complex natural log of the complex frequency response and the phase in radians is the imaginary part.

Just like radians is the mathematically natural unit for angle, so also are nepers the mathematically natural unit for relative change of magnitude.  A change of magnitude of 1% is nearly the same as 0.01 neper.  But, like dB, going up 0.1 neper followed by going down 0.1 neper will land you at exactly the original magnitude.  But going up 10% followed by going down 10% will land you at slightly less than the original magnitude.

So, the (unwrapped) minimum phase of a filter, given it's magnitude response is:

$$\begin{align}
 \phi(\omega) &= - \mathscr{H} \Big\{ \log \big| H(e^{j \omega}) \big| \Big\} \\
\\
 &= - \int\limits_{-\pi}^{\pi} \log \Big| H(e^{j \theta}) \Big| \, \Big( 2 \pi \tan \big(\tfrac{\omega-\theta}{2} \big) \Big)^{-1} \, \mathrm{d}\theta \\
\end{align}$$

Where $\mathscr{H} \big\{ \cdot \big\}$ is the Hilbert Transform.

If we were to evaluate that integral, we would need to do something called "take the Principal Value" to deal with the division by zero when $\theta = \omega$.  But we won't do the Hilbert transform that way.

Okay, so Step 2 is understanding that the Discrete-Time Fourier Transform (DTFT) of this sequence:

$$ g[n] \triangleq \begin{cases}
0 \qquad & n<0 \\
1 \qquad & n=0 \\
2 \qquad & n>0 \\
\end{cases} $$

is

$$\begin{align}
 G(e^{j\omega}) &= \sum\limits_{n=-\infty}^{\infty} g[n] \, e^{-j \omega n} \\
\\
 &= \frac{e^{j\omega}+1}{e^{j\omega}-1} \\
\\
 &= \frac{e^{j\omega/2}+e^{-j\omega/2}}{e^{j\omega/2}-e^{-j\omega/2}} \\
\\
 &= \frac{2\cos(\omega/2)}{2j\sin(\omega/2)} \\
\\
 &= -j \frac{1}{\tan(\omega/2)} \\
\end{align}$$

i'm running outa time.  i gotta return to this.


If the FM radio at the receiver mixes with a another frequency to bring down the carrier frequency, wouldn't that affect the message signal?

This is the fundamental thing that all DSP bases on:

Any band-limited signal is 100% accurately representable by a baseband signal of the same bandwidth, and vice versa.

So, no, this doesn't affect the information at all!

In fact, this is the central result of Nyquist's works on the theory of continuous-time signals (with a few additions, it's equivalent to his sampling theorem).

  In today's digital world, why the FM radio still uses frequency modulation (FM) and not FSK?

Because of roughly 80 years of legacy devices! FSK is something that has a discrete alphabet of symbols, so it's only applicable to digital transmissions, really. Audio, as is, is an analog signal. So you'd first be converting your analog audio to digital, then modulating it, and then converting it back.

Fact-provenly, the spectral efficiency of digital modes is higher than that of analog FM – but you need digital receivers, which somehow still aren't that widespread.

That's a really long, and partly sad, story of badly introduced standards – there's multiple standards that strive to replace the analog radio broadcasts, and as far as I know, none has so far reached the proliferation that was once envision (note that this applies to terrestrial broadcasts only. For satellite radio, digital modes are overwhelmingly more popular):

DAB / DAB+: High-quality FM broadcast audio alternative; practically unused by anyone who isn't sitting in a high-end car with the "all-incluse" entertainment system. The difference between being able to receive DAB or not usually is something that the car dealerships earn a couple hundred Euros with. So not attractive at all, although German government has spent a lot of money on pushing that standard. Part of the problem was that early DAB modes had weak error correction, making the reception often worse than what people were used to via FM, receivers are still a lot more expensive, and DAB+, although it solves the quality issue, hasn't been adopted by everyone selling digital receivers. Coverage is mediocre, at best, compared to traditional FM worldwide.
DRM / DRM+: Digital Radio Mondiale, meant to replace the AM broadcast. Actually stands way better chances at seeing widespread adoption, because it solves actual problems, and the Indian Government alone has a 1MW transmitter to reach rural areas (read: most of India's area), so there's a lot of potential audience. Also possible candidate for broadcasting in Australia. There's about one receiver you can buy in Europe for that (source: been sitting in the same lab as the guys implementing the only FOSS transmitter in Europe).
add-on standards to "normal" FM: Just as stereo is an add-on to classical mono FM using more bandwidth, there are digital sideband and in-band modem standards (HD Radio etc) that provide additional functionality, like fine-grained weather info etc. Pretty much stillbirths, because where FM SNR is good enough, chances are high that people just get 2G/3G/4G coverage, offering actual advantages over just FM RDS.

By the way, none (to my knowledge) of these standards use FSK – if you're already doing digital radio in the 21st century, you'd typically go for something with higher spectral efficiency (most of these standards are OFDM things with PSK or QAM subcarriers).

So, to conclude, reasons for not adopting digital radio are:

literally billions of legacy devices that would need to be replaced
lack of advantage of digital modes ("FM sounds just fine, and I don't need more channels") or need ("99% of the population can get FM radio, and listener numbers are declining, why do we need new technologies")
actual technical disadvantages (e.g. DAB (not DAB+) was often worse than FM reception, and power consumption is still higher)
cost of new receivers high
bad time to market (broadcast standards take much longer to adapt than cellular standards, which now make broadcasting obsolete in metropolitan areas)
segmentation of market


How to Quantify a Difference in Spectral Information of two Signals?

Quantifying a difference in the frequency domain can be useful if the information on how your signals differ is better expressed in the Fourier space. If this assertion is valid in your case, hereafter are hints on methodology, and basically the sort of you can find in Shazam to distinguish music from a short listening period.
By the way, I think that I have been trying to find an answer to a related question  in 50% of my research hours pertaining to spectral or harmonic  analysis (the rest was coding). I am not fully successful yet.  I would however rephrase it as:

How could I easily find quantitative or selective (Fourier) features for
relatively stationary signals?

By quantitative I mean that the feature set should be able to measure scalar values related (in a monotonic way) to a quantity I seek in the original signal. By separative I mean that features that could tell whether signals are "apart". Those two concepts could be related by proper distances or metrics between the feature sets of two signals.
In your proposal, there is a notion of location in the frequency axis, a sort of median of spectrum amplitudes : 50% of the area under the curve for the FFT. A notion of amplitude is skipped, as you may normalize spectra. There is no notion of uncertainty or spread either. I am often trying to consider, as first, feature  triplets of Location-Amplitude-Spread ($f^k_0,f^k_1,f^k_2$) to describe the morphology of signals, spectra, or other suitable sparsifying transformations on signals.

Those three parameters could be compared to higher-order cumulants or  moments of a distributions. They are related to integrals of the shape:
$$M_d=∫^∞_{−∞}(t−c)df(t)dt$$
or their ratios. A related question is Finding the right measure to compare sound signals in the frequency domain.


I think your understanding is correct. The steady state response of an LTI system is the part of the response that is caused by a steady excitation at the input, like a sinusoidal signal, or any other periodic signal. The transient response is caused by changes at the input, like switching on a signal, or changing the parameters of a periodic input signal, such as amplitude, phase or frequency.
As an example of a marginally stable system, take an ideal integrator with impulse response $h(t)=u(t)$, where $u(t)$ is the unit step function. You can show that the response to a switched sinusoid at the input
$$x(t)=\sin(\omega_0t)u(t)\tag{1}$$
is
$$y(t)=\frac{1}{\omega_0}u(t)-\frac{1}{\omega_0}\cos(\omega_0t)u(t)\tag{2}$$
The first term on the right-hand side of $(2)$ is a non-decaying contribution caused by switching on the sinusoid at the input at $t=0$, and the second part is the steady state response to the sinusoidal input. The non-sinusoidal part of the output signal is what would be a transient if the system were a stable system, in which case it would decay. Since for this example it is non-decaying due to marginal stability, the term transient would of course be a misnomer.


Are there standards/references for determining spatial scale in photographic images?

One possible answer to my own question is provided by the photographic scale reference detailed on this website.
https://smallpond.ca/jim/scale/

It doesn't quite provide the solution I was looking for, buy may help spur ideas. I really need something that can accommodate a good scale estimate even if the reference is rotated in space (about any of the three axes), and is of low pixel resolution.


Calculating the inverse DTFT of a signal

You can develop it relatively easily, anyways, here's a screenshot from wikipedia


One method that works if there's a relatively strong drum beat is to take the magnitude of the STFT of the waveform, and then auto-correlate it in only the time dimension.  The peak of the auto-correlation function will be the beat, or a submultiple of it.

This is equivalent to breaking up the signal into a lot of different frequency bands, finding the amplitude envelope of each, autocorrelating each envelope, and then summing them.  The noise and other parts of music are averaged out by the cross-correlation operation.

This is because drum beats produce short-lived sound at many frequencies (vertical lines), while other parts of music are long-lived at only a few frequencies (horizontal lines), and noise is long-lived but random at all frequencies.  You can see the beat repetition if you look at an STFT:

I came up with this for a school project to find a single BPM value for entire music files, but it could be adapted to a stream of audio with changing BPM, too.  You'd need to process chunks that are at least twice as long as the period of the BPM you're looking for.


Segmentation of high noise, high density and low resolution image

There are a few things you can try but it still is not going to be an easy task. In general, in the absence of data you would be trying to exploit extraneous information by models or assumptions.

1. Try to make sense of the particles: The reason I asked about the particles was to understand if they are are symmetrical. If they are, then you could try a Gaussian filter adapted to the particle size, which would give you "high outputs" for anything that resembles a blob at that size. If they are not and they are "tumbling", then this would need to be taken into account. A single frame output filtered at 6 pixels looks like this:

Not very useful. The basic idea here would be to create a model of the particle so that you can somehow discriminate it amidst the noise. The Gaussian here is trying to match the particle but other features in the image look like the particle too (blob-like density variation). Which brings us to:

2. Try to get rid of the background: There are various ways you can do that. One of them is to average (or median) all of your frames into one. The median is more powerful in terms of discriminating relatively static parts of the image. One "master" frame obtained in this way looks like this:

Which is great because now you can see some more detail that is not exactly visible from the individual frames.

Now we try to subtract this background image from each frame of the sequence to get rid of the background and (hopefully) see remaining blobs flying around. Not exactly what was expected, but the sequence now looks like this:

Maybe I am not getting the timing right here but I don't see any particles flying around (?).

From this point onwards, you can combine approaches. For example, you can now apply the filter on the background-subtracted sequence of frames and then try to track the filter output (the blobs).

Or, you could try subtracting two successive frames and then try to track the gradient blobs. Assuming that the frame rate is sufficiently high (compared to the velocity of the particles), subtracting two successive frames will give you a gradient-like blob towards the direction of movement for each particle. This is expected to be consistent for moving particles and affected (to some extent) by noise. Which would then lead us to:

3. Try to track particles along specific routes: Assuming that you know approximately the routes on which the particles are supposed to be moving, you can create a model of that and use it, during your tracking, to decide if you are tracking a possibly "positive target" or some noise blob that by chance looked like a particle. But that now is way high up in the complexity scale and it assumes that you have a relatively good estimate of what a particle is.

Of course, all of the above would benefit from an increase in the amount of data you are capturing.

Hope this helps.


For simplicity of notation I'll leave out the index $k$. You've already noted that the error is

$$E=R\cdot G-S=S(HG-1)+NG\tag{1}$$

The MSE can be computed by integrating the power spectrum of the error. Since the power spectrum is real-valued and non-negative, its integral (i.e., the MSE) can be minimized by making the power spectrum as small as possible. Assuming white data and white noise with powers $\sigma_s^2$ and $\sigma_n^2$, respectively, which are independent of each other, we can compute the error power spectrum from (1) as

$$S_E=\sigma_s^2|HG-1|^2+\sigma_n^2|G|^2\tag{2}$$

The expression for the equalizer response $G$ that minimizes (2) can be found in an elegant way by completing the square in (2):

$$S_E=(\sigma_s^2|H|^2+\sigma_n^2)\left|G-\frac{H^*}{|H|^2+\frac{\sigma_n^2}{\sigma_s^2}}\right|^2+\frac{\sigma_n^2}{|H|^2+\frac{\sigma_n^2}{\sigma_s^2}}\tag{3}$$

The expression (3) can be minimized by choosing $G$ such that the left-hand term vanishes:

$$G_{opt}=\frac{H^*}{|H|^2+\frac{\sigma_n^2}{\sigma_s^2}}\tag{4}$$

which is the optimal MMSE solution.

It is not immediately obvious that (2) and (3) are identical, but proving their equivalence is simple (yet tedious) by multiplying both expressions out and comparing them.


Identifying Bad Sectors on magnetic stripe, or noise in signal

Typically you wouldn't directly use time-domain samples of a signal to declare its integrity (or lack thereof). I'm not familiar with the magnetic-stripe specification, but in most systems, any encoded data is protected by some form of checksum. This checksum is usually located in the same stream of symbols as the data of interest. Many systems also include error-correction coding, to provide tolerance to some number of symbol errors while still being able to recover the encoded information successfully.

Once the signal has been demodulated and you have recovered the symbol stream, you would calculate the checksum of the recovered data and compare it to the checksum that you read from the medium itself. If they match, then you can assume that the data is correct. If they do not, then you can try re-reading the medium again, or after some number of tries, you might declare it as a bad or damaged region.

This all begs the question, though: what is your goal of marking these regions as damaged? What will you subsequently do with the damaged areas?


Your controller controls pitch and roll. However, in the graph you showed us we only see x,y and z.
Could you add the plot for pitch and roll? I'm not an expert in aircraft control, but I think you can have a constant pitch and roll and still move in the x,y,z plane.
I'll edit my answer as soon as you give more info.
Edit : You roll seems stable, you could try adding an integral component to your pitch roll controller.


Log-Polar DFT Based Scale-Invariant Image Registration

If you want something really robust, but that might be more computationally expensive, you might want to check out the algorithm I implemented here. It implements the paper, "Robust Image Registration Using Log-Polar Transform" (pdf). It also has the advantage of being rotation-invariant, in addition to translation and scale invariant. In my application (art), it was able to register even even similar-looking images, not just transformed versions of the same image.


In order to understand the effect of DFT length $N$ on the computed spectrum samples, you shall be aware of two things:

1- The interrelation between theoretical DTFT and practical DFT.

2- The windowing effect on the analysed samples

For the first relation, consider a causal, discrete-time signal $x[n]$ of finite length $L$, whose DTFT, $X(e^{j\omega})$ is:

$$  X(e^{j\omega}) = \sum_{n=0}^{L-1} x[n] e^{-j\omega n}$$

This is mainly a theoretical tool, it's a complex valued continuous-frequency function which can not be represented in a digital computer system which can only work with sampled data of finite bit length.

In order to solve this issue, an alternative definition of Fourier transform for discrete time signals have been proposed and called as the DFT (discrete Fourier transform) $X[k]$ given as:

$$  X[k] = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2\pi}{N} n k} ~~~,~~~ \text{for} ~~~ k=0,1,...,N-1$$

Where $N$ is the DFT-size and not necessaily the signal length (usually $N \geq L$ is taken). When $N=L$ it can be seen from visual inspection that:

$$ X[k] = X(e^{j \omega})|_{w = \frac{2\pi}{N} k}  =  X(e^{-j\frac{2\pi}{N} k})$$

That's to say, N-point DFT values $X[k]$ of finite length signal $x[n]$ are the uniform samples of the DTFT $X(e^{j\omega})$ of the signal $x[n]$ itself. Hence we maintain the view that DFT values are the samples of the DTFT evaluated at $w_k = \frac{2\pi}{N} k$ frequencies.

Now it's obvious that the more samples you evaluate (larger the $N$) on DTFT $X(e^{j\omega})$ the smoother representation $X[k]$ will have. In addition, the minimum number of DFT samples that could exactly represent DTFT is given by $N_{min} = L$. Furthermore at this point it's very critical to understand the fact that increasing $N$ alone and keeping $L$ fixed, does not provide a higher-resolution spectrum samples, it merely provides a smoother set of finely spaced samples of fixed DTFT.

The second relation is what gives the clues of expected spectral resolution of the DFT samples. Now from multiplication by a window, $w[n]$ of length $M$, property of the DTFT, it can be seen that a finite-length excerpt $v[n] = w[n] x[n]$ from a signal $x[n]$ of indefinite (or very long) length is providing you the DTFT spectrum of a limited-resolution version of the true spectrum $X(e^{j\omega})$ which is based on all the samples of $x[n]$

The limitation of the spectrum is fundamentally determined by the window length $M$ and window type. So in a practical setting, your FFT block length (or input buffer length etc) gives you the windowing effect. The effect can be roughly summarized as a blurring of the true spectrum. The result of which is that close-by (and detailed) frequency content is mixed up and get lost in each other, so that you cannot recover them back.

So if you want highly detailed spectral examination of a signal, you should observe it long enough (your window length should be long) and you should use a long enough DFT size $N$ so that you get a dense sample representation of the true (blurred) spectrum $X(e^{j\omega})$ after it was windowed in time domain whose effect is a convolution with a low-pass character frequency function $W(e^{j\omega})$ resulting in a smeared true spectrum.

Then coming to your practical setting: when you have a signal of $x[n]$ sampled at $F_s = 44.1$ kHz, and rectangularly windowed by $M=1024$ samples, the computed DFT samples are providing the values of blurred spectrum whose rought frequency reolution is about $2 F_s /1024 \approx 86.1$  Hz. Using a Hamming window would instead provide a resolution of 160 Hz. And you should expect a resolution something in between. Therefore, eventhough the DFT bins would sample the windowed DTFT at spacing of $f_k = F_s/N = 43$ Hz apart, you cannot expect meaningful information at less than $86$ Hz steps at best and $160$ Hz steps most usually.


The filter order equals the number of taps minus $1$. The filter order is the order of the polynomial corresponding to the filter's transfer function. E.g.,

$$H(z)=h[0]+h[1]z^{-1}+h[2]z^{-2}\tag{1}$$

has $3$ coefficients (taps) but it is a second order polynomial (in $z^{-1}$), so the filter order equals $2$ (i.e., it has $2$ zeros and $2$ poles).

For type I and type III linear phase FIR filters, the filter order is even, i.e., they have an odd number of taps, whereas types II and IV have an odd filter order (even number of taps).


filterbank: understand the different responses at the center frequency of each filter

Normalizing the filterbanks by their widths is optional and totally up to you (similarly to the warping scale Mel/Bark). Depending on your application, you can start without normalization and see what results you are getting. Personally I prefer to keep it fixed and have one knob less for turning. There are more important parameters to tune, such as warping scale, number of banks and coefficients.

For example in the case of DNN's normalization is redundant, since you will be normalizing the input features anyway (at least that's what you should do), so any way you will end up with the same distribution.


Binary classification of grayscale image with little texture

Two large classes of techniques that may be of interest to you:

Mumford Shah functional based techniques.
Active Contour techniques.

Some googling on these terms should open some new doors for you.


The DC bin of an FFT should be the mean of all the other bins in the spectrum, why does the FFT produce 0?

Depending on your normalization factor the DC bin is either the sum of the signal values in your frame (1), the average (1/N), a "representation" (1/sqr(N)),

In the DFT formula the exponential value is $e^0$, or 1, for each term.


How to correctly use sinc interpolation in Matlab?

I remember being asked to do something like this in a DSP course and it was confusing because I wasn't sure how to model the continuous variable $t$ in MATLAB where everything is digital. You can still make a nice demonstration of sinc interpolation by choosing a low sample rate and use interpolation to create the interpolated signal which is at a higher rate.
Here is a simple MATLAB script that does what I'm talking about. Hopefully it at least helps you see how to use $x[n]$ to create the interpolated signal.
% Create "continuous time" signal, Fc >> f
Fc = 1e6; % very high sample rate to simulate "continuous time"
Tc = 1/Fc; % sampling period
t = (-0.1:Tc:0.1)'; % time axis
f = 10; % signal frequency
xc = sin(2*pi*f*t); % "continuous time" signal

% Create sampled signal
Fs = 100; % sampling rate
Ts = 1/Fs; % sampling period
ratio = round(Ts/Tc);
tn = t(1:ratio:end); % sampled time axis
xn = xc(1:ratio:end); % sampled signal

% Plot the CT signal and sampled signal
figure
hold on
grid on
plot(t, xc)
stem(tn, xn, 'o')
legend('"Continuous time signal"', 'Sampled signal')

% Create and plot sinc train
sincTrain = zeros(length(t), length(xn));
nind = 1;
figure
cmap = colormap(jet(length(-floor(length(xn)/2):floor(length(xn)/2))));
ax = axes('colororder', cmap);
hold on
grid on

plot(t, xc, 'k', 'LineWidth', 3)
for n = -floor(length(xn)/2):floor(length(xn)/2)
   sincTrain(:, nind) = xn(nind)*sinc((t - n*Ts)/Ts);
   p = plot(t, sincTrain(:, nind), 'LineWidth', 2);
   stem(tn(nind), xn(nind), 'Color', p.Color, 'LineWidth', 2)
   nind = nind + 1;
end
xlabel('t')
ylabel('Amplitude')
set(gca, 'FontSize', 20, 'LineWidth', 3, 'FontWeight', 'bold')

xr = sum(sincTrain, 2); % sum(sincTrain, 2) is the interpolated/reconstructed signal, should be equal to xc

figure
hold on
grid on
plot(t, xc)
plot(t, xr)

error = mean(abs(xc - xr).^2);


Does a filter add oscillations to a signal?

A linear filter can increase the magnitude of frequency spectra already present in a signal.  Note that the Fourier Transform of any (reasonable) waveform or signal (including ones that do not look at all like they are oscillating) will decompose it into sinusoids (which do oscillate).  Amplify any one frequency, and an oscillation may start to stand out and make it look like it just appeared; but a lower amplitude form of that spectra was already hidden in the input signal to start with.

It wasn't added.  It was amplified.  Which is possible with an LTI filter.


What is an intuitive explanation of cepstrum for a 12-year-old?

“efficiency” comes from a-priori assumptions.  One assumption that makes a Cepstrum useful is that the signal contains something roughly periodic including a set of harmonics or overtones, possibly  inside some sort of frequency spectrum envelope.

So a Cepstrum is mostly just an FFT of an FFT.  with a few mods.  On mod is to take the log magnitude after the first FFT to keep big overtones from drowning out the little ones.

So a DFT (slow FFT) correlates periodic happenings against some sinusoidal basis vectors.  Unless the periodic time domain  happenings are a single pure sinusoid, those periodic happenings will often correlate well against multiple equally spaced overtone frequency DFT basis vectors.

Thus you end up where you started, with a bunch of roughly equally spaced happenings, but now in the frequency domain instead of the time domain.  The second FFT will again correlate these against sinusoidal basis vectors related to that spectral periodicity, which for audio signals turns out to be an interesting proxy for (overtone rich) pitch.  The spectral envelope, since it was forced to be positive (by the log magnitude), ends up surrounding DC, instead of becoming sidebands of the spectral peaks as with time domain modulation.  This helps separate envelope information (formant, etc.) from pitch information in the Cepstrum result.

A Cepstrum also works for a single pure sinusoid, as the two spectral lines (positive and negative frequency peaks) correlate against a negative cosine of the same width as the line separation.

Note a required trick: One often has to make zero and tiny FFT spectrum magnitudes bigger before the log  to avoid the log from producing gigantic negative spikes or NaNs.


Are There Good Methods for Converting Color Image to Gray Scale Image?

There are method which are called Contrast Preserving Decolorization.
Those methods are built to keep the contrast in the Color image in the converted image.
They are based on Optimization Problem Solving which creates a Contrast Model in the color space and in the gray scale space.

For instance, have a look on that:


learnable segmentation or learnable edge detection

There is a benchmark code and test and learning dataset to perform image segmentation with ground truth dataset at Berkeley website and the code in the same place.


Your signal is periodic with period $\frac{1}{2f}$ hence has a Fourier series expansion with spacing between the frequency components of $2f$. The output of your ideal filter then is the sum of the Fourier components that fall within the bandwidth.

Because your bandwidth is less than the spacing between the frequency components your output will contain at most one frequency component.


References on the MP3 encoding algorithm

A compression standard is a quite delicate thing, that took years to develop and tune. I suggest Analysis of the MPEG-1 Layer III (MP3) Algorithm Using MATLAB, 2012 by Jayaraman J. Thiagarajan and Andreas Spanias. After a general overview, it leads you step by step though the different blocks of the whole scheme, providing some Matlab code for each of them.


Understanding 2D alignment using least squares

The set of feature points $\boldsymbol{x}_i$ and $\boldsymbol{x}^{\prime}_i$ is given, and you want to optimize the parameters $\boldsymbol{p}$ of a given planar transformation $\boldsymbol{f}(\boldsymbol{x};\boldsymbol{p})$ such that it predicts the actual transformation from $\boldsymbol{x}_i$ to $\boldsymbol{x}^{\prime}_i$ with minimum error. In this way, the parametric transformation optimally describes (and predicts) the actual transformation.

One way to do this is to minimize the error between the actual difference $\Delta\boldsymbol{x}_i=\boldsymbol{x}^{\prime}_i-\boldsymbol{x}_i$ and the difference predicted by the chosen transformation. The predicted difference is formulated in terms of the Jacobian of the transformation $\boldsymbol{f}(\boldsymbol{x};\boldsymbol{p})$ in order to get a linear least squares problem (which is easy to solve): ${\Delta\widehat{\boldsymbol{x}}_i}=\boldsymbol{J}(\boldsymbol{x}_i)\boldsymbol{p}$.

So the error to be minimized is given by

$$\sum_i\left[{\Delta\widehat{\boldsymbol{x}}_i}-\Delta\boldsymbol{x}_i \right]^2=\sum_i\left[\boldsymbol{J}(\boldsymbol{x}_i)\boldsymbol{p}-\Delta\boldsymbol{x}_i \right]^2\tag{1}$$

Chosing the parameters $\boldsymbol{p}$ that minimize $(1)$ minimizes the average squared error between the actual feature point difference $\Delta\boldsymbol{x}_i$ and the one predicted by the transformation $\boldsymbol{f}(\boldsymbol{x};\boldsymbol{p})$.


Implementing short-time Fourier transform

admittedly, i have never used the buffer() function in matlab.  for doing STFT analysis and processing, i have just broken it into frames myself with code similar to this answer.  it's easier to see what's going on.

while you can use a gaussian window for analysis (i found it mathematically useful in this paper, email me if you want a pdf copy), you must use a complementary window, one that adds to 1 with the windows of its adjacent frames, for the reconstruction of the signal.  at least if you want this to pass the signal unchanged when the processing is set to "null".


Matching features points of rotate image using SURF

SURF Features are not completely rotation invariant.

In the paper, they try to make SURF more robust to rotations by following the following methodology.

In order to improve the performance for rotation, the descriptors are based on multi-orientations. The many-to-many tentative correspondences are determined with a maximum distance. Hough transform is used to reject the mismatches and the affine parameters are computed with a least-squares solution.

The image matching algorithm shows a better performance for image rotation than the standard SURF and it succeeds in matching the image including repetitive patterns which will deteriorate the distinctiveness of feature descriptors.


Obtain motion between image features by means of the homography matrix

First, calibrate your intrinsics: The focal lengths and the principal point. Homography is not really a rigid transformation and rather a mapping of a plane onto another one. What you really need is a 6DOF parametrization, which is the camera pose. If you have the intrinsics, you could transfer your coordinates to the normalized coordinates and estimate the essential matrix using 8 point algorithm - or a better one would be five point algorithm (check this, this and this). The operation is the same as solving for fundamental matrix in the uncalibrated setting. But, just because you have your $K$, you could recover essential matrix. Note that, such estimation doesn't assume planar scenes and can in fact operate on 3D structures.

You could then decompose your essential matrix as

$$E=R[t]_x$$
where $[t]_x$ is the skew symmetric form. To obtain the components, proceed as described in here. Note that, this problem has 4 possible solutions and identifying the best one is usually found via checking all solutions.

If your cameras are uncalibrated, you could maybe approximate the calibration using some known dimensions or a grid, within the scene.

If you are using simulated images, then you know the FOV and image dimensions. If this is the case, you could get the camera matrix. Or even better, if you have your $K$ matrix (you could artificially make up one), then you could compute the OpenGL equivalent. Check here and here for that.

To address your questions one by one:

They are in pixel coordinates and denote the principal point (computed center of the image).
Unfortunately, $E$ is defined up to a scale and you could only recover the scale, if you have a known object with size within the scene. You just need to know one distance to do this. It could be a marker, the length of some existing furniture, or the object you inspect and etc.
Use the explanation I already provided (above).

And additionally, it is not an affine transformation, it's a projective mapping (perspective transformation). In fact homography is nothing more than the transformation of one plane to another.


To add to Famous Blue Raincoat's comment: the process of upsampling a signal aims to increase its sampling rate without affecting the signal otherwise. In your example, outputting x or resx through your sound card should produce exactly the same sound. Upsampling is useful when you need to match sampling frequencies between different systems without affecting the signals otherwise.

Likewise, downsampling aims to modify the signal as little as possible. All frequencies below the new Nyquist frequency will be left untouched.


The bad news is that I disagree with the first assertion: there exist non-orthogonal DWT (biorthogonal for instance).

But the good news is that this is not important here, and the answer is yes. A Laplacian pyramid is an overcomplete expansion. It can not be treated as a basis, but as a frame. And the concept akin to orthogonal bases in frames is the tight frame.

You can get details on such constructions in Framing Pyramids (2003), whose abstract is:

  In 1983, Burt and Adelson introduced the Laplacian pyramid (LP) as a
  multiresolution representation for images. We study the LP using the
  frame theory, and this reveals that the usual reconstruction is
  suboptimal. We show that the LP with orthogonal filters is a tight
  frame, and thus, the optimal linear reconstruction using the dual
  frame operator has a simple structure that is symmetric with the
  forward transform. In more general cases, we propose an efficient
  filterbank (FB) for the reconstruction of the LP using projection that
  leads to a proved improvement over the usual method in the presence of
  noise. Setting up the LP as an oversampled FB, we offer a complete
  parameterization of all synthesis FBs that provide perfect
  reconstruction for the LP. Finally, we consider the situation where
  the LP scheme is iterated and derive the continuous-domain frames
  associated with the LP.

And yes, the pseudo-inverse can be less simple, but you will gain more optimality in the presence of noise.

Notably, such design allowed the construction of contourlets, a multiresolution directional tight frame designed to efficiently approximate images made of smooth regions separated by smooth boundaries.


If your code is working correctly, than there should be no aliasing.

Most likely this is a problem with your playback system. Most sound cards in computers are terrible and will create a lot of noise and other artifacts. Since the 19 kHz tone is inaudible for most people, you just hear the artifacts.

It's also possible that you turn up the volume too high. The 19 kHz may easily clip your playback hardware. You can't tell if it's "too loud" since your auditory sensitivity  at 19 kHz is basically zero. You may have also clipped the sine wave when sending it to the sound card.

There is also the software stack: On Windows, most sound will go through the kernel mixer, which often applies sample rate conversion and other processing. It's default behavior is quite bad, so this is also likely to create a lot of artifacts and noise.

Playing sine waves cleanly requires decent equipment. Typically a good external sound card and a dedicated SW driver that bypasses the operating system. You can try saving it as a wave file and play it over some "known good" hardware.


Image processing: how do I erode on boundary and dilate inside boundary?

You will want to use the morphological closing operation, which is the erosion of the dilation of an image. This means that you first dilate the image, removing the small 'holes' in the interior of the image, and then you erode the image, shrinking the boundary. If you want the final boundary to be 'smaller' than the original boundary, simply erode more than you dilate, ensuring that you dilate enough to completely fill the interior holes first.

The idea behind this is that by removing the holes via dilation, the information about those holes is lost from the image. When you then proceed to erode the image, there are no longer interior gaps to erode away from, so it only erodes away from the outer image. Therefore, despite being opposing operations, the erosion does not simply reverse the dilation, provided the feature sizes to be removed are 'smaller' than the size of the dilation.


Why does this image look so bad when dithered with a certain palette?

I think what's going on is something analogous to CCD sensor blooming.

The badly working palette does not cover the extreme color values of the original image. Floyd–Steinberg dithering works by error diffusion, where the error in color at a pixel being processed is spread to nearby pixels not yet processed, to the right and to some of the pixels on the next pixel row. With the kitten, you have near-black and near-white areas, and in those areas the algorithm – working like a bulldozer – will accumulate more and more error because the colors it can use will always fall short the same way, and when it finally comes out of the problem area, the accumulated error is so large that it takes many pixels to get rid of it, smudging the image.

In other words, if the image has out-of-gamut colors, you get problems. If a color is inside any polyhedron formed by any of the palette colors, then it is in gamut and can be reproduced by alternating colors from the palette. In other words, the color is in the convex hull of the palette colors. If it's not in the convex hull of the palette colors, it is out of gamut. Gamut mapping before dithering would resolve the problem.


As pointed out in PeterK.'s comment, the way it is written this does not make sense. What is meant here is that the insertion loss can be measured on a linear scale (i.e., as a factor) as

$$L(\omega^2)=\frac{1}{|H(j\omega)|^2}$$

or it can be given in decibels:

$$L_{dB}(\omega^2))=10\log_{10}(L(\omega^2))$$

The remaining question could be why $|H(j\omega)|^2$ can be replaced by $H(j\omega)H(-j\omega)$. The reason is the symmetry of the frequency response $H(j\omega)$ for real-valued filters:

$$H(j\omega)=H^*(-j\omega)$$

From this it follows that

$$|H(j\omega)|^2=H(j\omega)H^*(j\omega)=H(j\omega)H(-j\omega)$$


But is there a better bound for a general filter response?

No. The argument about it $H$ being a minimum phase filter doesn't make a difference. In essence you are asking a math question on some bound for
$$L = \sum X^2 \cdot \sum \frac{1}{X^2}$$
There is none. Let's assume a very simple signal $x = [a,1]$ . The product of the norms is
$$L = (a^2+1)(1/a^2+1) = a^2 + 2 + 1/a^2.$$
It's clear that $L >=4$ and the equation is easily solvable for $a$ for any given value of $L$. In other words: by choosing $a$ you can make $L$ as large as you want it to be, so there is upper bound.


Does modified Morlet wavelet function satisfy the admissibility condition?

No, it doesn't (at least not in theory). You know that the admissibility condition is

$$\int_{-\infty}^{\infty}\frac{|\Psi(\omega)|^2}{|\omega|}d\omega<\infty$$

where $\Psi(\omega)$ is the Fourier transform of the mother wavelet $\psi(t)$. If $\Psi(\omega)$ decays sufficiently fast (which is the case for the Morlet wavelet) then the admissibility condition reduces to $\Psi(0)=0$, i.e. $\psi(t)$ must satisfy a zero-mean condition. However, this condition is not satisfied for the Morlet wavelet. So the Morlet wavelet is not admissible, but in practice this is no problem because $\Psi(0)$ is very small.

[1] M. Vetterli, J. Kovačević: Wavelets and Subband Coding


Fast cross correlation with limited range

If the limited range is significantly more than log(N) shorter than the total signal of length N, then ordinary cross correlation over that limited range may be faster than FFT-based fast correlation.  You could just iterate away from the estimated/predicted lag until you find you have gone past a correlation peak (by a sufficient amount to account for expected noise).  But if you go more than log(N)/2 points away from your initial estimate, then perhaps switch back to fast convolution cross correlation.


Feature Extraction of FFT for One Class SVM

Assuming you have epochs/segments of data. For this kind of signals, it is a safe approach to extract features using wavelet representations.

Using FFT might work as well, but I dont know how problematic would be the stationarity assumption in this kind of applications. Besides, FFT estimates for this kind of signals is sometimes very noisy. If you insist on this kind of representation, I would suggest estimating the power spectral density using Welch or multitaper methods. These usually deliver a nicer spectral estimation.

In any case, using any spectral features as direct input to your SVM might not be the best idea. Irrespectively of what features you extract (FFT, Wavelets, etc..) I would include a data dimensionality reduction step (PCA would be the first thing I'd try)


Confusion regarding Nyquist Sampling Theorem

It could be interesting to start from an history of this theorem in Interpolation and Sampling: E. T.  Whittaker, K. Ogura and Their Followers, by Butzer et al.

Putting history aside, the main thing with the basic theorem is that one should refrain to say that the rate  "should" or "has to" be above something. This is could a sufficient, but not a necessary condition. One version is, with an "if": if $X(f)=0$ for $|f| > B$,  sampling at a rate above $2B$ theoretically allows you to recover the signal from the  regularly sampled sequence:

There are milder versions, but this one shows that you ought to be especially careful if you choose, dangerously, $B$ such that $X(B)\ne 0$. For instance with your sine, whose spectrum is not $0$ at $B=f_0$. But is is only an if.

Some non bandlimited signals can still be sampled perfectly, under additional conditions. Some signals can be sampled at $B$ even if they spectrum does not vanish at $B$. And some signals can be sampled at a lower rate. Especially when the signal is bandlimited. @JasonR already pointed to Undersampling, and is dealt with in The theory of bandpass sampling, Vaughan et al., and already discussed in Minimum Sampling Rate of Bandpass signal.

I am not so confident with my understanding of this literature. However, in favorable circumstances with a real bandpass signal, twice the effective bandwidth can suffice, which can be much lower than twice the maximum frequency. And there are other theorems for signals whose spectrum is made of unions of segments of non null frequency.

While mentioning such results, I believe the real issue is different. In some texts, authors assimilate (wrongfully) what they call "the bandwidth" with the distance from the $0$ frequency to the maximum frequency. Especially in texts that do not go into fine details such as the ones stated above and in other answers.

So in your case, I assume that the difference between the two courses can be explained by laziness or mundane talking.

In practice though, you cannot rely solely on this theorem: it requires you to know in advance the spectrum of a continuous signal, which in many cases you have not idea about. Physical modeling  and analog filtering can help you, but in real life, finite precision quantization, signal jitter, noise and especially finite length signals (which cannot have bounded spectra) forster more precaution and a choice of a sampling frequency sufficiently  higher than twice the maximum frequency (or bandwidth if that applies).

A potential additional lecture: Sampling: What Nyquist Didn’t Say, and What to Do About It.


Yes. This looks like it's a typical first order low pass. The time constant can be determined by looking at the time it takes for the falling edge to drop to 37% of the max amplitude ($e^{-1}$). The cutoff frequency of the low pass filter is $1/(2 \cdot \pi \cdot t_c)$ where $t_c$ is the time constant


With respect to the dirac delta function, what does $\delta [x + 2A, y - A]$

A (two-dimensional) finite-impulse-response filter has an impulse response $g(x,y)$
that is nonzero only for $(x,y)$ in a region of finite area (continuous parameters)
of for a finite number of values of $x$ and $y$ (discrete parameter). For example,
$g(x,y)$ might be nonzero only for points inside a circle of radius $r$.  Specifically,
consider
$$g(x,y) = \begin{cases}1, &x^2 + y^2 \leq r^2,\\0, &\text{otherwise}\end{cases}$$
describes an image filter whose response to a point source (impulse function)
at the origin is a disc of radius $r$ in the filtered image.  For input $f(x,y) = \delta (x + 2A, y - A) - \delta (x - A, y + A)$, the filter output would be the sum of a positive
disc centered at $(-2A,A)$ and a negative disc centered at $(A,-A)$ with cancellation
occurring where the discs overlapped. (Disc ovrlap would be determined by the
relationship between $A$ and $r$). More formally,
$$g(x,y)\circledast f(x,y) = \begin{cases}
1, &(x+2A)^2 + (y-A)^2 \leq r^2, (x-A)^2+(y+A)^2 > r^2,\\
-1, &(x+2A)^2 + (y-A)^2 > r^2, (x-A)^2+(y+A)^2 \leq r^2,\\
0, &\text{otherwise.}\end{cases}$$
where we have combined the outer darkness of $0$ with the $0$ produced
by the cancellation of the discs in the region of overlap (if any).

Similar calculations can be done for the discrete case as well.


Conceptually, it's straightforward. Let's consider a more general case, involving a grid:

First run ridge detection, box blurring, then line detection:

The purpose of this preprocessing is to make line detection more successful. Finally, partition the image vertically (only first partition shown):

then horizontally:

In Mathematica it can be done as follows:

img = Import["http://englishwalle.files.wordpress.com/2013/08/9-image-grid-no-logo1.jpg"]
{W, H} = ImageDimensions[img]
edges = ImageConvolve[Binarize@RidgeFilter[img, 0.001], BoxMatrix[2]]
lines = ImageLines[edges, 2, 0.1]
hlines = Select[lines, Norm[#[[1]] - #[[2]]] == W &]
offsets = Flatten@{0, #[[1 , 2]] & /@ hlines, ImageDimensions[img][[2]]}
hparts = Partition[H - # & /@ offsets // Sort, 2, 1]
vsections = ImageTake[img, #] & /@ hparts

(* Part 2 *)
hmask = ImageAdjust[ImageAdjust@RidgeFilter[Blur[#, {100, 0}], 5], 1] & /@ vsections
lines = ImageLines[#, 0.075, 0.25] & /@ hmask
vlines = #[[1, 1]] & /@ # & /@ lines
vlinesUnique = W*DeleteDuplicates@MeanShift[#, 1] & /@ vlines/W
vparts = Partition[Flatten[{0, Sort@#, W}], 2, 1] & /@ vlines_uniq
hsections = Outer[ImageTake, {#[[1]]}, {{0, W}}, #[[2]], 1] & /@ Transpose[{vsections, vparts}]

I leave it to you to translate this into your language of choice.

edit: With your image, I had better results with a derivative filter than a ridge filter (though that worked too, but the seam was not exactly at the right place). This is what directional blurring and derivative filtering (hmask = DerivativeFilter[Blur[img, {100, 0}], {0, 1}]) looks like on your image:

Line detection works trivially after this.


There are rare situations when processing is done with integers (embedded systems with no floating point units), but usually, audio signal processing is done with floats.
(Handwaving) Convolution is the operation through which a filter is applied to a signal. If the authors of the paper you're looking at found that applying a specific filter to the data (for example to remove noise, any information not processed by the human auditory system, etc.) yielded good results, then you have to trust them on this!
Since audio signals are not stationary, they are processed in short blocks, referred to as "windows" or "frames". The analysis/processing is thus done on short overlapping segments. For example you will process segment 0.00s:0.02s ; then segment 0.01s:0.03s ; then 0.02s:0.04s etc. In this example, the window size is 20ms and the overlap 50%. There is no "universal" window size for audio applications, look at what is used in the paper describing the technique you're implementing.
Not really, but going from the time domain to frequency domain is "one step closer" to what is ultimately happening in the human auditory system.
No operation could do that. If there was one, we would send all our files with the lowest quality and apply this magical operation to recover the highest quality. What this means is that the features used for fingerprinting must be robust enough to capture only the overall "shape" of the sound rather than focus on details - which would be lost by an increased amount of MP3 encoding.
Look at a signal processing textbook. These are very basic signal processing operations and can be used for many things inside or outside the context of a fingerprinting system - reducing sample rate, smoothing temporal sequences, emphasizing details, removing noise...
Let's say you have 800000 samples. You apply the convolution, you get 800000 samples out of it. Now you extract every 8th sample, you get 100000 samples. Then you apply your FFT on samples 0..2047 to get your 33 coefficients ; then on samples 64..2111 etc... so in the end you have a list of roughly 100000 / 64 = 1562 vectors of 33 coefficients.
(Handwaving). I am assuming you are referring to the Ke et al system (later refined in the works of Baluja et al), which is itself inspired by the Viola-Jones face detector. Once you have converted the signal into a sequence of vectors of 33 coefficients, you compute weighted sums of the coefficients in adjacent bins or along the time axis (a filter), and you compare the result to a threshold. This gives you a binary "feature". There are 32 such binary features, each of them computed with a different filter and a different threshold. They are "learned", because they were obtained the following ways: the authors built a database with millions of pairs of songs, some of them identical (but maybe at difference bitrate, volume, noise level etc.), some of them different. Using machine learning techniques (boosting on decision stumps), they evaluated which filters were the best at binning together identical content and separating unrelated content, and which threshold were being used for that. You can find the list of those filters/thresholds in the source code of an open-source implementation of lastfm's Ke-style audio fingerprinter.


The problem you might have is that when bright lights are present, the camera will be using a fast shutter speed to keep the image from saturating too much.  All the detail you are interested in is in the bottom few bits of the data.

So even if you remove the glare and then try and pull the rest of the image up in level, it's very noisy.  For example, if you take your image and play with the levels to pull up the darker sections and just allow the glare to clip, this is the result:

If that's OK for onward processing (I don't know what you have in mind) then you could maybe remove the glare by looking at distributions of saturated pixels, and expand into the "glow" around them.

A better solution (if you have the option) can be to use a high-dynamic range camera which has an non-linear response in the pixels, and also often 10 or 12 bits per pixel of usable resolution, which means you can keep 7-8 bits for the darker sections.


It looks okay to me. If you define empirical mean $\hat{\mu}_d = \frac 1N \sum \tilde{d}[n]$ and empirical second moment $\hat{\gamma}_d = \frac 1N \sum \tilde{d}^2[n]$ then you effectively have an equation of the form $$-N/x + A/x^3 + B/x^2 \stackrel{!}{=} 0,$$ where $A$ and $B$ directly depend on $\hat{\mu}_d, \hat{\gamma}_d, \mu_A, \sigma_A$. If you multiply this equation with $x^3$ it becomes a quadratic equation in $x$, which you can solve for $x$...
I get something like $$x^2 - \frac{\hat{\mu}_d \mu_A}{\sigma_A^2} x + \frac{\gamma_d}{\sigma_A^2} = 0$$ but that was scribbled down in a rush so I would have to check it again later. Can you confirm or correct?


MATLAB IIRs -- extracting and scaling the coefficients of second order sections

These scale values are applied before and after each section, so for $N$ sections you have $N+1$ scale values. Check the Matlab documentation, there must be some related property that you can set to 'scalar', and this will give you a single total scaling value. In any case, the total scaling is just the product of the three values in your case. So you have to multiply the b coefficients of one filter section by the product of all scale values, then the total scaling will be OK, i.e. the maximum frequency response magnitude will be 0 dB.


Sampling interval $T$ as a multiplier in digital processing of a continuous-time signals

I don't think it's a typo, but it has to do with the sampling process. As shown in Eq. $(4.30)$ in the 3rd edition of Discrete-time Signal Processing by Oppenheim and Schafer, (ideal) sampling involves scaling by $1/T$, where $T$ is the sampling period. If the continuous-time signal $x_c(t)$ is ideally band-limited, the spectrum of the discrete-time signal $x_d[n]=x_c(nT)$ after sampling is given by
$$X(e^{j\Omega T})=\frac{1}{T}X_c(j\Omega)\tag{1}$$
So if you divide both sides of Eq. $(155)$ by $T$ you see that on the right-hand side you get the scaled spectrum of the continuous-time signal, which equals the spectrum of the discrete-time signal.


How does decimation work without losing samples?

Your understanding is correct. You need to make sure the information content fits into the output rate, and that implies your signal needs to be band-limited to 1/16384. Otherwise you'll lose information.
Making sure that's the case with a filter is called anti-alias filtering.


I'm going to use a grayscale image for simplicity. Suppose you have a 3D matrix and each slice is the result of filtering the image with the Gabor filter for that angle.

I = randn(256, 256); % Just some random image
G = gabor_filter; % Assume G is 256x256xNum_Angles, e.g. 256x256x6 if you had 6 angles
Y = zeros(256, 256, size(G, 3)); % Preallocate a temporary matrix

for idx = 1:size(G, 3)
    Y(:, :, idx) = conv2(I, G(:, :, idx), 'same');
end

Ymin = min(Y, [], 3); % This takes the minimum along the 3rd dimension


area under the phase encoding gradient in MRI

Actually, it is as simple as the frequency encoding gradient. Both encoding gradients look terribly different, but in the end they work the very same way.

The outer $k$-space position $k_\textrm{max}$ is related to your resolution, e.g. $1\,\textrm{mm}$, or in terms of the FOV: $FOV/N$ where $N$ is the number of image pixels.
For the steps between two $k$-space lines, $\Delta k = 1/FOV$

With $k(t) = \frac{\gamma}{2\pi}\int_0^t G(t^\prime)\,\textrm{d}t^\prime$ you can find that for the largest gradient moment $M_\textrm{max} = \frac{\pi(N-1)}{\gamma FOV}$. The steps in between scale this down linearly.

Pay attention to the line numbering, because the central line in $k$-space should be acquired with zero gradient moment.


"Does the interpolation of a signal alter the FFT?"
Yes, it does.
The key point is that in a system that uses sampling, you simply do not know what does the signal "do" between two sampling instances. For a high enough sampling frequency Fs, it is just assumed that the corresponding sampling period is high enough to be approximated by a series of sinc pulses. For more information please see: http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem

By using interpolation, you only assume that the signal behaves in a certain way during the time that it is not observed. Therefore, the result of the FFT will be different, even for different interpolation methods.

Furthermore, interpolation has a "natural" smoothing action because all it can do is "degrade" the existing data, it cannot add anything more to it. For more information about this notion, please see the Data Processing Theorem (https://books.google.co.uk/books?id=F4j8EHBKS1IC&lpg=PA313&ots=tHO2C1ZJs1&dq=the%20data%20processing%20theorem&pg=PA313#v=onepage&q&f=false)

This is why, interpolation will tend to represent lower frequencies (for which you have more data anyway, due to the sampling frequency) better.

By the way, 'interp1' by default uses linear interpolation, that is, intermediate values between sampling instances are assumed to lie on a straight line that is defined by the two known points. An alternative to this would be spline interpolation which takes into account the slope by which the interpolated curve should enter and depart from the interpolating area. For more information about this, please see http://en.wikipedia.org/wiki/Cubic_Hermite_spline

(Please also note Trigonometric Interpolation, a technique very closely associated with the Fourier Transform. For more information please see: http://en.wikipedia.org/wiki/Trigonometric_interpolation)

From the Time Signal, the FFT on the complet signal is closer to reality. but it still doesn't show all frequencies. Why?

Because it is not exactly the same signal, as far as the FT is concerned.

A key property of the Discrete Fourier Transform is "Periodicity". According to this, the signal at the input of the Fourier Transform is assumed...periodic. That is, you may be passing a voice recording of "One Two Three" to the FT to recover its harmonic component but the FT assumes that the signal, in reality sounds like "...One Two Three One Two Three One Two Three One Two Three..." extending to infinity both forwards and backwards. For more information about this please see: http://en.wikipedia.org/wiki/Discrete_Fourier_transform (Periodicity property).

The signal that you call "Time Signal", contains a brief period of silence in the beginning and ends abruptly, that is, it would be discontinuous if it were to be a periodic signal. By the way, you may also be interested in the technique of windowing, for more information please see: http://en.wikipedia.org/wiki/Window_function#Symmetry_and_asymmetry

I almost can't extract any information fro the phase plot. What am I missing?

The phase spectrum is simply the angle between the imaginary and real components of the signal spectrum. However, having seen how the imaginary and real components of the FT look like, their phase content is not expected to be a smooth and well behaved function, except perhaps for very simplistic cases.

By the way, you may want to look up the concept of phase unwrapping, to improve your phase plots. For more information, please see https://ccrma.stanford.edu/~jos/fp/Phase_Unwrapping.html

Hope this helps.


Should I Use a High Pass Filter (HPF) or Low Pass Filter (LPF) to Enhance the Features for Image Splicing Detection?

Since you're after a filter which emphasizes abrupt changes you should use High Pass Filter.

The issue is you'd be also sensitive to noise.
Hence one way to do it is to apply High Pass Filter on slightly blurred image (Which actually results in a Band Pass Filter).

One easy choice would be using the Gradient of The Gaussian Filter.
Why is it?
Because it is easy to calculate (Analytically) and efficient to implement because of its separability property of the filter.


I tried to make my derivation match up exactly with embeddedrelated's. The filter you cited uses an approximation, which further obfuscates intuitively understanding it.

embeddedrelated does cite the frequency transfer function that the IIR filter is based off of:
$$
H(s)=\frac{1}{1+s\tau}
$$

The transfer function for a simple RC low pass filter is:
$$
H(s)_{RC}=\frac{1}{1+sRC}
$$
$$
\tau=RC
$$

Schematic:

We can derive a difference equation from the above schematic.

$$
I=C\frac{dV_{out}}{dt}
$$
$$
I = \frac{V_{in} - V_{out}}{R}
$$
$$
C\frac{dV_{out}}{dt}=\frac{V_{in} - V_{out}}{R}
$$
$$
C\frac{V_{out}[n]-V_{out}[n-1]}{ts}=\frac{V_{in}[n] - V_{out}[n]}{R}
$$
$$
V_{in}[n]+\tau f_s V_{out}[n-1]=\tau f_s V_{out}[n] + V_{out}[n]
$$
I have tested this form:
$$
\frac{V_{in}[n]+\tau f_s V_{out}[n-1]}{\tau f_s + 1}=V_{out}[n]
$$
$$
\tau f_s >> 1
$$
I have not tested this form:
$$
\frac{V_{in}[n]}{\tau f_s} + V_{out}[n-1]\approx V_{out}[n]
$$


Is there a closed form expression for main-lobe width increase given a window?

If I understand you correctly, you are in fact asking for a closed-form solution for the position of the first zero of a given window function. Such a solution cannot exist, because what you would need is an analytic solution to a polynomial equation, which does not exist for polynomials of degree 5 or higher (that's the Abel-Ruffini Theorem).


[My path was erroneous, I made a confusion, but kept a modified answer for one-sided $\mathcal{Z}$-transforms]

With $\mathcal{Z}$-transforms, formulae with powers of $n$ often aim at making you handle basic operations on $\mathcal{Z}$-transforms: linearity, time-product/differentiation, integration, time-shift, see for instance: Properties of Z-Transform:

$$ \mathcal{Z}(nx[n]) = -z\frac{d}{dz}X(z)\,,$$
and
$$ \mathcal{Z}(x[n-n_0]) = z^{-n_0}X(z)\,.$$

Hint: if your case were a one-sided $\mathcal{Z}$-transform, one could have focused on derivation, and on the simple geometric series:

$$\sum_{n=0}^{+\infty}z^{-n} = \frac{z}{z-1}\,.$$

As of today, as pointed out by Matt and Robert, I could not make sense of the
$$\sum_{n=-\infty}^{+\infty}z^{-n}\,,$$
and I can see but a zero radius of convergence.
For what it's worth, I will anyway try to add a common trick to try to derive something. One can use a stabilization with  some $0 < a <1$:

$$\sum_{n=-\infty}^{+\infty}a^{|n|}z^{-n}\,,$$
with the secret thought that, if $a$ is just close to one, we could get almost what we wanted.
I may have made mistakes again, so far this yields the formal form:

$$ \frac{(1-a^2)z}{(z-a)(1-az)}\,.$$

This still have a singular ROC, but one can play with  $\mathcal{Z}$-transform properties.