The advantage of the complex DFT or complex Fourier transform or
complex Fourier series is that linear
systems have the nice property that the response to $A\exp(j\omega t)$ is
$H(\omega)A\exp(j\omega t)$.  (Here $A$ can be a complex constant).
So the output is just a scalar multiple of the input.
More importantly, if we have a representation of the input
as a weighted sum of complex exponentials, the output is
just another weighted sum of the same exponentials.  Different
weights, but same set of exponentials. Furthermore, each new weight
is obtained by multiplying the old weight by an appropriate number.

Of course, no physical system has complex-values
signals going in and coming out; at least, not as of today though
one can always hope for better things in the future.  In the mean time,
we take real parts of the complex signals, or get the response to $\cos(\omega t)$ or
$\sin(\omega t)$ via linearity and superposition and liberal use of
$$\begin{align*}
\cos(\omega t) &= \frac{\exp(j\omega t) + \exp(-j\omega t)}{2}\\
\sin(\omega t) &= \frac{\exp(j\omega t) - \exp(-j\omega t)}{2j}
\end{align*}$$

In contrast, the response to $\cos(\omega t)$ is of the form
$B(\omega)\cos(\omega t) + C(\omega)\sin(\omega t)$.  So,
while linearity and superposition etc all work, the output
might well need the use of different basis functions than the
input does.  Very closely related, of course, but still possibly
different and maybe more  basis functions might be needed. For
example, input $\cos(\omega t)$ is represented by one basis function,
output $B(\omega)\cos(\omega t) + C(\omega)\sin(\omega t)$ by
two basis functions.  It can be argued that complex functions
require twice as much work as real functions and so any savings
are purely imaginary (pun intended), but complex representations allow
uniform treatment while sin/cos representations do not.  Quick!
Given the response to $\cos(\omega t)$ is
$B(\omega)\cos(\omega t) + C(\omega)\sin(\omega t)$,
what is the response to $\sin(\omega t)$?  You have to work at it
a bit, you may need to invoke formulas such as
$$\cos (\alpha + \beta) = \cos(\alpha)\cos(\beta) - \sin(\alpha)\sin(\beta)$$
and so on.  With complex exponentials, life is a lot easier.

But, as in real life, your mileage may vary, and if you feel
that sin/cos representations are the way to go and complex
exponentials should be eschewed, you are free to follow your heart.
If you have difficulty communicating your ideas to colleagues,
bosses, clients or consultants, that will be their loss, not yours.


FFTs of a complex signal - separating the real and imaginary parts

Below is a chart I had of "Universal Fourier Transform Properties", that apply in either direction (going from time to frequency or going from frequency to time).  For example, a signal that is periodic in one domain, will be discrete in the other: a digitized time domain signal becomes periodic in frequency (such that we can concern ourselves with the spectrum from $-F_s/2$ to $+F_s/2$ only. Similarly a signal that is periodic in time will only have discrete frequencies (at the fundamental repetition rate and its harmonics in frequency).

Specific to your case I draw attention to the properties when a signal is only real, and when a signal in only imaginary:

This may help you to see what is going on, specifically the FT of a signal that is only real in the time domain, and the FT of a signal that is only imaginary in the time domain. A signal that is only real will have a spectrum that is conjugate symmetric: the positive half of the spectrum is equal in magnitude to the negative half but opposite in phase (Meaning the real portion of the FFT is "even" and the imaginary portion of the FFT is "odd", in the same fashion that a cosine in time is an even function and a sine in time is an odd function. Similarly a signal that is only imaginary will be have the FFT real components to be odd and the FFT odd components to be even.

So hopefully you see the difference between taking the fft of b, and taking the fft of jb, where to be clear b is real, so jb is imaginary (we could complicate this by allowing a and b to be complex numbers but my guess in your case they are indeed real).

This is also interesting in showing how causal time domain signals MUST be complex in frequency, and the imaginary and real components are related by the Hilbert transform: if you consider the real and imaginary components of the spectrum seperately, we see via the odd/even relationship that everything for t<0 cancels while everything for t>0 adds.

This last point is further elaborated by Sir Robert Bristow-Johnson in this post: What is the easiest, most straight-forward way to prove this about minimum-phase filters?


I think they made a mistake. $2^{1/3} \approx 1.259921$ is not close enough to $10^{1/6} \approx 1.467799$ to say that there is any coincidence. Even the half-octave factor $2^{1/2} \approx 1.414214$ is closer.


The paper that you are quoting assumes that the motions along the different coordinates are independant functions, i.e., $$P(t) = \left(x(t),y(t)\right) = (f_1(t),f_2(t)) \neq f(x,y)(t),$$ where $P$ stand for position.
The various 1D functions are obtained in this way.

A bit of smoothing is generally helpful.
However, when to comes to gait recognition, you will probably have to test different smoothing procedures: my intuition is that gait recognition relies on high frequency / low amplitude clues that can easily be cancelled out (but I may be wrong for this).


defining the frequency range with fast fourier transform

There are a few things going on here.  First, one of the fundamental properties of the Fourier transform is that it conserves energy, i.e. the energy of the time domain is exactly equal to the energy of its Fourier transform.  What is the power of a sinusoid?  Assuming we restrict the calculation to integer multiples of cycles it is as follows (I suggest working out this equation for yourself if you are not sure about it):

$$
P(Asin(ft)) = \frac{A^2}{2}
$$

The energy of a sinusoid is thus $\frac{A^2}{2}T$, or, in the case of discrete data, $\frac{A^2}{2}N$.  The purpose of the following line of code-

X = fft(x)./N;

is to change the FFT output from an energy-based measure to a power-based measure ("based" because you have to square the transform results, just like you would in the time domain, to get the actual energy or power).  Power is usually more interesting since it is an inherent characteristic of the signal, whereas energy is dependent on the number of samples that you are using.

So, the power of your sinusoids is, according to the equation, 4.5, 2, and 0.5.  You got 1.5, 1, and 0.5, which, if you square them, become 2.25, 1, and 0.25.  You will notice that those results are exactly one half your expected power.  The reason for the halving is that you are only showing the positive frequencies from your Fourier transform.  That is what the following lines of code are doing-

cutOff = ceil(N/2);
X = X(1:cutOff);

They are cutting off the negative frequency part of the spectrum.  If you plotted the entire set of results it would look like the following-

The spikes on the left are the positive frequencies that you saw before.  The spikes on the right are the negative frequency components of the sinusoids.  Note that they have the exact same magnitude as their positive frequency partners, and thus comprise the missing factor of two in the power.

Lastly, you ask the following-

  Also, could someone explain the section

% define the frequency components of the series
freq = 0:N-1;
freq = freq*fs/N;

Samples sampled at frequency $f_s$ can represent frequencies from $-\frac{f_s}{2}$ to $\frac{f_s}{2}$ (see Nyquist Sampling Theorem), or, equivalently, frequencies from $0$ to $f_s$ (see aliasing).  The code implicitly assumes that the samples represent from $0$ to $f_s$.


Inverse DTFT of phase shifted complex exponential

HINTS:

$e^{j\pi/2}=j$, so $X(e^{j\omega})=j\omega e^{-j10\omega}$
The term $e^{-j10\omega}$ corresponds to a delay of $10$ samples in the time domain, so you can just ignore it for now, and then add a delay of $10$ samples. So you compute the IDTFT of $j\omega$, and then replace the variable $n$ of the result by $n-10$.
The IDTFT of $j\omega$ is $$\frac{1}{2\pi}\int_{-\pi}^{\pi}j\omega e^{jn\omega}d\omega$$ That integral is most easily computed using integration by parts.

Note that $X(e^{j\omega})$ is the frequency response of an ideal discrete-time differentiator with a delay of $10$ samples.


Amplitude of a signal. [-1, 1] vs variants and MFCC implications

What do these values actually mean? Are they arbitrary?

Yes, pretty much: they will just be the scale of numbers your system works with. In floating point systems, we often see that samples get normalized to [-1,+1], whereas in fixed-point system, it's often things like [-2⁻¹⁵,+2¹⁵-1], depending on the bit width of the samples to begin with.

So, this is purely conventional. Most algorithms don't care about the absolute scale of things.

So, if your algorithm does depend on amplitudes, well, normalize your signal before you process it.


A filter in IOS to eliminate a certain frequency range

You can think of a discrete signal as an array of values, a filter as another array of values (not necessarily of the same size), and filtering as performing an operation called convolution on both the filter and input arrays.

You can get the filter array by using programs such as Matlab, scientific python or octave (the last two are free). You can implement the filtering operation by coding the convolution operation in your own program.

Be aware that convolution can be cpu-intensive. The 'better' the filter (or, the higher its order), the larger the filter array. A better filter rejects more of the signal outside the 16-18 kHz range, and distorts less the signal inside the range. However, having more elements, a convolution using it requires more operations (basically multiplications and additions). This may be a concern on a mobile device.

See this page for detailed descriptions (including source code) of scientific python commands that do what you need. Look in particular for convolve and firwin2.


Iterative Blind Sinus Signal Suppression

I assume you do everything in discrete time. If I understood correctly then the observed signal is

$$x(n)=s(n)-A\sin(\omega_0n+\phi)$$

where $s(n)$ is a sinusoidal signal with frequency $\omega_0$ and unknown amplitude and phase. The parameters $A$ and $\phi$ are to be optimized such that the power of $x(n)$ is minimized. In a real-time algorithm you can define the cost function as the sum of $x^2(n)$ over a fixed length $N$:

$$F(A,\phi,n)=\sum_{k=n-N+1}^{n}x^2(k)$$

where $n$ is the current sample. For a steepest descent algorithm you need the gradient:

$$\frac{\partial F}{\partial A}=-2\sum_{k=n-N+1}^{n}x(k)\sin(\omega_0k+\phi)\\
\frac{\partial F}{\partial \phi}=-2A\sum_{k=n-N+1}^{n}x(k)\cos(\omega_0k+\phi)$$

The parameters are now updated in the following way:

$$A_n=A_{n-1}-\mu\frac{\partial F_{n-1}}{\partial A}\\
\phi_n=\phi_{n-1}-\mu\frac{\partial F_{n-1}}{\partial \phi}$$

where $\mu$ is the step size. Note that in every iteration you need to recompute the gradients with the current values of $A$ and $\phi$.

The figure below is a simulation of this gradient algorithm with the unknown signal $s(n)=2\sin(\pi n/3+\pi/2)$, $N=50$, $\mu=0.01$, and with zero initialization of both amplitude and phase (amplitude $A$ in blue, phase $\phi$ in green).


If humans can only hear up to 20 kHz frequency sound, why is music audio sampled at 44.1 kHz?

The sampling rate of a real signal needs to be greater than twice the signal bandwidth. Audio practically starts at 0 Hz, so the highest frequency present in audio recorded at 44.1 kHz is 22.05 kHz (22.05 kHz bandwidth).
Perfect brickwall filters are mathematically impossible, so we can't just perfectly cut off frequencies above 20 kHz.  The extra 2 kHz is for the roll-off of the filters; it's "wiggle room" in which the audio can alias due to imperfect filters, but we can't hear it.
The specific value of 44.1 kHz was compatible with both PAL and NTSC video frame rates used at the time.

Note that the rationale is published in many places: Wikipedia: Why 44.1 kHz?


Suitable sampling rate for triangle wave

The 5kHz is the repetition frequency and it's not the maximum frequency of your signal or its bandwidth and your signal isn't band-limited but it has small frequency content over very high frequencies which lead to the Fat32 answer.

but there are other sampling and reconstruction scheme for non band-limited signals like Finite Rate of Innovation (FRI) which could perfectly reconstruct certain families of non band-limited signals like piecewise polynomials (which your signal exactly fits in) with finite sampling rate (in your case 6 sample per each repetition).


Omega-K algorithm implementation for synthetic aperture radar

I am using this link for getting the reference by Richard Balmer :

There are no references 30 and 32


Intuitive explanation of cross-correlation in frequency domain

The concept is based on the convolution theorem, which states that for two signals $x(t)$ and $y(t)$, the product of their Fourier transforms $X(f)$ and $Y(f)$ is equal to the Fourier transform of the convolution of the two signals. That is:

$$
\mathcal{F}\{x(t) * y(t)\} = \mathcal{F}\{x(t)\}\mathcal{F}\{y(t)\}
$$

You can read more on the derivation of this theorem at the above Wikipedia link. Now, convolution is a very important operation for linear systems in itself, so the theory on its properties is well-developed.

However, what you're looking for is the cross-correlation between $x(t)$ and $y(t)$. Here's the key: you the cross-correlation integral is equivalent to the convolution integral if one of the input signals is conjugated and time-reversed. This allows you to utilize theory developed for evaluating convolutions (like frequency-domain techniques for calculating them quickly) and apply them to correlations.

In your example, you're calculating the following:

$$
\mathcal{F}\{x(t)\}\left(\mathcal{F}\{y(t)\}\right)^*
$$

Recall that in the Fourier domain, complex conjugation is equivalent to time reversal in the time domain (this follows directly from the definition of the Fourier transform). Therefore, using the first equation given above, we can state that:

$$
\mathcal{F}\{x(t) * y^*(-t)\} = \mathcal{F}\{x(t)\}\left(\mathcal{F}\{y(t)\}\right)^*
$$

If you then take the inverse Fourier transform of this equation, the signal you're left with is the cross-correlation between $x(t)$ and $y(t)$.


How do I find the peak frequency of a wav file in python?

import sounddevice as sd
import scipy
import numpy as np
from scipy.io import wavfile as wav
from scipy.io.wavfile import write ,read
from scipy import fftpack as scfft
from matplotlib import pyplot as plt

fs = 44100
seconds = 3
myrecording = sd.rec(int(seconds * fs), samplerate=fs, channels=1)
sd.wait()
write('output.wav', fs, myrecording)

fs_rate, signal = read("output.wav")
l_audio = len(signal.shape)
N = signal.shape[0]
secs = N / float(fs_rate)
Ts = 1.0/fs_rate
t = scipy.arange(0, secs, Ts)
FFT = abs(scipy.fft(signal))
FFT_side = FFT[range(N//2)]
freqs = scipy.fftpack.fftfreq(signal.size, t[1]-t[0])
fft_freqs = np.array(freqs)
freqs_side = freqs[range(N//2)]
fft_freqs_side = np.array(freqs_side)

volume=np.array(abs(FFT_side))
audible=np.where(volume>5)

HighestAudibleFrequency=max(freqs_side[audible])
print(HighestAudibleFrequency)

The above code records audio for 3 seconds from the microphone and outputs the highest frequency detected in the 3 seconds. Sorry for the code format, i dont know how to make it look better


Signals and Systems - LTI - Transforms - Impulse Response

$$H(jw) = \frac{output}{input} = \frac{Y(jw)}{X(jw)}$$

$$\frac{Y(s)}{X(s)} = \frac{ab}{(a+s)(b+s)}$$

$$s=jw$$

Then you can use some of the notions from this example to solve for your impulse response / differential equation. Hope that helps!


Hydrophones have a weird way of specifying sensitivity. $1V/ \mu Pa$ simply means $10^6 V/Pa$ which is frankly an idiotic unit.

You scale this by -200dB and  you get $ S = 10^{-4} V/Pa$ or $S = 0.1mV/Pa$ which is a much more reasonable representation.

This has NOTHING to do with dBSPL which is referenced to the human threshold of hearing. There is no reference needed here, since you are using Pascals directly.


How to uniquely identify phase difference between harmonics?

i am assuming this is a periodic (or quasi-periodic) function.  also seems like the signal is a musical tone.  periodic signals have a period, $P$, and the reciprocal of that period is the fundamental frequency; $f_0 \triangleq \frac{1}{P}$

or, as angular frequency; $\omega_0 = 2 \pi f_0 = \frac{2 \pi}{P}$ .

the periodic signal is $x(t)$ where

$$ x(t+P) = x(t)  \qquad \forall t \in \mathbb{R} $$

it can be represented as a Fourier series

$$ \begin{align}
x(t) \ &= \ a_0 \ + \ 2 \sum_{k=1}^{\infty} a_k \cos(k \omega_0 t) - b_k \sin(k \omega_0 t) \\
 &= \ a_0 \ + \ 2 \sum_{k=1}^{\infty} a_k \frac{e^{j k \omega_0 t} + e^{-j k \omega_0 t}}{2} \ - \ b_k \frac{e^{j k \omega_0 t} - e^{-j k \omega_0 t}}{2j} \\
 &= \ a_0 \ + \ \sum_{k=1}^{\infty} (a_k + jb_k) e^{j k \omega_0 t} \ + \ (a_k - jb_k) e^{-j k \omega_0 t} \\
 &= \ \sum_{k=-\infty}^{\infty} c_k \ e^{j k \omega_0 t} \\
\end{align} $$

where

$$ c_k =
\begin{cases}
a_{-k} - jb_{-k}, & \quad \text{for } k < 0 \\
a_0, & \quad\quad k = 0 \\
a_k + jb_k, & \quad\quad k > 0
\end{cases} $$

and, going in the other direction,

$$ \begin{array}{lcl}
a_0 & = & c_0 \\
a_k & = & \Re\{c_k \} \quad\quad \text{for } k>0\\
b_k & = & \Im\{c_k \} \quad\quad\quad k>0
\end{array} $$

note that for real $a_k$ and real $b_k$, then

$$ c_{-k} = c_k^* = \text{"complex conjugate of } c_k \text{ "}$$

it's easier to get $c_k$ from $x(t)$ and then compute $a_k$ and $b_k$ from the $c_k$ which is:

$$ c_k = \tfrac{1}{P} \int\limits_{t_0}^{t_0+P} x(t) \, e^{-j k \omega_0 t} \, \mathrm{d}t \qquad \text{for any } t_0 \in \mathbb{R}$$

you can uniquely define the phase of overtones (or harmonics apart from the fundamental) with reference to the fundamental if there is any energy at the fundamental.  (it is not always the case that a periodic signal has non-zero energy in its own fundamental frequency.)

now it's also true that:

$$ \begin{align}
x(t) \ &= \ \sum_{k=-\infty}^{\infty} c_k \ e^{j k \omega_0 t} \\
 &= \ a_0 \ + \ 2 \sum_{k=1}^{\infty} a_k \cos(k \omega_0 t) - b_k \sin(k \omega_0 t) \\
 &= \ a_0 \ + \sum_{k=1}^{\infty} r_k \cos(k \omega_0 t + \phi_k) \\
\end{align} $$

where

$$ \begin{align}
r_k &= 2 \big| a_k + j\,b_k \big|  \\
 &= 2 \sqrt{a_k^2 + b_k^2}  \\
\\
\phi_k &= \arg \{a_k + j\,b_k \}  \\
 &= \operatorname{atan2}(b_k, a_k)  \\
\end{align} $$

now, the phases of all of the harmonics are relative to some time $t=0$, but that reference time can be defined in such a way that the phase of the fundamental is zero.

let

$$\begin{align}
 \tilde{x}(t) \ &\triangleq \ x(t - \tfrac{\phi_1}{\omega_0}) \\
 &= \ a_0 \ +  \sum_{k=1}^{\infty} r_k \cos \big(k \omega_0 (t - \tfrac{\phi_1}{\omega_0}) + \phi_k \big) \\
 &= \ a_0 \ +  \sum_{k=1}^{\infty} r_k \cos(k \omega_0 t + \phi_k - k\phi_1) \\
 &= \ a_0 \ +  \sum_{k=1}^{\infty} r_k \cos(k \omega_0 t + \tilde{\phi}_k) \\
\end{align} $$

where $ \tilde{\phi}_k \triangleq \phi_k - k\phi_1 $.

you can see that $\tilde{\phi}_1 = 0$, that we defined our origin of $t$ so that the phase of the fundamental (when $k=1$) is 0.  then all of the other harmonics have their phases referenced to that "first harmonic" or fundamental.  the only problem is that the amplitude of that fundamental must not be zero.  $r_1 \ne 0$ so that the reference phase, $\phi_1 = \arg \{a_1 + j\,b_1 \}$, can be determined.


I'm thinking the signals could be compared by

$D = \mathcal{F}^{-1}(\mathcal{F}(A) - \mathcal{F}(B))$

Where A and B are the two audio files to be compared.

I guess after that you could come up with some criteria for deciding that the voices are different. Perhaps when $|D|$ goes over some threshold.

You could also use Cross correlation

Using some other library to help with this would save time. But I haven't done this with audio files.

AForge Looks like it has FFT although I haven't used it.


STFT multiplied with Gaussian window vs STFT multiplied by rectangular window

Using a Gaussian window has one advantage over all other windows-- it yields the most compact time-frequency view of the spectrum. In other words, the product of the uncertainty in the time and frequency domains is minimized. This is due to the uncertainty principle of the Fourier transform.

Incidentally, in physics this property of the Fourier transform is the mathematical underpinning behind the Heisenberg uncertainty principle, which states that the position and momentum cannot simultaneously be known to infinite precision. Replace "position" with "time" and "momentum" with "frequency" in that statement and you can see how this applies to spectral analysis.

(There are other reasons not to use a Gaussian window as well, depending on what your application requires.)

Also, incidentally, Gabor was a physicist who first applied many mathematical tools from physics to signal processing problems, such as the analytic signal representation. It would not surprise me if his use of the Gaussian window in signal processing was carried over from quantum mechanics.


Is there a separate word for the time domain when working with digital signals?

I use sampled-time domain or discrete-time domain.  I think I've heard other people use other terms, so they can chime in (discretely).


Is there any literature discussing PDF after quantization?

any discrete random variable has a p.d.f. that is a summation of dirac delta functions.

$$ p_\mathrm{y}(\alpha) = \sum\limits_i P_i \ \delta(\alpha - y_i) $$

where $\sum\limits_i P_i = 1 $.

if $y[n]$ is the quantization of $x[n]$:

$$ y[n] = \Delta \bigg\lfloor \frac{x[n]}{\Delta} \bigg\rfloor $$

where $ \lfloor \cdot \rfloor$ is the floor() function and $\Delta$ is the quantization step size.

then $y_i = i \Delta $ for some set of integers $i$.  the p.d.f. weighting constants become

$$ P_i = \int\limits_{i \Delta}^{(i+1)\Delta} p_\mathrm{x}(\alpha) \, d \alpha $$


Very often not.  Your posted code is for finding an FFT magnitude, not an autocorrelation result, and thus likely to fail at estimating the fundamental pitch of sound/voice/music with typically high overtone or harmonic content.


Can two modulations be used in the same frame.

Yes. Nobody's stopping you....

  For qam/qpsk/psk/ modulation scheme for the data bits can i have a fixed preamble with psk modulation?

If you want, why not? It's actually pretty common that systems with a flexible modulation scheme (e.g. IEEE 802.11, "wlan") have a header symbol that's in a fixed modulation that tells the receiver what the rest of the frame is.


Understanding Communication Theory Jargon for Statistical Linear Models

The equation $y =  hs + n$ is a discrete-time model of a certain kind of wireless communications system.

$s$ is a complex number, drawn from a finite set called a "constellation", that represents the information being transmitted.
$n$ is a complex, circularly-symmetric Gaussian random variable with zero mean and variance usually denoted by $N_0$.
$h$ represents the effect of the communications channel. In this simple model of a wireless channel, $h$ is a complex Gaussian random variable with zero mean and variance one. This model is called Rayleigh-fading because $|h|^2$ is a Rayleigh-distributed random variable.
$y$ is the received complex number. From it, the receiver must make an estimate of $s$. This estimate is usually calculated by first estimating $h$, and then finding $$ y' = \frac{h^* y}{|h|^2} = s + \frac{h^*n}{|h|^2} $$ Here you can see the importance of $|h|^2$ to enhance the noise, which is why the channel is called Rayleigh. The estimate of $s$ is calculated by finding the constellation element that is closest (in Euclidean-distance sense) to $y'$.


Removing Bounded Noise from an Audio Signal

If you have some prior information about your noise e.g frequency band, you can simply implement an noise removing algorithm with emphasis on removing the noise on that band (I mean use a somehow weighted noise removing process).

One of nice methods is weighted multi-band spectral subtraction (this is not the only method and in fact these weighted methods have many variations).
The improvement of MBSS is due to the fact that the multi-band speech takes into account the non-linear effect of coloured noise on the spectrum of speech i.e. some frequencies are affected more adversely than others.
Multi-band spectral subtraction (MBSS) method provides a definite improvement over other methods such as the conventional power spectral subtraction method.

Speech spectrum is divided into $N$ non-overlapping bands and spectral subtraction is applied to each band independently. The estimation of the clean speech spectrum in $i$'th band can be obtained from:

$$|x(\omega)|^2=|y(\omega)|^2 – \alpha_i \delta_i |d(\omega)|^2$$

for

$c_i<\omega<c_{i+1}$

Where $ c_i $ and $ c_{i+1}$ are the beginning and ending frequency band (your specific band) respectively. In the $i$'th band, $\alpha_i$ is the over-subtraction factor and $\delta_i$  is factor to control noise removal properties in each band individually. So you can set coefficient to remove the bounded noise from an audio signal.

For more help, you can read this paper :

  Sunil D. Kamath, Philipos C. Loizou, “A multi-band spectral
  subtraction method for enhancing speech corrupted by colored noise,”
   Proc. IEEE Int. Conf. Acoust., Speech, Signal Process., vol. 4, July 2002.


Is it valid to calculate the magnitude, power, and phase of a real time-domain signal without converting it to the frequency-domain?

No, this doesn't make much sense. What you can do in the time domain is compute the analytic signal and derive the signal's instantaneous amplitude (envelope) and its instantaneous phase from it.

Take as an example

$$x[n]=A\sin(\omega_0n+\theta)\tag{1}$$

The corresponding analytic signal is

$$x_a[n]=-jAe^{j(\omega_0n+\theta)}\tag{1}$$

Its instantaneous amplitude (envelope) is $\big|x_a[n]\big|=A$, and its instantaneous phase is $\arg\{x_a[n]\}=\omega_0n+\theta-\pi/2$.


x = [0 0 0 0 1 1 1 1] is one cycle of a (not bandlimited) square wave, with a DC offset, with a period equal to the period of the window.

....''''....''''....''''

The DC component is the FFT value at index 1 (remember that Matlab's array indexing starts at 1 instead of the more sane 0-based indexing), and to get the true value you have to divide by the number of samples (4/8 = 1/2 DC offset).  The square wave is the components at index 2 and 4.  Because the signal is real, the components at 6, 7, 8 are just mirror images of 4, 3, 2.

A square wave is made up of odd harmonics with 1/n amplitude, so with fundamental f1

frequency   f1, amplitude   1
frequency 3*f1, amplitude 1/3
frequency 5*f1, amplitude 1/5
frequency 7*f1, amplitude 1/7
...

Since your number of samples is so low, you only have room for the first 2 harmonics before you hit the Nyquist frequency.  So the value at index 4 is the 3rd harmonic, which should be approximately 1/3 of the value at index 2, the first harmonic (though not exactly 1/3 because of aliasing).

If your sampling frequency is 1 MHz, and the value at index 2 represents the fundamental frequency, then you know the square wave fundamental frequency is (index-1)*fs/N = (2-1)*1MHz/8 = 125 kHz, or 1/8th of the sampling frequency, as you would expect.

Try experimenting with longer signals and it will be easier to see what's happening.


Why should I zero-pad a signal before taking the Fourier transform?

Zero padding allows one to use a longer FFT, which will produce a longer FFT result vector.

A longer FFT result has more frequency bins that are more closely spaced in frequency.  But they will be essentially providing the same result as a high quality Sinc interpolation of a shorter non-zero-padded FFT of the original data.

This might result in a smoother looking spectrum when plotted without further interpolation.

Although this interpolation won't help with resolving or the resolution of and/or between adjacent or nearby frequencies, it might make it easier to visually resolve the peak of a single isolated frequency that does not have any significant adjacent signals or noise in the spectrum.  Statistically, the higher density of FFT result bins will probably make it more likely that the peak magnitude bin is closer to the frequency of a random isolated input frequency sinusoid, and without further interpolation (parabolic, et.al.).

But, essentially, zero padding before a DFT/FFT is a computationally efficient method of interpolating a large number of points.

Zero-padding for cross-correlation, auto-correlation, or convolution filtering  is used to not mix convolution results (due to circular convolution).  The full result of a linear convolution is longer than either of the two input vectors.  If you don't provide a place to put the end of this longer convolution result, FFT fast convolution will just mix it in with and cruft up your desired result.  Zero-padding provides a bunch zeros into which to mix the longer result.  And it's far far easier to un-mix something that has only been mixed/summed with a vector of zeros.


Computing overall bitrate of a transmitter with a two chunks package

Don't use fancy calculations which merely muddle up the problem. The header is $K$ bits where $K = Lk_1$; the body is $M$ bits where $M = Nk_2$. You are using a $[n_1, k_1]$ code to transmit the header which you do by dividing the $K$ bits of the header into $L$ chunks of $k_1$ bits each, encoding them, and thus sending $L$ codewords of length $n_1$ bits each. Total number of channel bits $= Ln_1$. Similarly, you send a total of $Nn_2$ channel bits ($N$ codewords of length $n_2$ each) giving a grand total of $Ln_1+Nn_2$ channel bits to send the $Lk_1+Nk_2$ "data" bits in the header and body.  So, the overall rate of the scheme is
$$R = \frac{Lk_1+Nk_2}{Ln_1+Nn_2}.$$
But be aware that some people refuse to acknowledge the information in the header as data and insist that actually only $Nk_2$ data bits are being transmitted; everything else is overhead. For these folks, the rate would be just
$$R = \frac{Nk_2}{Ln_1+Nn_2}.$$


How can I find the transfer function of the following block diagram?

In such cases it can be useful to introduce auxiliary variables at the input of the integrators. For the given diagram you could use a signal $u(t)$ at the input of the second integrator. The equation for its Laplace transform $U(s)$ becomes
$$U(s)=\frac{1}{s}\big[X(s)-bY(s)\big]-aY(s)\tag{1}$$
You need another equation relating $U(s)$ to $Y(s)$, but that one is trivial. From those two equations you can express $Y(s)$ in terms of $X(s)$, and in this way you can obtain the transfer function $H(s)=Y(s)/X(s)$. The solution in your question is not correct. The correct solution must be a rational function, not a polynomial.


Strictly speaking you can't because without specifying the ROC, the inverse Laplace transform is generally not unique. However, in many contexts there is the implicit assumption of causality of the corresponding time function (i.e., $x(t)=0$ for $t<0$), which is equivalent to stating that the ROC is a right half-plane.


Calculating similarity between bit strings/signal

You should indeed look into sequence alignment methods used in bioinformatics. Those can penalize both for insertions/deletions and mismatches. The binary alphabet having a size of 2 and that of DNA having a size of 4 does not require the use of different algorithms if the element-wise comparison results in either "match" or "mismatch", and a mismatch is penalized for. Check out the Needleman–Wunsch algorithm, which is based on dynamic programming.

Some DNA/protein sequence alignment algorithms such as MAFFT use Fast Fourier Transform (FFT) for accelerated cross-correlation of residue properties to approximate match calculations.


Can some one explain what is the meaning ( in real life) of anti causal (precursor) part of the channel taps?

Actual channels are always causal (like everything else in the physical universe).

Actual (discrete-time) channels also sometimes have one tap that is considerably larger than the rest; an example impulse response would be h = [0.1, 1.5, 0.2]. Some authors prefer to define h[0] as the largest tap; in my example, we'd have h[-1] = 0.1, h[0] = 1.5, and h[1] = 0.2. Then, h[-1] would be "anticausal" since it occurs before the impulse is applied at n = 0.

As Marcus says, this is just a model. It has zero effects on the receiver, except on how you write the indices in your algorithms.

Personally, I never worry about this and always define the first tap as h[0] (or rather h[1], if implementing it in Matlab or Julia).


That's a simple process, where $y(k)=u(k-3)$, as the $z^{-1}$ blocks just perform a delay by one sample. $x_1(k)$, $x_2(k)$ and $x_3(k)$ are just arbitrary names. In other words: output will be exactly like input, only delayed by three samples.


Need help finding the frequency response

Is this for homework?

If so I won't directly give an answer, just a hint about the general approach that can be used to solve such questions without a single derivation:

Plot the impulse response. How does it look? Try building this by combining or modifying (through addition, dilatation, translation, multiplication, convolution, derivation, integration...) very simple functions the DFT of which is known (complex exponentials, constants, rectangles, diracs, Heaviside step). From then, use identities on the DFT to get the DFT of your impulse response.


Match Filter in Time-Frequency domain instead of just Time domain. Redundant, or better?

Think of the time - frequency ambiguity of your matched filter like so:

Frequency ambiguity means it will respond to a range of frequencies
Time ambiguity means the response will be 'smeared' around the spatial location.

If you you have 0% frequency ambiguity, the matched filter must look like a sine wave and go on for ever, which in the frequency spectrum looks like a dirac delta.

0% time ambiguity is a dirac delta in the time domain.

So if you have a matched filter that is more than 1 sample wide in the time domain, then it already is ambiguous in both time and frequency domains.

If you are doing matched filtering of the envelope, then you are just looking at the modulating signal, and there is no need to look at the 2D time-frequency spectrogram.

If you want to match the envelop (modulating signal) and the base frequency, then you need a quadrature filter with a bandwidth around the range of frequencies you expect. A quadrature filter is required because it will make the response invariant to the phase of the base signal.

If you don't know the base frequency, a 2D time-frequency spectrogram will be useful as it will show you what frequency is being modulated. Essentially, the spectrogram is response of the signal - time axis - to a bunch of different (centre) frequency quadrature filters - frequency axis.

TLDR:

The premise that the evelope-matched filter in the time domain is 100% localised is incorrect.


Take a (real-valued) low pass filter with impulse response $h_{LP}[n]$ and with cut-off frequency $(f_2-f_1)/2$, where $f_1$ and $f_2$ are the lower and upper band edges of the desired band pass filter. Define a center frequency $f_0=(f_1+f_2)/2$. Then a real-valued band pass filter is obtained by

$$h_{BP}[n]=2h_{LP}[n]\cos(2\pi f_0n)\tag{1}$$

Depending on the design of the low pass filter, the filter in $(1)$ can be identical to the one in your question.

If you want a complex-valued band pass filter with a pass band at positive frequencies (but not at negative frequencies), you simply have to modulate the low pass impulse response $h_{LP}[n]$ by a complex exponential instead of by a cosine:

$$h^c_{BP}=h_{LP}[n]e^{j2\pi f_0n}\tag{2}$$


Implementing wavelet transform for finding transients in the power supply

I would say that a matching Mother wavelet could be the best for detecting a transient; but both the selection and implementation would be much slower. The old adage: do you want quantity or quality :) or Bandwidth vs. noise.  Life doesn't come in our neat intellectual packets.
BTW: the simplest technique is an isolation capacitor and a threshold detector (sad to say I know a lot of ways to "cheat"  ).

Evaluation process:

I would look to octave or scilab for a simulator to use.
Construct a model of transients you think are likely or that you are looking for.
Pass the transients through various canned analysis routines and then
Find the transients in the output; the crisper the "find" the more applicable the analysis is.

Reconstruct the input and see how good the reconstruction is; i.e. apply an error criterion to achieve a faithfulness measure.

If this seems tedious it's because it can be: but at each point in the process you get an idea about the effort involved.

Selection time

Design time and complexity

Execution time

Accuracy.

And most importantly:

Knowledge!

Like the wavelet theory: start coarsely/crudely and refine the evaluation process.

I looked on google and found:
----http://www.soest.hawaii.edu/MET/Faculty/bwang/bw/paper/wang45.pdf
A real application of Waveform/Wavelet analysis of El nino and such; with all the hair:
--https://inst.eecs.berkeley.edu/~ee225b/sp14/lectures/shorterm.pdf
Which has a description of the process and history.

--http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.506.6298&rep=rep1&type=pdf
Which discusses real usage ambiguity.

--Why Wavelet developed when we already had Short-time Fourier transform
Which links to another discussion.

--https://en.wikipedia.org/wiki/Wavelet

Back to question:
Most of the techniques are "complete" in the sense that the transforms are invertable and can recover the original input.  In that sense all invertable transforms are equivalent.  But like I said they isolate different features.  Wavelet transforming is a mathematical process for utilizing Mother wavelets to isolate information in time and space. STFT is very specific whereas there are innumerable of "wavelet transforms" each having different sampling Mother wavelet shapes that are applied recursively to give many (time, freq) graphs; two dimensional pictures.
STFT and Gabor and such are (as far as I know) specific but closely correspond to normal fourier/z-transform analysis.

There are a lot of c wavelet programs in the open source domain.

An effort to clarify "disturbance" analysis.

First you have to  be able to define what a disturbance is and then the allowable type one/type two results; i.e. false positives/negatives and then do the system design around this.

Let me list a few "disturbances"; not necessarily related to power supplies.

For a nuclear containment vessel we tried to "listen" for what amounts to tearing at 100khz-1Mhz with transducers.  A time exponential increase in the power of the signal was taken as (and experimentally verified) an indication of a problem.

As far as power supplies go:
There are safety related issues such as over-voltage causing downstream failures.  For these I always put in absolute hardware; zeners, surge suppressors.  That I think will limit damage and fail-safe; shut things down.  Do the thermal calculations with respect to the calories produced by hardware shunts.

Under-voltage is also a problem in some cases, as I found out by almost smashing a co-workers hand.  Dropping the field voltage to a DC shunt drive motor causes it to speed up :)  So "pulling the plug" at random is not a great idea; and a drop in field voltage should probably be acted on.

Other expected cases are parasitic resonances/oscillations occurring in the output power stage of supplies.  These are typically produced by altered frequency response characteristics of bipolar output transistors during transitions.  Typically these are real transients and not worrisome unless some downstream equipment is sensitive.

You have specification errors where for some reason, say a passing cyclotron, the system is substantially normal but some specified tolerance is breached.  This can also happen slowly when equipment ages and drifts out of spec.

Of course in motor-generator situations transient signals might indicate mechanical wear.

Then you have fast raw supply transients that find their way through the regulators and occur at the output.  That's the reason I use hardware suppresors/zeners that can absorb some over-voltage attempts without allow ing the disturbance through and recover easily.


What is the significance of Butterworth filter in image processing?

Some important differences:

In Image Processing (IP), there is no causality like in Signal Processing (SP), hence there is not a tradeoff between filter quality and sampling sequence.
In IP, the FIR versions of SP are preferred instead of the IIR version (which are rare as you pointed). A possible relevant cause for this is FIR are designed as linear phase, unlike IIR which cannot be linear phase. Linear phase in FIR means the delay is constant (some pixels or fractions of pixels) while nonlinear phase in IIR means distorsion (like a blur) over the image toward the axis in which the filter is applied.
In IP there is not frequency spectrum, you do not seek to have a cutoff at x Hz and a flat passband. In a Butterworth this is the main design focus.
A butterworth or any other SP filter IIR version can be readily be extended to IP. The implementation is not a limitation.
A FIR like the Gaussian and others, dont have ripples in their impulse response. Butterworth, Chebychev, eliptical and ALL other IIR filters have an infinite number of ripples! This is another way to see the distortion which are produced by them.
The butterworth is ripple-less in their frequency response (the spectrum of the inpulse response) which we already agreed it do not have a direct sense in IP.
IIR means infinite impulse response: hence if you apply a Butterworth,  one convolution affects all other axis pixels.
Nevertheless, there are frequency figures on IP. You can define and extend 2D noises distribution based on an spatial frequency, and thus define frequency cutoffs such as the SP case. But this normally solved with a FIR filter, without paying the cost of the distortion an IIR implies.

Which is the significance? Most of their good properties in SP are not applicable in IP, with a major demerit added.


Converting a pitch bend (MIDI) value to a "normal" pitch value

At the core of MIDI is a representation of music as discrete note events, each of those having a static pitch. This is perfect for representing music as played on keyboard instruments. You can convert any frequency corresponding to a note on the tempered scale into a MIDI note number, using:

$69 + 12 \times \log_2 \frac{frequency}{440}$

Under the assumption that the MIDI receiver is calibrated for A4 = 440 Hz.

This representation is alright for piano music, but the problem is how to represent pitches which are not mapped to the tempered scale (non-western music, non musical sounds), and how to represent pitch variations over the duration of a note (glissando, vibrato).

This is done in MIDI by using "pitch bend messages" which instruct the synthesizer to shift the pitch of the currently played note by a small interval. Most synthesizers are calibrated by default for +/- 2 semitones over the course of the pitch bend message range (0 .. 16383). 8192 corresponds to no pitch bending - the emitted pitch is exactly that of the note value. The mapping between the pitch bend value and the frequency shifting ratio is given by:

$\frac{f_{emitted\_note}}{f_{note\_message}} = 2 ^ {\frac{pitchbend - 8192}{4096 \times 12}}$

You can thus get the frequency of a note played by a syntheszier from the following formula:

$440 \times 2^{\frac{note - 69}{12.0} + \frac{pitchbend - 8192}{4096 \times 12}}$

Where note is the 7-bits MIDI note number of the last received Note On message ; and pitchbend is the 14-bits value of the last received Pitch bend message. A synthesizer starts with its pitch bend register set to 8192, and this value is also reset during the reception of a "Reset all controllers" message.

Let us take the following example. You want to express a flute trill with the following frequency trajectory: 500 Hz, 510 Hz, 500 Hz, as MIDI messages.

The base note number is:

$round(69 + 12 \times \log_2(500 / 440)) = 71$.

So you send a "note on" message with note# equal to 71. This is equivalent to a pitch of:

$440 \times 2 ^ {(71 - 69) / 12} = 493.88$

Which is the nearest pitch on the tempered scale. You need to send a pitch bend message to raise the pitch by a factor of:

$\frac{500}{493.88} = 1.0124$

And get your 500 Hz. The corresponding pitch bend value is:

$round(8192 + 4096 \times 12 \times log_2 1.0124) = 9065$

To get your 510 Hz, the pitch bend value would be:

$round(8192 + 4096 \times 12 \times log_2 \frac{510}{493.88}) = 10469$

So your sequence of MIDI messages for 500, 510, 500 Hz would be:

NOTE 71
PITCH BEND 9065
...
PITCH BEND 10469
...
PITCH BEND 9065

You can think of the MIDI note number as the "integral" part of the pitch ; and the pitch bend as a redundant "fractional" part of the pitch.


How to apply Butterworth High Pass Filter

If you just want to do thresholding, then you don't have to use any filter. Just throw away all those samples whose values are below your threshold (in this case, $100 \ \mathrm{pT}$).

In MATLAB, you could do something like this:

data(data < 100) = [];

If you want to keep the array size, then you might consider replacing those values with NaN or 0 instead of [].


understanding discrete-time convolution in LTI systems

You know that $h(n)$ is the system's response to the input $x(n)=\delta(n)$. Because of time-invariance, the response to $x(n)=\delta(n-k)$ is given by $h(n-k)$. If you write the input signal as

$$x(n)=\sum_kx(k)\delta(n-k)\tag{1}$$

then, due to linearity, you can write the output signal as the sum of the responses to each individual term $x(k)\delta(n-k)$ in (1), which is $x(k)h(n-k)$. So you get the convolution sum

$$y(n)=\sum_k x(k)h(n-k)\tag{2}$$

To compute (2) for a given value of $n$ you need to multiply $x(k)$ with $h(n-k)$ for all relevant values of $k$. Since you have $h(n\color{red}{-k})$ in (2) (i.e. with a minus sign before $k$) you need to flip the impulse response (shifted by the current value of $n$), multiply it with $x(k)$ and compute the sum.


Basically you never want to use the Transfer Function representation (with b and a) and rather use the Zeros-Poles-Gain (z,p,k). This will allow you to avoid the numerical errors. In your case you might design your filter in following way:

fs = 44100;           % Sampling frequency
Wp = [30, 70]/(fs/2); % Pass band frequencies (as normalized frequency)
Ws = [20, 90]/(fs/2); % Stop band frequencies
Rp = 3;               % Ripple at pass band
Rs = 50;              % Ripple at stop band

[n, Wn] = buttord(Wp, Ws, Rp, Rs);     % Get order and omega vector
[z, p, k] = butter(n, Wn, 'bandpass'); % Design filter accordingly
[sos, g] = zp2sos(z, p, k);            % Convert to state matrix
Hd = dfilt.df2sos(sos, g);             % Create the filter object

Which for some dummy random signal:

x = rand(1,100000);
y = filter(Hd, x);

Will produce a stable output:

And here is the filter frequency response (everything looks as requested):


Threshold for peak detection of noisy signal in frequency domain

This is a pretty common problem in radar signal processing. Typically, we’ll use what’s called a CFAR detector (constant false alarm rate). This essentially boils down to a non-linear smoothing filter where to attempt to estimate the local average power of the signal (excluding the current cell under test (CUT)). Afterwards, a “bias factor” is applied; this is typically determined a priori. Anything above the CFAR level + bias factor is labeled as a detection, and everything else is ignored. I’m sure there is a wealth of publicly available information on the topic, and I’d highly recommend doing some googling, or even looking around more on this website.

EDIT:
Neyman-Pearson is indeed the standard technique.  There are loads of papers about this topic on IEEE Xplore, I highly recommend you check them out.  That "bias factor" I talked about earlier is the thing you'll want to calculate based on your PFA or PD curve, which is a function of SNR.


Response of Linear System to Stochastic Process

Consider the following LTI system with impulse response $h[n]$

$$ \{v[n]\} \longrightarrow \boxed{H(z)} \longrightarrow \{u[n]\} $$

From the analysis of LTI systems with WSS random inputs, the following results deduced. Given the input $v[n]$ which is zero mean white WSS random process, its ACS (Auto-Correlation Sequence) sequence is:
$$\phi_{vv}[m] = \sigma_v^2 \delta[m]$$ where $\sigma_v^2$ variance of the input.

Then the output ACS is given by:
$$\phi_{uu}[m] = h[m] \star h[-m]^{*} \star \phi_{vv}[m]$$

Similary the input and output PSD (Power Spectral Densities) are related by:
$$S_{uu}(w) = H(w) H(w)^{*} S_{vv}(w) = |H(w)|^2 S_{vv}(w) $$ which follows from the relation that PSD of a WSS RP is the DTFT of its ACS; i.e.,
$$ S_{uu}(w) =\sum_{m=-\infty}^{\infty} \phi_{uu}[m] e^{-j w m} \longleftrightarrow \phi_{uu}[m] = \frac{1}{2\pi} \int_{-\pi}^{\pi} S_{uu}(w) e^{j w m} dw$$

Then the variance of the output is given by the following:
$$ \sigma_{uu}^2 = \phi_{uu}[0] = \frac{1}{2\pi} \int_{-\pi}^{\pi} S_{uu}(w)dw$$

Now using the given information $H(z) = \frac{1}{1-a z^{-1}}$ and $\phi_{vv}[m] = 1 \delta[m]$  you can compute the variance in either way:

$$ \sigma_u^2 = \left( h[m] \star h[-m]^{*} \star \delta[m] \right) |_{m=0} $$ or
$$ \sigma_u^2 = \frac{1}{2\pi} \int_{-\pi}^{\pi} |\frac{1}{1-a e^{-jw}}|^2 dw $$

For example chosing the first approach with $h[n] = a^n u[n]$ yields:
$$
\begin{align}
\sigma_u^2 &= \left( h[m] \star h[-m]^{*} \star \delta[m] \right) |_{m=0} \\
\sigma_u^2 &= \left( h[m] \star h[-m]^{*} \right) |_{m=0} \\
\sigma_u^2 &= \sum_{k=-\infty}^{\infty} h[m-k]h[-k]^{*} |_{m=0} \\
\sigma_u^2 &= \sum_{k=-\infty}^{\infty} h[-k]h[-k]^{*} \\
\sigma_u^2 &= \sum_{k=-\infty}^{\infty} a^{-k} u[-k] (a^{-k})^{*} u[-k]\\
\end{align}
$$

$$ \sigma_u^2 = \sum_{k=-\infty}^{0} (a a^{*})^{-k} = \sum_{k=-\infty}^{0} |a|^{-2k} = \sum_{k=0}^{\infty} (|a|^{2})^k  = \frac{1}{1-|a|^2}$$


It's already answered in the comments so just to close it out: the answers are 3, 7 and 3. The periodicity is basically given by the denominator of the fraction in front of 2*pi*n after the fraction has been simplified (no common divisors between nominator and denominator).  The simplified fractions for the cases are 1/3, 5/7 and 7/3 so the periodicity is 3, 7 and 3 respectively.


Discrete-time Convolution Convergence Issue

In the third edition of Oppenheim and Schafer's Discrete-Time Signal Processing, the corresponding exercise is given with $h[n]=a^nu[-n-1]$, $a>1$ (note the positive sign in the power of $a$). In that case everything works out fine. So either there's a typo in your edition or you made a mistake copying the exercise.


Filtered signal does not follow my cut off frequency

did i make some mistake ?

Yes. The cutoff frequency is specified relative to the Nyquist Frequency and not the sample rate

From the documentation

  Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to half the sample
  rate.

Also: you design a 10th order butterworth and filtfilt() doubles that again to 20. That will give you significant domain "ringing".


Finding Distance between Cameras in Stereo Vision

Calibrate your cameras with the same world coordinate system. Use the same checkerboard and take images at the same time with both cameras.

After you have your extrinsic matrices for each camera, you can use this answer here where you need to solve the equation for the position of each camera center in world coordinate system.

Then you can calculate the distance between two coordinates.


Scaling the input vs scaling the impulse response for an LTI system

Aside from

Input that is only a constant value or
a System that does nothing to the signal, i.e. $H = \text{const.}$,

these things are not the same.
Sketch of a proof:

Scaling input shifts energy to frequency components that hadn't energy before,
Scaling the system means changes the Transfer Function, but keeps the system LTI, which means it can't add energy to any components that had none before.


How to implement a moving average in C without a buffer?

A moving average can be implemented recursively, but for an exact computation of the moving average you have to remember the oldest input sample in the sum (i.e. the a in your example). For a length $N$ moving average you compute:

$$y[n]=\frac{1}{N}\sum_{k=n-N+1}^nx[k]\tag{1}$$

where $y[n]$ is the output signal and $x[n]$ is the input signal. Eq. (1) can be written recursively as

$$y[n]=y[n-1]+\frac{1}{N}(x[n]-x[n-N])\tag{2}$$

So you always need to remember the sample $x[n-N]$ in order to compute (2).

As pointed out by Conrad Turner, you can use an (infinitely long) exponential window instead, which allows you to compute the output only from the past output and the current input:

$$y[n]=\alpha x[n]+(1-\alpha)y[n-1]\tag{3}$$

but this is not a standard (unweighted) moving average but an exponentially weighted moving average, where samples further in the past get a smaller weight, but (at least in theory) you never forget anything (the weights just get smaller and smaller for samples far in the past).


What is meant by "correlation" when referring to spectral coherence

The coherence function, as used in signal processing, measures the normalized correlation between to power spectra:

$$
C_{xy}(f) = \frac{|G_{xy}(f)|^2}{G_{xx}(f) G_{yy}(f)}
$$

where $G_{xx}$ is the power spectral density (PSD) of the signal $x(t)$, $G_{yy}$ is the PSD of $y(t)$, and $G_{xy}$ is the cross-spectral density (CSD) of $x(t)$ and $y(t)$.

The CSD is computed as the Fourier transform of the cross correlation of $x$ and $y$:

$$
R_{xy}(\tau) = \operatorname{E} \left[ x(t) y(t + \tau) \right]
$$
$$
G_{xy}(f) = \mathcal{F} \left[ R_{xy} (\tau) \right]
$$

You can think of the CSD as measuring "hidden periodicities" in the cross correlation, just like the PSD measures periodicities in the autocorrelation. (This is where the name for the "periodogram" PSD estimator comes from.)

Incidentally, if the cross-correlation has harmonic content at some frequency $f$, you can think of the signals being correlated at that frequency.


How to make the impulse response sparse? How does one know that the channel is sparse?

How you parameterize your sparsity will depend on your application.  The authors of that paper, in a paragraph on page 231 say:

which is why they clump the coefficients together in $P$ blocks of start time $t_{B_k}$ of duration $n_{B_k}$.

For one impulse response this is shown in their figure 1(b).

The overall sparsity will depend on the reflections and delays that they have when this impulse response is used by the direct signal and all its indirect reflections.

  How are the impulse response/coefficients made sparse?

You need to decide which coefficients in the overall impulse response are zero. To do that, you need to decide why they are zero (what is the cause?). Perhaps there is a long dead-time between the onset and the first secondary response.

  How does one do parameterization of impulse response to decide for sparse/non-sparse?

Once you've decided on the form, you decide on the structure of the non-sparse parts of the response and how they are separated from each other.

Suppose your non-sparse parts are all defined by $\alpha_k h[n]$, where $\alpha_k$ is the gain of the $k^\mbox{th}$ non-sparse segment and $h[n]$ is the FIR response of all segments (modulo a gain term). $h[n]$ is of duration $H$.  Then the sparsity will come about by how you space out the $\alpha_k h[n]$:
$$
g[n] = \sum_{k=1}^K \alpha_k h[n-t_{B_k}]
$$
where $t_{B_k} \gg H$ for sparsity to be true (so there is lots of space between the $\alpha_k h[n]$.

  In general, how is parameterization applied to make the impulse response to zero and in general how does one know which coefficients are sparse?

This comes down to how do you choose $t_{B_k}$ and $H$. That will depend on what system you're trying to model.


Why are zero values added in the FFT of a concatenated noise signal?

You're running into a property of the DFT that is usually used in the opposite direction: stuffing zeros between samples in one domain results in replication of the entire sequence in the opposite domain. Let's start in the frequency domain with the signal you plotted, $X[k]$, and look at the inverse transform to see what time-domain signal it corresponds to.

$$
x[n] = \sum_{k=0}^{N-1}X[k] e^{j 2 \pi k n/N}
$$

$N$ is the length of the DFT that you apply to the data. Let's say that $X[k]$ consists of sparse values with $M$ zeros between each one, so that $X[k] = 0\ \forall\ k \neq kM$. The above then looks like:

$$
x[n] = \sum_{k=0}^{\left \lfloor{\frac{N-1}{M}}\right\rfloor}X[k] e^{j 2 \pi k M n/N}
$$

Do some rearranging within the exponent to get:

$$
x[n] = \sum_{k=0}^{\left \lfloor{\frac{N-1}{M}}\right\rfloor}X[k] e^{j 2 \pi \frac{k n}{\frac{N}{M}}}
$$

This is nothing but a $\frac{N}{M}$-point inverse DFT of the nonzero samples in your frequency-domain signal $X[k]$. Due to the periodicity property of the DFT, you can then show that $x[n]$ is periodic with period $\frac{N}{M}$, which is exactly the form of signal that you shoved into the DFT to begin with.

Summary: Replication of the signal in one domain results in a zero-stuffed signal in the opposite domain. This is used in the theory behind interpolating filters and other multirate processing. There is a dual property where if you downsample a signal (by removing $M$ out of every $N$ samples), the data in the opposite domain collapses on top of one another and adds (resulting in aliasing).


How to make frequency axis for even and odd FFT length?

You can make a positive frequency spectrum quite simply (where fs is the sampling rate and NFFT is the number of fft bins). In the Matlab implementation of the FFT algorithm, the first element is always the DC component, hence why the array starts from zero. THis is true for odd and even values of NFFT.

%//Calculate frequency axis
df = fs/NFFT;
fAxis = 0:df:(fs-df);

If you need to wrap the frequency spectrum, you need to take into account whether you have an odd numbered NFFT. There always needs to be a DC component, so . .

df = fs/NFFT;
fAxis = (0:df:(fs-df)) - (fs-mod(NFFT,2)*df)/2;

Notice how the calculation of the positive frequency axis is identical to above, but the FFT-shifted term changes to accomodate even or odd FFT lengths.

These code snippets were taken from a long answer posted on SO (that you might find interesting) found here:
https://stackoverflow.com/questions/9694297/matlab-fft-xaxis-limits-messing-up-and-fftshift/9699983#9699983


Comparing similarity of two simultaneously recorded EEGs in freq domain

You can express correlation in the frequency domain but the closest thing to what you are describing would be Coherence.
It returns a number between 0 and 1 per each frequency bin of the spectrum. The closest this number is to 1 the more "similar" the signals would be and you can average this number over the whole spectrum or over a frequency range.
Hope this helps.


Identify Where Singing Starts in a Voice Only Recording

If the background noise is white-ish, you could measure spectral flatness and consider it to be voice when the amplitude is above some threshold and the spectral flatness is below some threshold.

Basically you just take an FFT of a chunk of the signal, then divide the geometric mean of the magnitude of the spectrum by the arithmetic mean.

You could also use a band-pass filter to emphasize only the frequency regions where the human voice usually sits (as simple as setting the unwanted regions of the FFT to 0 before measuring the spectral flatness)


Transforming RGB images including NIR to LAB

Let's start with the source of your RGB + near IR images. To break the process down:

Some camera must have taken those images.
The camera sensor is sensitive to light within two given wavelengths. The technical term for this is spectral response.
Alternating elements (Bayer pattern) in a sensor are etched with filters. Each filter type passes certain wavelengths while supressess the others.
The alternating sensor elements accumulate the filtered light from different wavelengths. This creates our very first raw signal. This raw space is then mapped onto a color space, most commonly the sRGB space.
Filter properties are precisely measured by the manufacturers. They use those properties (and assume a D65 white point, unless manually chosen to something else) to convert the captured raw space into sRGB which gets passed down to you (usually with a little compression and anti-aliasing filtering).

Answering your final question first, with the sRGB assumption in place, you are now armed with one very important source of information - your white point is D65 unless explicitly stated elsewhere. You may now safely proceed to convert your sRGB image into LAB. MATLAB has a handy function for the same.

To answer your questions:

  First question would be: is that correct, or did I just fall for
  single people claiming something?

A: People say a lot of things - largely from their personal experience. You must create your own. I personally exect LAB to help. You can read more about that on my Quora answer here.

  Is there a way to properly include NIR in the LAB color space? Or
  would it be better to keep it as separate information?

All traditional color spaces are designed to fall within the locus of human perception. Since near IR is out of bounds for our perception, it's usually modelled as a standalone entity (seperate channel). However, many imaging devices have a near IR spectral response below 750 nm. If that's the case, I would expect statistical correlation between your G or R channels with near IR. A simple PCA would answer that question.


Discrete-time signals increase in frequency and period

You're right that the periods of the two signals are very different. The period of a sinusoidal sequence
$$x[n]=\cos(n\omega_0+\phi)\tag{1}$$
with $\omega_0=\pi/2$ equals $N=4$, whereas for $\omega_0=999\pi/1000$ it equals $N=2000$.
But this is not what "frequency" refers to here. Note that in general the sequence $x[n]$ as defined in $(1)$ isn't even periodic. It is only periodic if $\omega_0$ satisfies
$$\omega_0=\frac{2\pi p}{q},\qquad p,q\in\mathbb{Z}\tag{2}$$
Maximum frequency for a discrete-time signal is achieved if the samples alternate in sign. This is the case for $\omega_0=\pi$, and it is achieved for almost all samples if $\omega_0$ is close to $\pi$. For $\omega_0=\pi/2$ there are always two samples with the same sign before the sign changes, so the frequency is half of the frequency of a sequence for which consecutive samples have alternating signs.


Design of efficient digital interpolation filter

Butterworth low-pass filters are not going to work for this, as demonstrated in the following. You'd need to use other types of filters, for example linear-phase finite impulse response (FIR).

The magnitude frequency response of an order $N$ discrete-time Butterworth low-pass filter can be approximated by that of an analog prototype (from Wikipedia, with a small notation change):

  $$G(\omega) = {\frac{1}{\sqrt{1+{\omega}^{2N}}}},$$

with normalized frequency $\omega$. The Butterworth filter magnitude frequency response is strictly decreasing. For a pass-band edge gain of $a = 10^{-0.001/20}$, that is, -0.001 dB, we'd employ a further gain factor of 0.001 dB to have a gain of 0.001 dB, that is, $\frac{1}{a}$ at 0 Hz. With this gain factor the magnitude frequency response is $\frac{1}{a}G(\omega)$. Then the transition band begins at the frequency $\omega_0$ found by solving:

$$\frac{1}{a}G(\omega_0) = a$$
$$\Rightarrow \omega_0 = \left(\frac{1 - a^4}{a^4}\right)^{\textstyle\frac{1}{2N}}$$

Then we plot the magnitude frequency response at the beginning of the stop band (the end of the transition band) at frequency $\omega_1 = \frac{0.55}{0.45}\omega_0$ as function of filter order $N$:

Figure 1. $\frac{1}{a}G(\omega_1)$ as function of $N$.

Based on the graph, -174 dB of stop band corner gain would require a filter order of about 120. I don't think such a high order is feasible for a recursive filter as it will lead to numerical problems. Filter types other than Butterworth will be able to meet the specification with a lower filter order.


Book suggestion for upsampling and downsampling

The answer to my question is:
Essentials of Digital Signal Processing 1st Edition
by B. P. Lathi
https://www.amazon.com/gp/product/1107059321/ref=ppx_yo_dt_b_asin_title_o00_s00?ie=UTF8&psc=1
where in Section 6.6 you can find detailed discussions and good examples.


why does my integrator perform noise shaping?

The dominant effect you see is from the differentiation which is a high pass function.

Observe the spectrum after differentiation alone:

As would be clearer on a log log plot the signal is going up at rate f versus frequency. (consistent with a differentiation).

Next observe the spectrum at the moving average output:

This is consistent with expected: The moving average approximates a Sinc filter response, and the nulls of a Sinc function go down at rate f versus frequency. With the differentiator going up at rate f and the moving average going down at rate f we see the flat response as observed with the nulls at $1/T$ where $T$ is the duration of the moving average.

Your final operation that you call an integrator is actually another average combined with a down-sample (decimation), and with that would result in the low frequency portion of your signal after another 16 point average; so another Sinc filter response of similar form and then the decimated (close in) spectrum out to half of the first main lobe:


FIR-Decimation and Low-pass filter (taps vs number of input points vs number of decimation stages)

I'll answer your questions first:

You should look at the final decimation factor when designing the earlier filters. Hence, for the first stage you can have a cutoff frequency of 1Hz and stopband frequency of 100000/10 - 1 = 9999 Hz, since this is what is folded/aliased into your band of interest. You would actually prefer to have several stopbands 9999-10001, 19999-20001, 29999-30001, 39999-40001, and 49999-50000 Hz /you'll have to normalize yourself...). This would require a rather low order filters, which is nice since it runs at the highest sample rate. For the next stage, 999-1001, 1999-2001, and so on. Higher order (about ten times), but ten times lower sample rate. And so on.
Roughly more stages => fewer multiplications per input sample (given that you do a proper design as above). This holds up to some limit which can be quite straightforwardly determined by searching some combinations.
After the first stage you will have $1000 + 2 \times 127$ samples (with 128 taps) and so on.

If you want to get the DC value you can just average the sequence. This is, apart from the possibly obvious observation that the DC level should be the average "offset" of the otherwise AC signal, also what the DC component of the DFT evaluates to.


It depends on your specification really.

When you're decimating, some degree of aliasing will occur. Your tolerance of this aliasing, signal delay, and changes in the phase response of your input signal will determine the best approach to take.

A very long FIR might provide a very steep cutoff and linear phase performance (depending on design), but will delay the signal (group delay), may cause pre-ringing and can be computationally demanding. FIRs however also allow for elegant 'polyphase' implementations (which reduce computational demand by effectively breaking the filter down into smaller component filters).

IIRs on the other hand have the potential to have a shorter delay, at the expense of a non linear phase response and perhaps a more shallow roll-off. You could achieve a steep rolloff with an IIR but it would likely be quite high order and hence susceptible to coefficient quantisation and/or stability problems.


Complex Samples (IQ) - Baseband Filtering

A real signal is complex conjugate symmetric, so any filter over $-fs/2$ to $+fs/2$ is really only unique over the range from $0$ to $fs/2$. This would be a primary motivation for working with a fully complex signal as the unique range is then truly extended over the full $-fs/2$ to $+fs/2$ range.
Another more dominant reason to implement a complex filter besides the bandwidth requirement outlined above (where the choice is really sample twice as much or carry two datapaths; similar complexity and certainly you could map either to be similar if that was the only goal) is when the passband itself is not symmetric: specifically as a bandpass filter when the positive frequency passband is not equal and magnitude and conjugate in phase to the negative frequency, or as a baseband filter- the passband shape itself is asymmetric or perhaps completely one sided. Another dominant reason is when real signal conditions would create images very close to the signal where further filtering would be difficult, such as if conditions require operating very close to the Nyquist boundary.  If the OP doesn't have these conditions, I don't see a strong motivation for complex filtering, which requires 4 real multipliers for every complex multiplication operation (assuming complex inputs and outputs).
Channel equalization is another example where a full complex filter is often necessary (as an implementation on a baseband signal) but not needed for a passband signal.


Designing the filtering for a time domain data acquisition application

Is this a good filter for my task, or would you recommend something else?

It might be given that you do not need to get rid of any other noise/signal in that band. The advantage will be the lack of multiplications and therefore most likely fewer cycles in your implementation.

  Should I just send the latest sample, or perhaps take the average of all samples between two "PC sample points"?

Taking the average is just like using another moving average filter, so this may or may not be a good idea. It is easily analyzed in a similar way. Given that your signal is "good enough" to downsample, filtering again will, while removing a bit of noise, primarily mess up with the passband.

It is worth noting though that it may be worthwhile to use a so called polyphase decomposition. You do not really want to spend time computing intermediate samples which are not transferred to the PC anyway. Hence, you do in fact have quite a bit more cycles per sample. The straightforward way to go about this to use an (or several) FIR filter. In the case of several compute the convolution of all. Only apply the filter such that it computes exactly the output you want. If you go with this approach you might just as well design a single FIR filter using general multipliers (since the convolution of the moving average filters will have general coefficients anyway).

Assuming you have the signal package for Octave (and that the call is identical to Matlab), you can try something like:

fs=1000
fn = fs/2 % To simplify things
delta = 3 % range in Hz around the power line frequencies, should be less than 5...
eps = 1e-7 % Some small number....
N = 50 % Filter order
W = 10 % How much more you'd like to attenuate the power compared to generic noise

h = firpm(N, [0 5 (50-delta) (50+delta) (50+delta+eps) (60-delta-eps)  (60-delta)  (60+delta) (60+delta+eps) fn]/fn,[1 1 0 0 0 0 0 0 0 0], [1 W 1 W 1])

The lesser you pick delta, the more likely that you will get a notch exactly where you want it...


Understanding Arrangement of second order sections

Different orderings of second-order blocks have different behavior with respect to quantization noise and overflow. Usually it is best to avoid very peaky transfer functions in second order blocks, because these peaks may cause overflow and quantization noise amplification. This is why it is common to pair poles and zeros that are close to each other, so the zeros cancel part of the peaks caused by the poles. The pole-zero pairing determines the second-order sections, but not their order. Concerning the order, there is always a best choice minimizing the quantization noise, but this order needs to be determined individually for the filter under consideration. However, there are two types of ordering that work well in most cases. Depending on the way scaling is performed, either ordering the sections with increasing or decreasing pole radii is best. In Matlab, these two common ways of ordering are referred to as 'up' and 'down', respectively. Have a look at this document for more detailed information.


How can I find the amplitude of the 20 Hz frequency component of a signal faster?

As pichenettes already pointed out in a comment above, one approach is to use the Goertzel algorithm. However, it is a somewhat common misconception that in order to deduce the amplitude of various frequency bands in a signal, one must use Fourier-transform-based techniques. There are other methods that may be more appropriate to your application.

In your case, for example, you would do equally well by just creating any bandpass filter centered at 20 Hz and watching the amplitude of its output over time. The DFT (and by extension the Goertzel algorithm) can be viewed as just a bank of equally-spaced, critically-sampled bandpass filters. If you only need one output, then a DFT is overkill, and the Goertzel algorithm can be limiting, as you don't have much control over its frequency response about 20 Hz.


I just did some research on existing implementations of the MPEG-7 audio descriptors and found the following implementations:

The official implementation of the 17 low-level descriptors (part 1 and part 2) (MATLAB).

Alghough the descriptors are implemented in MATLAB, it should not be a big issue to implement them in Python, based on the MATLAB code. Maybe the following resources are also of interest:

The Timbre Toolbox, which implements a number of features described in this JASA article, which include some of the MPEG-7 audio descriptors (MATLAB) and
An an online computation of all 17 audio descriptors with the MPEG-7 Low Level Descriptor Extractor of the TU Berlin.


Can we choose a sampling frequency to remove unwanted noise at a specific frequency?

Or another way, visually, with the plot of the magnitude response of your filter shown below you see that it goes to zero at $\pi/2$.

So what you could do is use the relationship between the angular frequency $\omega$ in [radians/sample] and the relative frequency $f$ in [cycles/sample]
$$
\omega = 2\pi f = 2\pi \left(\frac FF_s\right)
$$
Then compute $F_s$ for $F = 80 \ \rm Hz$ and $\omega = \pi/2$. You then get
\begin{align}
\frac \pi 2 &= 2\pi \left(\frac{80}{F_s}\right)\\
\implies \frac 12 &= 2\cdot \left(\frac{80}{F_s}\right)  \\
\iff 80 &= \frac 12\cdot \left(\frac{F_s}{2}\right)
\end{align}
Giving you the $320\ \rm Hz$.
In MATLAB
You can plot the magnitude response from the frequency response using MATLAB's freqz function where you can specify the frequency grid from $0$ to Nyquist (i.e. $\pi$). With your resulting $320\ \rm Hz$, you can run the function with this $F_s$ and see that $80\ \rm Hz$ is indeed suppressed.


Bandlimiting and OFDM spectrum - what is considered "bandlimited"?

This isn't a question specific to OFDM:

All transmissions need to be bandlimited by law. Hence, law defines band-limited here.

Typically, this happens in the shape of spectrum masks (often expressed in textual form in regulations); these say something like (this is just an example, actual legislation varies across regions and is definitely different):

"not more than -40 dBm/Hz within 10 MHz of the carrier frequency, not more than -90 dBm/Hz between 10 and 20 MHz of the carrier frequency, and -120 dBm/Hz more than 20 MHz away"

Regarding the specialties of OFDM:

You indeed can't pulse shape OFDM without losing its benefits. (You can do a bit of a tradeoff, and trade a bit of orthogonality for out-of-band radiation suppression, but you'd typically avoid that)

That's why you typically have "empty" carriers at the edges. You simply leave room for the sinc to decrease.


If I read the plot correctly, I see tones repeating every 1 MHz, 3MHz, 5MHz, etc under the envelope of a $Sinc^2$ function with nulls every 200 MHz and an overall sampling rate of 2 GHz (You further clarified 10 GHz, so the full digital spectrum must not be shown?). That envelope can be the result of a 10 ns triangular pulse, and given the repetition rate, consistent with the plus minus pattern you described repeating every 1MHz.

Please include a zoom-in of the time domain of your waveform, specifically over a 50 ns time span showing the discrete samples. This will further clarify if this theory is correct of if you are also being affected by a zero-order hold.

Here is a plot of showing the frequency response of [1 1 1 1 1 1 1 1 1 1] convolved with itself, thus creating a triangular pulse 10 ns wide at a 2 GHz sampling rate), which is a $Sinc^2$.  The same result can be achieved at the 10 GHz sampling rate; being simply discrete samples of a 5 ns rectangular pulse - which result in a Sinc function in frequency with nulls at 200MHz spacing; since we see from the frequency roll-off in your plot what would be a $Sinc^2$ function, then two rectangular pulses each 5 ns wide convovlved in time (resulting in 10 ns triangular pulses) would then give the $Sinc^2$ function achieved:

Generally, you can determine the Fourier Transform of any repeating pattern by using the Fourier Transform of the base pattern as the envelope (in your case the triangle waveform which results in the $Sinc^2$ response as the envelope) and then within that envelope, due to the signal repetition, the only non-zero frequency components that can only exist will be at integer multiples of the repetition rate of the base pattern (I explain this in detail in another response and will update this with the link to that if I can find it). In your case it is further complicated by a positive triangle wave followed by a a negative triangle wave as your "base pattern" that repeats at the 1 MHz rate. To get the resulting envelope in this case, this can be easily conceptualized by first passing your triangle wave through a comb filter response as shown in the plot below (The comb filter is a delay of 500 ns and a subtraction resulting in nulls at DC, 2 MHz, 4 MHz etc).  Thus the original $Sinc^2$ envelope now has nulls at DC, 2 MHz, 4 MHz. When you repeat this base pattern at the 1 MHz rate, this would result in tones every 1 MHz, except the even ones are nulled due to the effect of the comb filter.

In the plot below, this process is shown both in the time domain and the associated response in the frequency domain for each portion of this process as previously described.

So with all that said, as long as your triangle pulses are 10 ns and not the 50 ns you first stated, the FFT you show is how it should look and not a problem with your settings. However you may want to consider windowing your signal prior to taking the FFT to reduce the sidelobes shown on each of the tones (that are an artifact of the FFT when windowing is not used). This will reduce the width of each "tone" but substantially reduce the sidelobes.


What are some alternatives to Spectrogram charts?

The ltfatpy 1.0.16 package

  is a partial Python port of the Large Time/Frequency Analysis Toolbox
  (LTFAT), a MATLAB®/Octave toolbox for working with time-frequency
  analysis and synthesis.

Among linear and quadratic time-frequency methods, there is a large number of options for sharper analysis tools converting a 1D signal into 2D data. You can even get more information by decomposing onto 3D or 4D spaces (coming back on that later). Quadratic (or higher-order) representations are usually more difficult to invert than linear ones.

And here is a documentation for ltfatpy.


Influence of image size to edge detection in OpenCV

You should play around with the many parameters of the Canny() and HoughLinesP() functions before resorting to changing the image size.

The prototype for Canny() is:

void Canny(InputArray image, OutputArray edges,
           double threshold1, double threshold2,
           int apertureSize=3, bool L2gradient=false )

You should play with the threshold1, threshold2, aperturSize and L2gradient parameters until you get an output from Canny() that accurately marks the pixel involved in the lines in your image.  The threshold2 parameter controls the intensity above which a gradient value will be considered definitely part of an edge.  Setting it higher will give you fewer edges.  The threshold1 parameter controls the intensity above which a pixel will be considered to continue an edge from a neighboring pixel.  Setting it lower will give you longer edges.  The ratio threshold2/threshold1 should be at least 2.0 and can be as large as 3.0, and maybe larger.  Setting L2gradient to true will make the algorithm slightly slower but will use a more accurate gradient calculation from the first derivatives.  Setting apertureSize to 5 or 7 might smooth out some noise in the first derivative calculation.

Only once you have Canny() accurately finding the edges in your image should you move on to running HoughLinesP().  HoughLinesP() also has a bunch of knobs that you need to tune to get decent results.

void HoughLinesP(InputArray image, OutputArray lines,
                 double rho, double theta,
                 int threshold,
                 double minLineLength=0,
                 double maxLineGap=0 )

rho and theta can be set larger and smaller depending on how many lengths and angles you want to allow.  Setting rho and theta larger will give you fewer line segments, but they will be less accurate.  threshold can be increased to reduce the number of line segments that you get.  minLineLength can be increased to reject short line segments and maxLineGap can be increased to make it easier for short lines to get connected together.

Setting Canny()'s apertureSize parameter larger is somewhat like working with a smaller image, because a larger aperture is doing more blurring/smoothing of the details of the image.  Setting HoughLinesP()'s rho and theta parameters larger is somewhat like working with a smaller image because you are lumping together small differences in angle or length that would be indistinguishable in the smaller image.

The OpenCV documentation for HoughLinesP() has an example with:

Canny( src, dst, 50, 200, 3 );
...
HoughLinesP( dst, lines, 1, CV_PI/180, 80, 30, 10 );


Difference between PSD estimate and variance of DFT

So, as usual, I start by pointing out what the PSD actually is:

The power spectral density is the Fourier Transform of the autocorrelation function of a wide-sense stationary random process.

  Always assuming  $\hat{X}=0$, is it generally correct  to interpret the PSD as an biased estimate of the DFT variance? If not, can someone  explain  the difference, please?

The Wiener-Chintschin¹ theorem points out that this happens to coincide with the expectation value of the absolute of the Fourier transform of a signal.  In other words, to answer your first question:

The PSD is a stochastic (hence, "invisible") property of a random process. The DFT can be used to estimate that; not the other way around!

To answer your second question, a bit of juggling the term of variance is required:

  Would Bartlett's method become an unbiased estimator for the DFT variance if (3) was scaled by $\frac{1}{L-1}$ instead of $\frac{1}{L}$?

Variance is formally defined as the expectation of the square of the deviation of any observation to the mean:

$$Var(X) = \mathbf E\left[\left(X-\mu\right)^2\right]$$

Now, applying Parseval's theorem, the power in a time and frequency domain should be identical. Considering a white random process, where the covariance of two samples is 0, it's easy to see that variance is kind-of-linear:

$$Var\left(\sum\limits_{i=1}^L X_i\right) = \sum\limits_{i=1}^N Var(X_i)$$

Considering the special case of L=1, it becomes clear that $(4)$ isn't applicable here – you kind of get the standard deviation of the process, but not the variance.

¹ there's more latin-lettered spellings to Chintschin than there's letters in in his cyrillic name...


DFT conjugate of $X^*[k]$, how to prove its formula in terms of $x^*[n]$?

So you want to show $X^*[k]$ ?

Actually it follows very simply from the fundamental properties

$x[n] \longleftrightarrow X[k]$
$x[-n] \longleftrightarrow X[-k]$
$x[n]^* \longleftrightarrow X^*[-k]$

and by combining the last two you have:

$x[-n]^* \longleftrightarrow X^*[k]$

which shows that $X^*[k]$ is given by the forward DFT of the signal $x[-n]^*$.

Note that, due to periodicity of DFT sequences, the negative index can be repaced by $-n = N-n$, and hence $x[-n]^* = x[N-n]^*$ too.

Direct derivation is as follows:

$$ X[k] = \sum_{n=0}^{N-1} x[n] W_N^{nk} $$

where $W_N= e^{-j \frac{2\pi}{N} }$.

Take conjugate of both sides:
$$ X^*[k] = (\sum_{n=0}^{N-1} x[n] W_N^{nk})^* $$

$$ X^*[k] = \sum_{n=0}^{N-1} x^*[n] W_N^{-nk} $$

replace $n$ with $-n$ so that forward DFT appaers:

$$ X^*[k] = \sum_{n=0}^{N-1} x^*[-n] W_N^{nk} $$

which again shows that $X^*[k]$ is given by DFT of $x[-n]^*$. Note again that periodicity of $x[n]$ is utilized to rearrange the limits.


RF-Chain Signal Delay for Sensor Switching

The OP clarified in comments under the question that the intention is to compensate for group delay distortion introduced in the hardware. The typical approach to optimally compensate for this (in a least squared sense) with processing alone is to use the Wiener-Hopf equations to determine the coefficients of an equalizer that can be implemented with an FIR filter (meaning a difference equation with only feed-forward terms). I detail the full approach of doing this in these posts so will provide the links below, but to bottom line the process, the channel to be equalized (the receiver) is "sounded" with a known waveform that is spectrally rich (pseudo-random noise or frequency chirps are great choices as they also offer high average power or SNR; an impulse is a poor choice since it is a challenge to do that with high SNR), and with that and the received signal after the channel the reverse deconvolution can be computed in a least squared sense to determine the effective inverse channel, but importantly applicable to mixed phase systems (systems with both leading and trailing echoes) which on their own can't be inverted due to having zeros in the right-half plane. This approach would be ideally suited for distortions introduced in hardware implementations that are not changing with time (within an acceptable tolerance) since the channel can be sounded once and fixed coefficients determined relatively easily with pre-processing that can then be used without further modification, in contrast to the iterative algorithms that are needed to provide an adaptive solution for time-changing channels.
Details of the Wiener-Hopf Equations and shows application to determine the transfer function of the channel:
Compensating Loudspeaker frequency response in an audio signal
Equalizer Implementation Example by swapping Tx and Rx in the previous case, we can instead solve for the causal equalizer for a mixed phase system instead of the channel itself:
How determine the delay in my signal practically


DSP filter. Is prewarping performed when discretizing using Forward/Backward?

(After some helpful comments, I think I see it a bit more clearly. Here is an attempt at answering my own question, please let me know if I am completely off track.)

Answer:

In my question I must have misunderstood from where the top expression $H(z)$ comes from.
I thought that such an expression was the result of a direct derivation from an analog circuit using the bilinear transform.

As Robert comments to the original question, $H(z)$ as I've written it, is the transfer function for the Hal Chamberlin SVF. And is not the result of a direct derivation from an analog transfer function using the bilinear transform (also pointed out by Matt).

Also, to answer my other question as far as I know there is no pre-warping done for Euler forward/backward and the term $\sin(\pi f_c T)$ does not come up because of this.

Where does the $\sin(\pi f_c T)$ term come from?

Bluntly, I think from an approximation. An approximation that is baked into to what we call the Hal Chamberlin SVF (presented in "Musical Applications of Microprocessors").

By using Roberts fantastic EQ cookbook, which has transfer functions and expressions for the coefficients, but which are derived from analog prototypes, I can compare and see where the approximation is being made.

(I received a nice guide that walks through the Hal Chamberlin SVF its transfer function and how it compares to transfer functions derived using the Bilinear transform and later the Impulse Invariant mapping. Link here: "Exact Bilinear Transform or Impulse Invariant mappings of Q and resonant frequency from analog prototype to Hal Chamberlinʼs State Variable Filter (SVF)" by Robert Bristow-Johnson. I'll also paraphrase the part which explains the $\sin{()}$ term.)

Starting from (note: not derived via bilinear transform)

$$
H(z) = \frac{k_f^2 \, z^{-1}}{1-(2-k_f(k_f+k_q))z^{-1} + (1-k_f k_q)z^{-2}} \, ,
$$

I can compare the coefficients in the denominator with the coefficients that were derived from an analog lowpass biquad which are presented in the EQ cookbook. From this we see that

$$
  \frac{a_2}{a_0} = \frac{1 - \frac{\sin(\omega_c)}{2 Q}}{1+\frac{\sin(\omega_c)}{2 Q}} = 1 - k_f k_q
$$
$$
  \frac{a_1}{a_0} = \frac{ - 2 \cos(\omega_c)}{1+\frac{\sin(\omega_c)}{2 Q}} = k_f^2 + k_f k_q - 2
$$

where $\omega_c = 2 \pi f_c T$.

Using these two and solving for $k_f$ gives us,

$$
k_f = \frac{2 \sin(\frac{\omega_c}{2})}{\sqrt{1+\frac{\sin(\omega_c)}{2 Q}}} \, ,
$$

If we assume $Q > \frac{1}{2}$ and that $\omega_c \ll \pi$ (cutoff frequency is small compared to Nyquist) we can approximate,

$$
k_f \approx 2 \sin\left(\frac{\omega_c}{2} \right) = 2 \sin(\pi f_c T) \, ,
$$

and we have the $\sin(\cdot)$ term I was looking for.


Good for what purpose?  There are always trade-offs.

You can play with the color palette or function, or do local non-linear color selection.

You can try compositing/stacking multiple resolutions of spectrograms (both time and frequency) for higher contrast and perhaps more interesting appearances, but that process can add its own artifacts.


the plots are likely from a line or what is often called a ULA and there is a mirror image response.  It is sometimes called a left/right ambiguity.

Your question’s answer depends on the purpose of your simulation.  The plots are probably equivalent, the response is over 360 degrees in both plots but people often drop the redundant half circle.

if you have a problem where objects are distributed over 360 degrees, simulate that, otherwise 180 is fine.

I usually plot in $\cos(\theta)$ coordinates which range from -1 to 1 for most of my work. For a ULA, the pattern is the same regardless of the steering direction using this convention. If the array weights are symmetric, 0 to 1 is sufficient.

It all depends on what your application is and what you are investigating and the conventions of your community which isn’t specified in your question.

The only caveat is that in some applications you might have some structure blocking the back half of the array. In SONAR it is called a baffle. In RF there are other caveats like a ground plane on one side.

Actual arrays have more complicated responses because what you showed are probably Array Factor responses where each element is assumed ideally “isotropic”.  Real array elements have some directionality. The Array Factor also ignores element to element coupling.

In summary, it could be warranted but in many applications it might not matter.


The Hough transform seeks to replace finding colinear points by  finding concurrent lines.  The two problems are  mathematically equivalent. But dealing with all discrete pixels and a large number of digital lines is complicated and computationally expensive.
Most implementations of the Hough transform discretize the set of distances and angles. With respect to some origin, one can describe a line in polar coordinates by:

a positive distance and an angle in $[0,\,2\pi[$
a signed distance and an angle in $[0,\,\pi[$

This is explained in Polar coordinate system/Uniqueness of polar coordinates:

any polar coordinate is identical to the coordinate with the negative
radial component and the opposite direction

History: the rho-theta parametrization universally used today was first described in Duda, R.O.; Hart, P. E. (1972). "Use of the Hough Transformation to Detect Lines and Curves in Pictures".


What is the betwwen [i,j] and [u,v] in explinations of correlation and convolution (picture attached)?

$[i,j]$ are your spatial coordinates and $[u,v]$ are summation indices.

The difference between correlation and convolution is simply the sign: correlation has $i+u$ and convolution has $i-u$.

They are close cousins: it's easy to express one as a function of the other. For example: correlation is the same as convolution with a properly mirrored/flipped version of the kernel. If the kernel is sufficiently symmetric, both correlation and convolution will yield the same result.


How to estimate the period of a pulse train?

If all you need to estimate is the period of your pulse train then you can threshold the original waveform and work on the clean one by detecting the rising and falling edges.

An FFT based approach would return an average estimate of its "frame of observation" anyway whose accuracy would depend on the length of the frame. (that is, an average estimate of the buffer that is passed at its input)

If you are concerned about the noise, use a short moving average filter (for instance a length 3 or length 5) to smooth the curve prior to estimating and applying the threshold. If absolute timings are essential for your application don't forget to take into account the group delay of the filter.


I wrote this program long ago and published, but removed it later, because could be an abuse. It is possible to write a plugin for the browser, that extracts hidden images and show them in the browser in place of originals, so some visitors will be seeing art gallery and others, who downloaded plug in, something else, may be illegal. I don't mean regular pornography, it is legal. It may be used by terrorists.
The way I did it was follows: I compressed one image and redistributed it in two least significant bits of each byte of pixel data. I used my own compression algorithm. In steganography the compressed image is called payload and where it is embedded is called carrier. My payload was half of the area of carrier. And I hid even multiple images within each other like Russian dolls. The carrier can be compressed and hidden within new carrier. However each payload was half of carrier area. I tested those program that claimed to determine steganographic content. Neither of them detected payload. However, there is a way to find that out. The carrier must be lossless, I used PNG. And that raised suspicion. Why image is posted in PNG format. Steganography is not illegal, I decided to remove my program from internet myself, there was no pressure.


In other words, $h[n]$ is ideal but its frequency response $H(e^{jω})$ is practical and realistic.

This is not true. To compute $H(e^{j\omega})$, you must evaluate the summation $\sum_{n=-\infty}^{\infty}h[n]e^{-j\omega n}$. A real system can't do this because $h[n]$ is not finite length.
Say you choose to use a rectangular window $w[n]$ (with length $2M+1$) to force $h[n]$ to a finite length and the summation becomes $\sum_{n=-\infty}^{\infty} h[n]w[n]e^{-j\omega n} =\sum_{n=-M}^M h[n] e^{-j\omega n}$. There is still the problem of $\omega$ being a continuous variable. For example, a computer can represent $h[n]w[n]$ in an array, but how is a computer supposed to store $H(e^{j\omega})$ in memory? It should have to store all function values between $-\pi$ and $\pi$ of which there are infinite.