Estimating the Nyquist rate when simulating data

You must know the method of computing Q(t), and if Q(t) = Q[f(t)], meaning that Q is a function of another function, like the output of a filter. Then, knowing the autocorrelation or power spectrum of f(t), one can use various methods to find the power spectrum of Q. If the function Q is linear with respect to f(t), then the solutions are simple.

However, it seems like you should clear up your needs. An FFT of the signal is simply its spectrum, not is power spectrum or its power spectral density. They are slightly different.

To answer your last question, it seems like you are formulating this in a probability/statistics sense when you say independent? If they are truly independent (i.e. one value of Q gives absolutely no clue about the next value), then your power spectral density represents white noise and it is flat (a constant). Again, it would be best if you clarified your definition of "independent". For example, if Q is a relatively smooth function of time, then it is not statistically "independent" as you say, and to find the spectrum or power spectrum density, you should analytically look at other functions that Q depends on first.


Any guidance on converting Time Domain to Frequency Domain in Matlab

@Nimit. I recommend that you study and learn to design and implement what are called 2nd-order IIR lowpass filters. Once you thoroughly understand such filters you should design a 2nd-order IIR lowpass filter, that has unity gain at zero Hz, and experiment by applying your noisy signal to that filter and see how much noise reduction is achieved. Then you should apply your noisy signal to a cascade of two identical 2nd-order IIR lowpass unity-gain filters and see how much noise reduction is achieved. Finally you should apply your noisy signal to a cascade of three identical 2nd-order IIR lowpass unity-gain filters and see how much noise reduction is achieved. Nimit, you have much experimenting to do.


"Does this mean discrete time systems are clocked so that the values only at discrete time are available to the system?"

yes, but the system can do some mathematics to interpolate values in between the discrete times of the samples.

"How are these discrete time systems implemented physically?"

well, it's called "electronics" and the discipline to learn is "electrical engineering".

the most essential part is a computer of some form.  DSP chips are little computers of a form that are optimized to perform the most common DSP algorithms.  you will likely also need something like an A/D converter and D/A converter, but without more specific information about what the discrete-time system is supposed to do, i wouldn't know that for sure.


Help me with BM3D algorithm's aggregation scheme

BM3D intuition

The intuition behind BM3D is the same as other non-local algorithms: an image is made by the juxtaposition of many, self-repetitive small image elements called patches. Thus, when looking at small enough pieces, it is possible to find enough similar patches that can be combined and denoised together.

The different non-local algorithms are different in the way used to combine these similar patches.
BM3D relies on sparsity. This principle has become widely accepted now, and basically states that for many useful signals (including images) there exists at least one frame where this signal can be decomposed into a small number of meaningful components.

For images, DCT and wavelet basis are obvious candidates for the sparse frame: JPEG / JPEG2k compression rely on the transformation of the input images into 8-by-8 blocks followed by retaining the strongest coefficients, which is indeed a sparsifying transform.

When noise is added to a signal, its frequencies get mixed or added to the frequencies of the noise-free signal. For example, white noise adds non-zero values to the whole spectrum. Thus, finding a sparse approximate of a noisy signal (in a good basis of course) yields usually a good denoising performance.

BM3D mechanics

The 3D transform of BM3D is actually a 2D+1 transform.

Algorithmically, this is what happens (I'll describe the mathematical intuition afterwards):

make a stack of similar patches. All the patches of the stack are supposed to be small variants of the same ideal (noise free) patch. The variant term is important, you'll see why later;
take a 1D transform (DCT / wavelet) along one of the horizontal axes;
take a 1D transform (DCT / wavelet) along the other. What you have now is a stack of transformed patches;
now, take the transform of choice (DCT / wavelet / Wiener filtering) along the vertical axis. What you are doing here is taking a transform (DCT / wavelet) or directly computing a denoised value (Wiener filtering) for each point in the vertical column;
(optional) if you used DCT / wavelet transform, apply a sparsifying transform, usually hard thresholding.
Then, compute the denied patches by inverting the different transforms (in reverse order);
now, back project the result to the original patch locations. Patches can (and will) overlap when you denies the whole image. This raises an important practical question: how do you fuse the results on overlapping areas? Which denoised patch is the most reliable? BM3D uses Eq. 5 to fuse the results by a weighted averaging: the weight of a back-projected patch is inversely proportional to the variance in the patch stack. When the variance of the stack is high, it means that the patches did not really come from the same ideal one, and thus the result is marked as less reliable.

Mathematically, the intuition is that most patches should be sparse in some good domain (typically, DCT domain for $8 \times 8$ patches).
BM3D also supposes that patches in a stack are variants of a same patch, and thus the variations between these transformed patches and the ideal patch (in the transformed domain) should also be sparse. This assumption is enforced by the 3rd transform and the hard thresholding step.


Normalization of signal against reference

I assume your reference image data y_reference is acquired when there is no light on the image sensor - that's the standard thing to do. Then the usual normalized signal is
y_normalized = (y_signal - y_reference) / y_max
where y_signal is the signal image data and y_max is the maximum possible pixel output. y_max can be measured by increasing the light level on the image sensor until the output no longer increases.


Implementation of Block Orthogonal Matching Pursuit (BOMP) Algorithm - Fix Given Code

Answer taken from Implementation of Block Orthogonal Matching Pursuit (BOMP) Algorithm.

The Block Orthogonal Matching Pursuit (BOMP) Algorithm is basically the Orthogonal Matching Pursuit (OMP) Algorithm with single major difference - Instead of selecting single index which maximizes the correlation we chose a set of indices, sub set of columns of the matrix and the solution vector.

A good reference for the algorithm is given in:

An Optimal Condition for the Block Orthogonal Matching Pursuit Algorithm.
Block Sparsity: Coherence and Efficient Recovery.

The code is given by:

function [ vX ] = SolveLsL0Bomp( mA, vB, numBlocks, paramK, tolVal )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveLsL0Omp( mA, vB, paramK, tolVal )
% Minimizes Least Squares of Linear System with L0 Constraint Using
% Block Orthogonal Matching Pursuit (OMP) Method.
% \arg \min_{x} {\left\| A x - b \right\|}_{2}^{2} subject to {\left\| x
% \right\|}_{2, 0} \leq K
% Input:
%   - mA                -   Input Matirx.
%                           The model matrix (Fat Matrix). Assumed to be
%                           normlaized. Namely norm(mA(:, ii)) = 1 for any
%                           ii.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numBlocks         -   Number of Blocks.
%                           The number of blocks in the problem structure.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - paramK            -   Parameter K.
%                           The L0 constraint parameter. Basically the
%                           maximal number of active blocks in the
%                           solution.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - tolVal            -   Tolerance Value.
%                           Tolerance value for equality of the Linear
%                           System.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range [0, inf).
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  An Optimal Condition for the Block Orthogonal Matching Pursuit
%       Algorithm - https://ieeexplore.ieee.org/document/8404118.
%   2.  Block Sparsity: Coherence and Efficient Recovery - https://ieeexplore.ieee.org/document/4960226.
% Remarks:
%   1.  The algorithm assumes 'mA' is normalized (Each column).
%   2.  The number of columns in matrix 'mA' must be an integer
%       multiplication of the number of blocks.
%   3.  For 'numBlocks = numColumns' (Equivalent of 'numElmBlock = 1') the
%       algorithm becomes the classic OMP.
% Known Issues:
%   1.  A
% TODO:
%   1.  Pre Process 'mA' by normalizing its columns.
% Release Notes:
%   -   1.0.000     19/08/2019
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

numRows = size(mA, 1);
numCols = size(mA, 2);

numElmBlock = numCols / numBlocks;
if(round(numElmBlock) ~= numElmBlock)
    error('Number of Blocks Doesn''t Match Size of Arrays');
end

vActiveIdx      = false([numCols, 1]);
vR              = vB;
vX              = zeros([numCols, 1]);
activeBlckIdx   = [];

for ii = 1:paramK

    maxCorr         = 0;

    for jj = 1:numBlocks
        vBlockIdx = (((jj - 1) * numElmBlock) + 1):(jj * numElmBlock);

        currCorr = abs(mA(:, vBlockIdx).' * vR);
        if(currCorr > maxCorr)
            activeBlckIdx = jj;
            maxCorr = currCorr;
        end
    end

    vBlockIdx = (((activeBlckIdx - 1) * numElmBlock) + 1):(activeBlckIdx * numElmBlock);
    vActiveIdx(vBlockIdx) = true();

    vX(vActiveIdx) = mA(:, vActiveIdx) \ vB;
    vR = vB - (mA(:, vActiveIdx) * vX(vActiveIdx));

    resNorm = norm(vR);

    if(resNorm < tolVal)
        break;
    end

end

end


Nyquist Frequency on semi-unevenly sampled data

the lomb scargle in matlab, (plomb) returns a frequency vector as the second output. i would be surprised if the routine you are using doesn’t do the same.  there is a common tendency for python signal processing libraries to be functionally equivalent to matlab conventions.

your plots show symmetry around a “center” frequency which is the same as what a FFT does for real inputs, except a FFT would not extend to start a third interval.  The high half frequencies(above $\approx$ .017 to .033)  are the negative frequencies in FFT speak. you can ignore them when you have a real signal.

There really isn’t a simple “Nyquist” frequency with nonuniform sampling.

The documentation of the library you are using should detail how frequency is treated. If it doesn’t, you can find a routine that does or do your own calibration by using test sines of known frequency.


how can I generate SIMO underwater acoustic channel

While the literature is small, the Oceans are large.  There is a large variation in physical detail.  Unlike RF, there are no ITU standard models. A serious effort really requires real data, You can spend 20 years refining a model.

You could start by looking at,

  Stojanovic, Milica, and James Preisig. "Underwater acoustic
  communication channels: Propagation models and statistical
  characterization." IEEE Communications Magazine 47.1 (2009): 84-89.

At 100Khz, a ray tracing model should be adequate and you can get ray models at Porter's Ocean Acoustic Library at:

http://oalib.hlsresearch.com/

Ray shooting can require a lot of trial and error to get eigenrays. At short ranges, an iso-velocity range independent model might be adequate.  Most models are Fortran, not Matlab, but you can call them from Matlab. If you look around the web, particularly those of Australian establishments , there is code with tighter MATALB integration.  Most prop code is deterministic.  Ira Dyer is a pioneer is random propagation, so papers that cite him, might be useful, and you can do that using Google Scholar. If you can read Russian, a lot of their literature is very good.

Noise models are local but at 100KHz, noise sources are going to be near.

Projector and forward volume scattering effects are hard to model. Without data, you need to make some guesses.


FFT implementation with decimation in frequency vs decimation in time?

Well, that article equates "FFT" with "Radix-2 FFT", which isn't really true.

Anyway, since both decimation in frequency and decimation in time need exactly the same components running at the same speeds, the only difference I'd expect for radix-2 FFTs would be if you can have a parallel input on one side, and serial output on the other (or vice-versa); in that case, you could do the "large" butterfly on the side where you need to have the data in parallel to avoid latency, and could save on components by re-using a smaller butterfly for the serial end without much of a downside. Whether you want that or not very much depends on space, power, latency and throughput constraints.

For mixed-radix FFTs, where you're not only dealing with a single prime factor in the FFT size, you might not have the freedom to rearrange these layers.


Why do we do oversampling in raised cosine?

A (baseband) pulse amplitude modulated (PAM) signal is given by

$$s(t)=\sum_{k}A_kp(t-kT)\tag{1}$$

where $A_k$ are the data symbols, $p(t)$ is the transmit pulse, and $1/T$ is the symbol rate. Note that $(1)$ describes a continuous-time signal. If you want to simulate a digital communication system, you need to use sampled versions of continuous-time functions, such as the transmit pulse $p(t)$. For this purpose, you need to define a sampling rate greater than the symbol rate. The ratio between the sampling frequency and the symbol rate is the oversampling factor.

Assume you choose a sampling rate an integer (oversampling) factor $m$ greater than the symbol rate:

$$\frac{1}{T_s}=\frac{m}{T}\tag{2}$$

From $(1)$, the sampled signal $s(nT_s)$ is then

$$\begin{align}s(nT_s)&=\sum_{k}A_kp(nT_s-kT)\\&=\sum_{k}A_kp\left(\left(n-k\frac{T}{T_s}\right)T_s\right)\\&=\sum_{k}A_kp((n-km)T_s)\end{align}\tag{3}$$

Or, in notation for discrete-time signals

$$s[n]=\sum_{k}A_kp[n-km]\tag{4}$$


How Is the Bandwidth of an Image Reduced by an Octave in Gaussian Pyramid Procedure?

The model is simple.
Assume we have data I0 on the frequencies range of $\left[ 0, 1000 \right] $ [Hz].

We have 2 Black Boxes:

BlackBox 1 - LPF which always has a cut off at the middle of the spectrum.
BlackBox 2 - Decimation by 2.

What happens if we apply BlackBox 1 on the data?
We'll get an LPF version of it we'll call it L0 with data on the range $ \left[0, 500 \right] $ [Hz]..
We can also have H0 = I0 - L0 which is vasically all data on the range $ \left[500, 1000 \right] $ [Hz].

Now, be decimating I1 = decimate(L0) we have new data on the range $ \left[ 0, 500 \right] $ [Hz].
Applying the same tricks as above (Using the same black boxes which is what's nice about it) can yields the following series:

L# - Always contains the lower half frequency of the # Level.
H# - Always contains the higher half frequency of the # Level.

If you see the amount of data is reduced by factor of 2 each level (Octave).
Another nice property that you can reproduce the data using the lowest level LPF and the HPF of all other levels.


@Yatekii The linked you shared helped solve the problem. I needed to take into account (1) the gain of the window used and (2) summate over a short range of frequencies. Below is the updated code:

% Create signal
    fs = 48000;t = 0:1/fs:5-1/fs;
    A = 0.5;
    x = A*cos((2*pi*1000)*t);
    % Run FFT
    [pxx, f] = pwelch(x, 1024, 512, 1024, fs, 'one-sided');

    % Spectral parameters
    df = f(2) - f(1);
    CG = 0.54;
    NG = 0.3974;

    % Extract peak value
    isig = round(1000/df)+1;
    Aout = sqrt(sum(pxx(isig-5:isig+5)*df)) * sqrt(2);


What is the best way to represent audio visually? (x-post from UX)

What a human (or their ear-brain) perceives in sound is a psychoacoustic phenomena, and may or may not be exactly related to the actual audio as recorded.  e.g. the exact notes, beats and instruments that a human "hears" may be influenced by visual cues, memory of other similar music, and the musical context around the actual sound of the note in question.


It really depends, how will you measure the quality of the solutions? what are your requirements, real time, high accuracy? how large is the point cloud?

You mentioned valid, yet fancy signal processing based methods to address the problem.

Let me add three methods you have not mentioned that are classic and more statistical in nature: least squares, ridge regression and lasso.


How to compute the main lobe width of generalized Hamming windows (i.e. Hann, Hamming, etc.)?

(This used to be a partial answer, but I have edited it to include the requested result at the bottom.)

First of all, you can view a Hamming family window as the sum of a Hann window and a rectangular window: If $W_N$ is the fourier transform of a Hann window, then $W_H(\omega)=(\alpha-2\beta)W_R(\omega)+4\beta W_N(\omega)$. If, for simplicity, we assume unit amplitude ($\beta=\frac{1-\alpha}{2}$), then $x:=4\beta=2-2\alpha$ determines the weighting of the Hann window against the rectangular window ($\alpha-2\beta=1-x$ in that case). In the general case, you just need to factor out the coefficient $\alpha+2\beta$ first and let $x:=\frac{4\beta}{\alpha+2\beta}$.

For example, here is a plot of $W_H$ (green) and its two components for $x=0.5$ (assuming $M=2$, I guess, so that $W_R(\omega)=\frac{\sin \omega}{\omega}$):

Obviously, giving the Hann window more weight by increasing $x$ also increases the main lobe width. However, plotting this relationship reveals an interesting effect:

Since $W_R(2\Omega_M)=W_N(2\Omega_M)=0$, as soon as the Hann window component completely cancels out the first side lobe of the rectangular window component, the main lobe width stays constant at $4\Omega_M$. By calculating the derivatives at this point, it is easy to find out that this is the case for $x\geq\frac{6}{7}$ (equivalently, $\alpha\leq\frac{4}{7}$ in the unit amplitude case, or $3\alpha\leq 8\beta$ in the general case).

Unfortunately, from the plot, I haven't successfully guessed anything about the remaining $\frac{6}{7}$. ;-) Hopefully someone else will figure this out, if it is possible at all.

Edit: Don't know why I didn't go for the simplest route first. Maxima (open-source CAS) gives $\Omega_M \sqrt{2\frac{x-2}{x-1}}$ for the remaining part. So, in general, the main lobe width is:

$\begin{cases}
2\Omega_M \sqrt{\frac{\alpha}{\alpha-2\beta}} & \text{if } 8\beta<3\alpha\\
4\Omega_M & \text{otherwise}
\end{cases}$


The peak in frequency spectrum from sampled data is higher than the true amplitude

10 Hz aliasing down to two Hz which is the Nyquist frequency.
When you sample a Nyquist frequency signal, you only get two samples per period and the two samples are exactly 180 degree apart, i.e. you always get the same numbers with alternating sines.
Which number you get depends on the relative phase fo your sample clock and the signal.If you sample a consine you samle at the peaks and you get [1 -1 1 -1 ...] .  If you sample a sine wave , you sample at the zeros and you get [0, 0, 0, 0 ...]. For any other phase, you get something in between.
Keep in mind that the sampling theorem requires the signal the be bandlimit to "less than half the sample rate" not to "equal or less than half the sample rate".


What's a good physical description of the following matlab code?

The input is audio and the code first performs an STFT/spectrogram (sliding window FFT of size 1024 with 50% overlap, hanning window applied). The output is then filtered to the range (using Matlab notation) to the range (2:N/2). In other words, the first bin is dropped and the negative half frequencies too.
. . .
The value filt_y is then summed and that becomes the output signal:

What you are describing thus far is equivalent to a boxcar bandpass filter followed by integration similar to the way a moving average filter works.

So, for every window's FFT we get a single output value. The output frequency is audio_frequency/512 because we slide the 1024 sample window along by 512 samples at each step.

It would not be straightforward to define "output frequency" here. Most of the times, when you refer to "frequency" you might implicitly be referring to sinusoids. However, here, there is no sinusoidal output. But, because of properties of the DFT, your output will indeed fluctuate at a frequency that depends on the length of the FFT window.
This is because, the FFT assumes that the signal at its input is periodic. In other words, it assumes that whatever lands in its input frame is played again and again in a loop to infinity on either sides (both positive and negative). But this is not what is happening in reality. Because of this, your "running sums" from frame to frame have discontinuities. If you were to listen to that signal, IT!WOU!'D.CO!NTAI!NDIS!CON!TIN!UIT!IES!.
There are however two elements in your code that reduce the effect of those discontinuities. The first, is the use of the hann window and the second is the fact that the summation drops the DC component (and low frequencies depending on NFFT). With the use of the window, the discontinuities between the end of the frame and its beginning are reduced (remember, that is the start and end of the "loop"). By dropping the DC, you drop the result of partial summation that could end up having very high values.
Imagine for example that you are summing a zero mean sinusoid with a "frame" whose length just so happens to be half its period. Then in that case, during the positive phase of the sinusoid you would get a high positive value and then during the negative phase you would get a negative value. If the length of the "frame" was one period of the sinusoid, you would get a zero mean.
If you left the DC in and tried to "listen" to the signal you would be getting a soft "flutter" sound in the background, like holding a microphone close to a fan because of these discontinuities that are produced at each frame.

My question is: what would you call that output signal?

Some sort of band-limited integration. As mentioned previously, this is an integrator with a crude band-pass filter at its input.

What sort of physical significance does it have, or what use is it for analysis?

It returns the total amplitude of the frame of the signal minus its average value. It would be related to the average amplitude in a "moving average filter" sense.
It is difficult to say what its use would be for analysis without knowing what else is going on in the processing pipeline, but something like this could be used to get a sense of how "loud" is the current frame.
If this was a compressor for example, this signal could be used to trigger it at a sharp transition from soft sound to louder sound. In fact, because of the crude band-pass filter it contains at its input, this would then become a de-esser.
Hope this helps.


Converting cascaded biquads to parallel biquads

You need to do partial fraction expansion.
Most of the time, the result will be 3 parallel sections with the same poles as the original filter but only with a real zero in each parallel section. The first section is just a single real scale factor.
$$H(z) = r0_0 + \frac{r0_1 + r1_1 \cdot z^{-1}}{a0_1 + a1_1 \cdot z^{-1} + a2_1 \cdot z^{-2}}  + \frac{r0_2 + r1_2 \cdot z^{-1}}{a0_2 + a1_2 \cdot z^{-1} + a2_2 \cdot z^{-2}} $$
Note that the "a" coefficients are the same as in the cascaded form, i.e. you only need to calculate the "r" coefficients.
This is all relatively straight forward if your poles are unique and complex conjugated pairs. Repeated poles and single real poles require special treatment.
The Matlab function residuez() can typically help with this although it can have numerical problems if the poles are too close to the unit circle. If you need to do it by hand take a look at https://ccrma.stanford.edu/~jos/filters/Partial_Fraction_Expansion.html
Since you only have two sections, the following approach should work

Write both versions (parallel and serial) as a single fraction
Discard the denominators, they are the same
Both numerators are a 4th order polynomial in $z^{-1}$, i.e. they have 5 coefficients.
Each power of $z^{-1}$ needs to have the same coefficient in both numerators.
That gives you five linear equations with 5 unknowns that you can solve for $r0_0, r0_1, r1_1, r0_2, r1_2$


Which physical device accomplish the task of pulse shaping?

It's a combination of two processes: first, you create a discrete-time signal that is pulse-shaped. Then, the DAC converts it to an anlog signal. (And, in a radio, the analog front-end following the DAC converts it to an RF signal).

The discrete-time part is usually done as follows (using Matlab here for convenience). Assume a BPSK data vector data. First, the vector is upsampled to get the desired symbol rate. For example, if the sampling frequency is $f_s = 100$, and the symbol rate is $R_s = 10$, then you must output one symbol every $100/10 = 10$ samples, which means you need to insert 9 zeros between samples:

data_up = upsample(data, 10);

Then, assuming the pulse shape is stored in vector p, the pulse-shaped signal is

data_ps = conv(data_up, p);

Then, data_ps is the input to the DAC, which converts it into a continuous-time pulse-shaped baseband BPSK signal.


Please see my answer to this question: Detecting and isolating part of an image

You can use circular shortest path to do this quite elegantly. You are using OpenCV and and not MATLAB so some modification of my code is required.


I think the system is always non-causal if it has more zeros than poles. @Hilmar, your implementation might be possible digitally but if you see a transfer function like this:
$H(z) = \frac{(z-1)(z-2)}{(z-3)}$, taking the inverse Z-Transform,
$y(n-1)-3y(n-2) = x(n)-3x(n-1)+2x(n-2)$ is clearly non-causal and cannot be implemented real-time(of course you can store the signals and do processing offline, but real time implementation is impossible). So I guess Non-Causality is the answer, apart from the infinite high-frequency gain.

Update by Val
Ok, I understand now. In wikipedia, we see that $$a_0 y_n + a_1 y_{n-1} + a_2 y_{n-2} + \cdots = b_0 + b_1 x_{n-1} + b_2 x_{n-2} + \cdots$$ is identical to $$H(z) = {b_0 + b_1 z^{-1} + b_2 z^{-2} + \cdots \over a_0 z^{0} + a_1 z^{-1} +\cdots}.$$ I have finally noted that there are no positive powers. That is, $a_0$ corresponds to the current output, $z^0 = y[n]$, and $\cdots = a_{-2} = a_{-1} = \mathbf{0} = b_{-1} = b_{-2} = \cdots$. If we had more zeroes then we have more power in the nominator and, after normalizing $H(z)$ so that $a_0$ has $z^0$, we still have a positive power of z in the nominator, which means that output depends on future input.

In the simplest case, $Y(z) = zX(z)-x_0$ corresponds to $y[i] = x[i+1]$, that is output depends on future input (we can also drop $x_0$ assuming our sequences are causal, i.e. $0 = y_{-1} = x_0$ aka $z^{-1}Y = X$ or $Y = zX$). It is curious that causal is an example of non-causal. Could anybody clarify it further?

It also worth noting that $(1-az^{-1})(1-bz^{-1})$ is not a zero of order two. It has equal number of zeroes and poles, since it is equal to $(z-a)(z-b)\over z^2.$


One way to solve such a difference equation is to use the $\mathcal{Z}$-transform. Transforming the given difference equation yields

$$X(z)\left(1-z^{-1}\right)=Y(z)\left(-10z^{-1}+2+2z\right)\tag{1}$$

which results in the following transfer function:

$$H(z)=\frac{Y(z)}{X(z)}=\frac{1-z^{-1}}{2z+2-10z^{-1}}=\frac12\frac{z^{-1}-z^{-2}}{1+z^{-1}-5z^{-2}}$$

The poles of $H(z)$, i.e. the zeros of the denominator polynomial in $(1)$, are both real-valued:

$$1+z^{-1}-5z^{-2}=(1-az^{-1})(1-bz^{-1})\tag{2}$$

with $a=\frac{-1+\sqrt{21}}{2}\approx 1.79$ and $b=\frac{-1-\sqrt{21}}{2}\approx -2.79$. Note that $1<|a|<|b|$ is satisfied, so both real-valued poles are outside the unit circle of the complex plane.

With $(2)$, $H(z)$ can be rewritten as

$$H(z)=\frac12\frac{z^{-1}-z^{-2}}{(1-az^{-1})(1-bz^{-1})}=\frac{1}{2(b-a)}\left[\frac{(1-a)z^{-1}}{1-az^{-1}}-\frac{(1-b)z^{-1}}{1-bz^{-1}}\right]\tag{3}$$

This transfer function can have three possible regions of convergence (ROCs), each of them corresponding to a different impulse response:

$$\begin{align}(1)\quad &|z|>|b|>|a|\\
(2)\quad & |b|>|z|>|a|\\
(3)\quad & |b|>|a|>|z|\end{align}$$

ROC $(1)$ corresponds to a right-sided impulse response, ROC $(2)$ corresponds to a two-sided impulse response, and ROC $(3)$ corresponds to a left-sided sequence. Only ROC $(3)$ includes the unit circle $|z|=1$, so only the impulse response corresponding to this ROC results in a stable, yet anti-causal, system. The three impulse responses corresponding to the three ROCs given above are

$$\begin{align}(1)\quad & h_1[n]=\frac{1}{2(b-a)}\left((1-a)a^{n-1}-(1-b)b^{n-1}\right)u[n-1]\\
(2)\quad & h_2[n]=\frac{1}{2(b-a)}\left((1-a)a^{n-1}u[n-1]+(1-b)b^{n-1}u[-n]\right)\\
(3)\quad & h_3[n]=\frac{1}{2(b-a)}\left(-(1-a)a^{n-1}+(1-b)b^{n-1}\right)u[-n]\end{align}$$

These impulse responses describe three different systems, all of them characterized by the same transfer function given by Eq. $(3)$. $h_1[n]$ corresponds to a causal unstable system, $h_2[n]$ to a non-causal unstable system, and $h_3[n]$ to an anti-causal stable system.


What is meant by 20 dB signal-to-noise ratio?

A discrete signal is often interpreted as an analog voltage signal. The signal x(n) should be considered as a sampled version of x(t) with a unknown sampling frequency.

Q1:

So when we're talking about the energy of the signal, it is presumed that the voltage signal is obtained in a circuit with an 1Ω resistor in series, so that v(t)^2/R = P.

Q2:

Signal-to-noise ratio is used in many areas and is defined as the logarithmic power ratio of the signal and noise, i.e. SNR = 10*log(Ps/Pw). The other SNR = 20*log(Ps/Pw) is not true, because the Ps and Pw are power values. The latter is used when the signals are represented as "voltage" signals, i.e. SNR = 20*log(s(n)/w(n)) = 10*log(s(n)^2/w(n)^2).


Designing a feature vector for discriminating between different sonic waveforms

You want only objective criteria to seperate the signals or is it important that they have some sort of similarity when listened to by someone? That of course would have to restrict you to signals a bit longer (more than 1000 samples).


How does the intuitive notion of causality fit in with control systems?

As near as we can tell by experiment, causality is nature's way of doing its thing.
Causality says that if you have a system $y(t) = h\left(x(t), t\right)$, and it is causal, then $y(t_0)$ is dependent only on values of $x(t)$ for $t < t_0$*.
Causality doesn't say you can't combine systems into larger systems (as a rather pertinent example, you can certainly combine systems in a way that makes a control loop).
Causality doesn't imply stability, or instability.  There's a rich set of behaviors that causality allows -- it just doesn't allow the output to follow future inputs.
So in your exponential growth example, if you have a system described by $\dot x = a x + u$, and $a > 0$, the output will tend to grow exponentially (unless $u$ magically or through feedback damps the tendency).  We can say that the growth is caused by $x$ feeding back on itself -- but that causation still happens without violating causality.  Any change in $x$ or $u$ at time $t_0$ will affect the derivative of $x$ at time $t_0$, but (in the absence of infinities) it won't change the immediate value of $x$ -- so, the system is still causal.
* I may be overly strict here -- I'm ruling out $x = u$ as a causal system.  I'm not sure what the popular definition is in control system terms, but I'm lazy so I'm going to let it stand.


How do I get the most accurate camera calibration?

You have to take images for calibration from different points of view and angles, with as big difference between angles as possible (all three Euler angles should vary), but so that pattern diameter was still fitting to camera field of view. The more views are you using the better calibration will be. That is needed because during the calibration you detect focal length and  distortion parameters, so to get them by least square method different angles are needed. If you arn't moving camera at all you are not getting new information and calibration is useless.
Be aware, that you usually need only focal length, distortion parameters are usually negligible even for consumer cameras, web cameras and cell phone cameras. If you already know focal length from the camera specification you may not even need calibration.

Distortion coefficient are more present in "special" cameras like wide-angle or 360°.

Here is the Wikipedia entry about calibration.
And here is non-linear distortion, which is negligible for most cameras.


How do I interpret the result of a Fourier Transform?

The Fourier transform of a function $x(t)$ is defined by

$$X(\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\tag{1}$$

Since $e^{-j\omega t}$ is complex-valued, the Fourier transform $X(\omega)$ is generally a complex-valued function, even if $x(t)$ is real-valued.

$\omega$ is the frequency in radians per second (rad/s), and it is related to the frequency in Hertz $f$ by $\omega=2\pi f$.

Since $X(\omega)$ is generally complex-valued, it can be represented by its real and imaginary parts

$$X(\omega)=X_R(\omega)+jX_I(\omega)\tag{2}$$

or by its magnitude and phase

$$X(\omega)=|X(\omega)|e^{j\phi(\omega)}\tag{3}$$

Clearly, the real and imaginary parts are not the same as magnitude and phase.

Note that for real-valued $x(t)$, the real part $X_R(\omega)$ and the imaginary part $X_I(\omega)$ of the Fourier transform $X(\omega)$ are even and odd functions, respectively. Consequently, the negative frequency components are redundant, and knowledge of the spectrum at positive frequencies is sufficient. If you are interested in power, you should rather plot $|X(\omega)|^2=X_R^2(\omega)+X_I^2(\omega)$ (which, by definition, is non-negative, as it should be the case for energy or power). The function $|X(\omega)|^2$ shows the distribution of the frequency components of $x(t)$.

As you've noted, the DC component of the signal is zero, because its average value over time is zero:

$$X(0)=\int_{-\infty}^{\infty}x(t)dt=0\tag{4}$$


$$1 - e^{-4j\omega} = e^{-2j\omega}(e^{2j\omega}  - e^{-2j\omega}) \tag{1}$$

Now,

$$ \sin(2\omega) = \frac{e^{2j\omega}  - e^{-2j\omega}}{2j} \tag{2}$$

Equation 2 is a consequence of Euler's formula. Multiply and divide by $2j$ in (1) and use the identity (2) in equation 1 we have:

$$1 - e^{-4j\omega} = 2je^{-2j\omega}\sin(2\omega) \tag{3}$$

Now $j = e^{\frac{j\pi}{2}}$, substituting in 3 we have

$$1 - e^{-4j\omega} = 2e^{-2j\omega}e^{\frac{j\pi}{2}}\sin(2\omega) \tag{4}$$

The term $e^{-2j\omega}$ in (1) is represntative of the delay of the envelope in time domain, the term $(e^{2j\omega}  - e^{-2j\omega}) $ is reprsentative of the magnitude response and actual shape of the filter in frequency domain.

A general trick in simplifying expressions of FIR responses type is to tap expressions around the mid point of the impulse response taps, you could use this strategy to pair terms in responses with more than 2 terms as well


How to generate colored Gaussian noise and adding it to a ODE system - Do I need Euler-Maruyama method?

Question 1) How to generate colored noise in general?

Produce white noise and filter it.

  Question 2) How to add the colored noise in the form of process noise to the Rossler system or any continuos time system

This is a little bit more complicated. The addition of noise is a simple matter of $y = x + q$ where $q$ is the noise signal. The key problem here is how much noise to add? This is important for linear cases (such as $y$ previously) because random number generators are producing numbers in some fixed interval ($[0 \ldots 1]$, $[-1 \ldots 1]$ or other) which will be entirely different than the interval the signal $x$ occupies.

In this case, assuming that the signal is stationary, then given some $x$, generate a $q$, scale $q$ so that it has the same power as $x$ and then mix them together at the desired ratio. For example $y = (1 - a) \cdot x + a \cdot q$. Here, the contribution of noise is directly proportional to $a$.

This is nice because you can now use $a$ later on in your experiments as a variable that shows how much noise was added to the signal.

But, when you are dealing with non-linear systems, this becomes even more important (and difficult). The reason for this is that, depending on the system, you cannot assume that it is stationary. Or, you can, but only at a narrow range of operation and furthermore, you are going to have to make sure that the system stays in this region.

The added difficulty is in two places:

1) You are dealing with differential equations that get integrated

Therefore small spurious excitations here and there can accumulate to large numbers quickly

2) You are dealing with non-linear systems

Therefore, combined with #1 it becomes likely that you might push the system into a region that its behaviour changes. For example, you might be interested to use the system when it produces nice smooth waveforms but inadvertently you might tip it over to the region where it starts producing spiking oscillations.

The reason I am mentioning this is because you refer to "process noise". So, you can take the Rossler model and integrate it and look at its signature attractor as an experiment or, you might want to use it as an actual model of a process, in which case, $x,y,z$ now have units, $dt$ is tied to your sampling frequency or other physical limit and also your noise now has a physical meaning too.

I don't know what these are in your case, but you have to understand them in order to be able to add the noise, especially if you have to add the noise in the derivatives (are you sure about this?).

The derivative is not the signal but rather the rate of change of the signal with respect to time. This is incredibly powerful. If you were to produce some $x$ by integrating the $\frac{dx}{dt} = b$, you would get a ramp in $x$. If your derivative now was a ramp itself with $\frac{dx}{dt} = c \cdot t + b$, your $x$ would vary exponentially.

Imagine then, how much you are impacting $x$ by specifying some random factor in its $\frac{dx}{dt}$.

If you don't have to add the noise at the derivatives, it would be slightly easier, because you simply integrate your differential equations to produce some output $z$ and then add noise to this specific $z$ as if it was our $x$ from above.

  ...and which method to use for solving the system -- ODE or Euler Maruyama? Is the noise added Rossler system stochastic or not?

If you have to add the noise to the derivatives then the system becomes stochastic and it would be better to use the Euler - Maruyama. The difference of Euler - Maruyama from the plain Euler integration is in the inclusion of a derivative term for the random process as well. Of course, if you differntiate white noise, you still get white noise, so from a practical point of view, there is no difference there because you still just have to add a random number. BUT this random number has to be appropriately scaled for its amplitude and then be multiplied by the integrating step too.

Hope this helps.


Typically both are used, coming from an 802.11p perspective. The short preamble (STF) is used for packet detection and coarse frequency offset. The long training preamble (LTF) is used for timing synchronization and fine frequency offset. The LTF can also be used for channel equalization. Both are used from my understanding to accomplish total synchronization.


Doppler Model, how does it affect the fading impulse response?

This is not a full answer, but I'll try to provide some intuition. First, imagine a static channel -- neither the antennas nor the reflectors around them move. In such a channel, the impulse response is always the same -- it's time invariant. The channel characteristics are fixed.

Now, introduce movement of the antennas and/or the reflectors. The channel impulse response now varies with time. The channel coherence time indicates how fast the channel impulse response changes, and for how long can we consider the channel to remain relatively static.

If there is movement of the antennas and/or the reflectors, we also get a Doppler shift in the transmitted signal. The more movement there is, the higher the shift.

The channel coherence time and the Doppler shift are, then, two manifestations of the same phenomenon -- movement. If one increases, the other decreases. The easiest way to introduce this into the channel model is to assume that the channel is fixed for a certain amount of time, after which it changes to a completely different, independent channel.


In the context of transfer functions, what is the relationship between the terms "proper", "causal", and "realizable"?

Causality is a necessary condition for realizability. Stability (or, at least, marginal stability) is also important for a system to be useful in practice.

For linear time-invariant (LTI) systems, which are fully characterized by their transfer function, we get realizability constraints on the transfer function. For continuous-time LTI systems, if we work at frequencies for which the lumped element model is valid, we require the system's transfer function to be rational for the system to be realizable. Also for discrete-time LTI systems we require rationality of the transfer function, which implies that the system can be realized by adders, multipliers, and delay elements.

For an LTI system to be causal and stable, its poles must lie in the left half-plane (for continuous-time systems), or inside the unit circle (discrete-time systems). From this it follows that the rational transfer function of an LTI system must be proper, otherwise you would get one or more poles at infinity, causing the system to be unstable (or non-causal).


Determine the stability of a system without using the $\mathcal Z$-transform (described by a difference equation)

It boils down to the same thing, with or without $\mathcal{Z}$-transform. You determine the roots of the characteristic equation, which in your case is

$$\lambda=a\tag{1}$$

For stability, the magnitude of the roots of the characteristic equation must be smaller than $1$. So for your example you get the condition $|a|<1$ for stability.


Learner level information on Kalman filtering for different input kinds

I am most convinced you should try Hidden Markov Chains (which i am not explaining here), just justifying its use:

Kalman Filter: continuous state space, discrete observations

https://stats.stackexchange.com/questions/183118/difference-between-hidden-markov-models-and-particle-filter-and-kalman-filter


Energy per bit ($E_b$) in BPSK with non-unity channel gain

The energy per bit $E_b$ means one thing and one thing only: how many joules is the transmitter spending, on average, per information bit transmitted. Note that the channel has absolutely nothing to do with $E_b$.

However, your question involves the SNR too. In a channel with non-unit gain, you need to decide where you'll measure the SNR.

You can define the SNR at the transmitter: $\text{SNR}=2E_b/N_0$ for BPSK. Obviously, the BER depends on both the SNR and the channel gain.

You can also define the SNR at the matched filter's output: $\text{SNR}=2G^2E_b/N_0$ for BPSK, for a channel with gain $G$. In this case, the BER depends exclusively on the SNR.

When $G$ is random with known probability density, it is common to average the SNR over $G$, in order to obtain the expected BER over many channel realizations.


Unexpected results for fourier resampling (Python 3)

If you perform the full DFT of a signal with odd number of samples, the frequency values in the spectrum are a DC coefficient and (N-1)/2 conjugate pairs of sinusoid coefficients.

If you perform the full DFT of a signal with even number of samples, the frequency values in the spectrum are a DC coefficient, (N-1)/2 conjugate pairs of sinusoid coefficients, and a single coefficient corresponding to frequency 0.5.

rfft only keeps the positive frequency part of the spectrum, and the inverse, irfft, presumably recreates the negative part of the spectrum before performing an inverse-DFT.

For example, for an odd size signal with rfft of [0 1 2 3], the full spectrum would be [3 2 1 0 1 2 3].

For an even size signal with rfft of [0 1 2 3 4] it presumably has a flag set somewhere in the code that indicates an even sized signal, and so the full spectrum would be [3 2 1 0 1 2 3 4]. Note the 4 only appears once.

When you pad the rfft result though, e.g. [0 1 2 3 4 0 0 0 0], it probably recreates the spectrum as [0 0 0 4 3 2 1 0 1 2 3 4 0 0 0 0]. That is, the 4 is now repeated twice, and therefore that frequency component is twice as strong.

Therefore, before padding you must halve the last rfft coefficient iff the signal has an even number of samples.


Detecting a specific pop in a real time audio signal

One of the simplest methods to optimally (minimising the error energy while restricted to LTI filters) detect the existance of a known signal inside another is to use a matched filter. A matched filter will have its impulse response equal to the signal being searched in time-reversed form. This resulting filter equivalently computes the correlation of the known signal with the measured signal. And the peaks of this correlator yields the position of the best possible matches.

Let the searched signal be $s[n]$ with $s= \{-2,-1,0,0,1,2\}$ with $s[0]=-2$ and $s[5]=2$ and length of signal is $L=6$.

Assume that this signal is embedded into various locations of another longer and uncorrelated signal x[n] with its length M > N. Then the resulting combined (by adding sifted versions of the test signal s[n] into x[n] at various locations) signal is called as y[n] of length M.

Now to optimally find the index of the locations into where the signal s[n] was added, we can use an LTI filter with impulse respose h[n] which is the time reversed version of s[n] as: $$h[n] = s[-n]$$

(you can refer to Simon Haykin Communication Systems Theory for a simple and accessible derivation of this optimal filter which is called as the matched filter. Note that the impulse response h[n] of this optimal filter is the same (except time reversal) as the signal being searched, therefore it's called as "matched": optimal detector's impulse response matches the signal being searched)

Following the example of $s[n] = \{-2,-1,0,0,1,2\}$ then $h[n] = \{2,1,0,0,-1,-2\}$
Now however $h[-5]=2$ and $h[0]=-2$. This is an non-causal signal. In practice (using computer system) it is wiser to shift it right, enough to make it causal. It is just a convenience and will not effect the operation performance. Then we will have the following same impulse response $h[n]= \{2,1,0,0,-1,-2\}$ with $h[0]=2$ and $h[5]=-2$ for mathematical convenience.

Now the output of the LTI filter is $$z[n] = h[n]*y[n]$$ This is the convolution operator, which enables the computation of the output of the any discrete time LTI filter and is expressed as $$z[n] = \sum_{k=0}^{k=M-1} {h[n-k]x[k]}$$ The dummy index k ranges within the intersection of the valid signal ranges which is k=0 to k=M-1 for this particular example. Also the output index n will also range from 0 to M-1, which is by convention...

So the output sequence z[n] is expected to inlcude some peaks, these peak locations after being thresholded indicates a possible detection of signal s[n] inside y[n]. (practically the exact location will be L samples back)

What if the signal s[n] is not only shifted but also scaled? Then the peak value will be reduced as much as the scaling. Hence the threshold must be selected very carefully. Or you must find a way to estimate the scaling a priorily.


Why does slowing down audio introduce aliasing?

That's because you are interpolating the signal. Take a look at this figure:

This is zero interpolation, i.e. inserting zeros between samples, but applies to your case too. The replicas at 2*pi are distorting your signal. You can see them in the third graph, at -2pi/L, +2pi/L. So, you need a low pass filter with a cut off frequency pi/L, after expander.

In your case, the expander's output dtft will be X(w*2) + exp(-jw)*X(w*2), assuming input signal's dtft X(w). So, your extrapolator output will look a bit different, but the idea is the same.


Mathematics / Signal theory behind billiard ball 'wave pendulum' effect

As was pointed out in a comment above, this is a simple consequence of the dynamics of a pendulum. There's nothing particularly signal-processing-related about this problem, just some simple physics and trigonometry.

For a pendulum that is displaced from its angular equilibrium point by an angle of $\theta_0$ with an angular velocity of $\omega_0 = 0$ at $t = 0$, for small initial angles, the motion of the pendulum can be approximated by simple harmonic motion:

$$
\theta(t) = \theta_0 \cos\left(\sqrt{\frac{g}{l}}t\right)
$$

where $g$ is the magnitude of the force of gravity and $l$ is the length of the pendulum. In the video, a number of pendulums of varying lengths are given an approximately equal initial displacement $\theta_0$ and released simultaneously. Once they are released, each pendulum traces out a path whose displacement angle is a sinusoid with frequency related to the pendulum's length. The lengths were carefully chosen to give the interesting visual effect that you see in the video; the frequencies of the sinusoids are multiples of some fundamental frequency, so for particular values of $t$, they line up in interesting ways.

You'll notice toward the end, all of the balls "line up" again; this is because each of the sinusoids that describe their motion have moved through an integer number of periods, so the balls are back at a common angular displacement. If there were no effects of friction, air resistance, etc., the oscillation would continue, with the same patterns repeating indefinitely.


How many bytes would I require to store few seconds of fan controller PWM signal?

The easiest answer is: find the maximum frequency in the signal, and double that. In the case of PWM, call $T_p$ the duration of the shortest pulse possible. The signal bandwidth is in theory infinite, but in practice, frequencies above $10T_p$ will contain very little energy. So, you can safely sample the signal at $20T_p$ samples per second.

If you need better fidelity, then sample faster, accordingly.

Now, a PWM signal typically has only two amplitudes, so you can use one bit per sample, or $20T_p$ bits per second.


what sampling rate and bit depth should i use?

The typical audio sampling rate is 44 kHz, so you would need between 130 000 and 175 000 samples.
In practice though, it depends on the sampling rate you use on your FPGA. You can lower, but you will need to low-pass your signal before sampling it to avoid aliasing.

Note that if you put a threshold on z[n]=x[n]*y[n], it will not be invariant to volume differences.


Minimum Number of Errors required for BER Curves in Matlab Simulation

I have done a fair number of BER simulations and have not found a definitive answer on this. It depends on the error behavior of your channel, Gaussian or some other pdf. My process for AWGN was as follows. First I determine the approximate BER I am expecting. Say that it is 10-3. The inverse of that is 1000 bits. So I would run at least 10x that to start, or a length of 10,000 bits as a minimum. For coded signals, or those with fading, I might use a larger multiplier, such as 50 to 100 times. By watching the error pattern we can often tell if steady state has been achieved. For coded channels with BERs of 10-8 etc., the run lengths quickly become an issue.

This applies to Monte Carlo simulations. You can do semi-analytic methods that do not require such large runs. However semi-analytic methods are generally not that good for non-AWGN channels. Importance sampling is one such semi-analytic method which is used with coded signals. All of these are estimates.


No transform is required for that kind of thing except making sure that the rays of light hit the lens head on.

This is only possible when the subject is far away enough from the camera and for a narrow range of its field of view.

In other words, the end result image we are talking about here is composed of vertical projections of 3D points right on the back of the camera. As if you had a one pixel lens which you scanned across the full apperture of your scene to sample exactly what is right at the focus point of that lens. Head on.

But what really happens is that for every picture you take, the light goes through a single lens. The only pixel that is vertically sampled (or projected, from the world to the sensor) in this case is the one at the centre of that lens. Everything else is imaged at an angle that is proportional to how far away from the center of the lens the sampling pixel is located.

Take the picture of a very tall multistorey building for example. At eye level, straight ahead, the windows appear rectangular. But as you raise your gaze upwards, you "sample" the windows higher up at an increasing angle to the point that some of the windows might even start to appear square (their perceived height, approximately equal to their perceived width).

If you do not have a 3D mesh of the scene, there is no transform that will correct for this change in perspective.

So, what you do is to "fake it".

You crop a patch from the acquired images that given the size of the subject, its distance from the camera and the focal length of the lens, it appears to be imaging the subject head on. And you take lots of images.

No further reprojection is required except of course to correct for the relative camera poses between takes as the camera is not on a leveled rail but that's scaling and rotation only.

Depending on the application this is sometimes called the Pushbroom camera model.

Provided that you have taken care during acquisition, Hugin can do that type of stitching.

Hope this helps.


Output gain when upsampling and downsampling

All of the answers above are good in their own way. The reason that the Oppenheim-Schafer book and other similar resources use a gain of $L$ for upsampling and a gain of $1$ for downsampling is to maintain a constant amplitude for the signal before and after the operation.

  Do the processes of upsampling and downsampling affect the magnitude of the transform of a signal?

The short answer is that it depends on whether we are upsampling or downsampling. Upsampling does alter the magnitude of the transform (more on this in a minute). Downsampling doesn't alter the magnitude of the transform (assuming a perfect anti-aliasing filter).

For me, this is easiest to understand with a concrete example. So, let's consider a simple sinusoid, $x(t) = A \cos ( 2 \pi f_0 t + \phi )$. Let's sample this at some sampling rate $f_s$ higher than the Nyquist rate (i.e. $f_s > 2 f_0$). This results in a sample sequence

$$ x[n] = A \cos \left( \frac{2 \pi f_0}{f_s} n + \phi \right). $$

The process of upsampling consists of inserting zeros in between every other sample. A mathematical way to express this is

$$ x_{\uparrow 2}[n] = \left\{\begin{array}{cr} A \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right) & \mbox{$n$ is even} \\ 0 & \mbox{$n$ is odd} \end{array}\right. $$

and another equivalent way is

$$ x_{\uparrow 2}[n] = \frac{A}{2} \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right) + (-1)^n \frac{A}{2} \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right). $$.

If we note that $\cos(\pi n) = (-1)^n$, then we obtain

$$ x_{\uparrow 2}[n] = \frac{A}{2} \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right) + \frac{A}{2} \cos\left( \pi n \right) \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right). $$

From here, we apply the identity $\cos(\alpha) \cos(\beta) = \frac{1}{2} \left[ \cos(\alpha + \beta) + \cos(\alpha - \beta) \right]$ and use the fact that $\cos(\pi - \alpha) = \cos(\pi + \alpha)$. This yields the final form that we care about

$$ x_{\uparrow 2}[n] = \frac{A}{2} \cos \left( \frac{2 \pi f_0}{2 f_s} n + \phi \right) + \frac{A}{2} \cos \left( \frac{2 \pi \left(f_s - f_0\right)}{2 f_s} n + \phi \right). $$

This is the signal prior to passing through the interpolation filter. Notice that now we have $2$ sine waves. In general, if we upsample by $L$ then we have $L$ sine waves. Note that each one is weighted by $1/L$, this ensures that the summed amplitude of all sine waves has the same value as the original sine wave at the correct sampling instants and is zero elsewhere. After applying the interpolation filter, consider the transform of the sequence. It is the same obtained by sampling the original function at the higher rate, but it is scaled by $1/L$. Thus, to preserve the signal amplitude while filtering out the aliased images (in this case the additional sine wave component), we must scale the passband of our interpolation filter by a factor of $L$.

For the second part, assume that we are downsampling by $2$. We start with the same $x[n]$, but this time we choose $f_s > 4 f_0$ so that it isn't rejected in the anti-aliasing filter. The decimated signal is now given by

$$ x_{\downarrow 2}[n] = A \cos\left( \frac{4 \pi f_0}{f_s} n + \phi \right).$$

Unlike the upsampling case, this remains a single sinusoid within the Nyquist bandwidth. So, the magnitude spectrum is unscaled compared to the spectrum we would have arrived at by sampling the original signal at a slower rate.


MR Segmentation to draw marker over the MR Dicom Images

Typically research papers handling DICOM standard images would probably extract the ROI using pixel processing methods. However, a lot of medical imaging companies like to implement their own DICOM tags(public/private) to keep track of such information for various purposes.

For example, General Electric keeps the following tag enabled on most MR images:

VerticesOfTheRegion found at (0028,9503), which is according to the DICOM standard, a private tag(odd location tags are private)

In short: Its implausible to tell you,how was the ROI extracted without having access to the DICOM metadata of the image in question.


Discreteness and periodicity in Fourier transform

This is, essentially, what the sampling theorem is about.  Uniform sampling in one domain (e.g. the time domain) causes periodic extension in the reciprocal domain (e.g. frequency domain).
The reason why is that the sampling function is a periodic function which means it can be represented as a Fourier series
$$\begin{align}
\operatorname{III}_T(t) \ &\triangleq\  \sum_{k=-\infty}^{\infty} \delta(t - k T) \\
 &= \sum_{n=-\infty}^{+\infty} c_n e^{j 2 \pi n \frac{t}{T}}  \\
\end{align}$$
where the Fourier coefficients are
$$ \begin{align}
c_n\, & = \frac{1}{T} \int_{t_0}^{t_0 + T} \operatorname{III}_T(t) e^{-j 2 \pi n \frac{t}{T}}\, \mathrm{d}t \quad ( -\infty < t_0 < +\infty ) \\[4pt]
& = \frac{1}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} \operatorname{III}_T(t) e^{-j 2 \pi n \frac{t}{T}}\, \mathrm{d}t \\[4pt]
& = \frac{1}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} \delta(t) e^{-j 2 \pi n \frac{t}{T}}\, \mathrm{d}t \\
& = \frac{1}{T} e^{-j 2 \pi n \frac{0}{T}} \\[4pt]
& = \frac{1}{T} \ .
\end{align} $$
Using the definition of Fourier Transform most common with electrical engineering
$$ X(f) \triangleq \int\limits_{-\infty}^{+\infty} x(t) e^{-j 2 \pi f t} \ \mathrm{d}t $$
$$ x(t) = \int\limits_{-\infty}^{+\infty} X(f) e^{+j 2 \pi f t} \ \mathrm{d}f $$
when sampling $x(t)$
$$\begin{align}
x_\text{s}(t) &\triangleq x(t) \cdot T \cdot \operatorname{III}_T(t) \\
 &= x(t) \cdot T \sum_{k=-\infty}^{\infty} \delta(t - k T) \\
 &= T \sum_{k=-\infty}^{\infty} x(t) \ \delta(t - k T) \\
 &= T \sum_{k=-\infty}^{\infty} x(kT) \ \delta(t - k T) \\
\end{align}$$
which shows how $x(t)$ is converted to samples $x(kT)$ and this is also true:
$$\begin{align}
x_\text{s}(t) &\triangleq x(t) \cdot T \cdot \operatorname{III}_T(t) \\
 &= x(t) \cdot T \cdot \sum_{n=-\infty}^{+\infty} \frac1T e^{j 2 \pi n \frac{t}{T}} \\
 &= \sum_{n=-\infty}^{+\infty} x(t) \ e^{j 2 \pi n \frac{t}{T}} \\
\end{align}$$
the resulting Fourier transform is
$$ X_\text{s}(f) = \sum_{n=-\infty}^{+\infty} X\left(f - n\tfrac{1}{T} \right)
 $$
which is periodic in the frequency domain with period $\frac{1}{T}$.
Because of the Duality of the Fourier Transform, if can be also shown that sampling in the frequency domain causes periodicity in the time domain, which is essentially all that Fourier series is about.


What is $F(0)$ is "dc" component in the context of image processing?

The zero frequency of an image DCT is the mean gray value of the pixels of the input image (for graylevel images).
$F(0)$ is not an image: it is a single coefficient.


Algorithm for implementing an IIR digital filter, Chebyshev type I low pass

okay, this is, or can be, stuff straight outa a textbook. by "SOS", you mean "2nd-order sections"? i usually call those "biquads".  maybe that's not the best term for it in the s-plane.  i dunno.

anyway, defending on your passband ripple, you should have the resonant frequency and Q for each LPF biquad.  you can use those two parameters and get digital biquads directly out of that using the Audio EQ Cookbook.  google will find it fast for you.

otherwise convert your Tchebyshev from s-plane to z-plane however you want.  the Cookbook uses the bilinear transform, which means every bump or feature in the analog frequency response will have a corresponding bump or feature in the digital frequency response.  if you're more interested in matching the impulse response, you would transform H(s) -> H(z) using the Impulse Invariant method.

someone will probably vote me down for not answering the question completely (they've done it in the past).  big deel.


You actually did a two dimensional Fourier transform there. The 0 Hz frequency is the ,,DC offset'' of your signal, that is, its average value. In computers, it's no wonder that there's a high average value, because signals are often defined using positive numbers.

In a 2D Fourier transform, the coefficients of the transformed correspond to products of 1D Fourier basis elements, that is, signals with a period equal to the length of the signal. You could call these overtones.

2D Products of two such signals are not ,,orientation agnostic''. In fact, they have obvious orientation along the edges of the signal. It is not entirely clear what you are asking for here.


How many times does a certain sound appear in an audio clip?

How can I count these hit sounds ( I should get six in this clip).

In this case you can simply rectify the signal (i.e. run it through an absolute function, then through a simple running average low pass filter and at the output of this add a threshold. Essentially, when the power of the input signal exceeds the threshold you will be receiving a pulse whose width is approximately equal to the duration of the "hit". Alternatively, you can do away with rectification and filtering and simply pass the whole thing via a threshold operator. That is, when the power of the input signal exceeds the threshold, you get a "pulse". In either way, count the "pulses" and you can count the number of "hits" that you have in that wave file.

  What is the best tool/software/app for this problem.

This is not an incredibly challenging Digital Signal Processing task, so you could simply write a Python script that does it in 3-5 lines of code (I am not implying that Python is not capable of serious processing here). Alternatively, you could install something like scilab and again do it in 3-5 lines of code, using one of its toolboxes (or, this one).

  Where can I find a good source that explains the theory behind this.

This "rectify and threshold" technique is a sort of a "standard" in counting pulses, alternatively, you might want to have a look at the moving average filter as an example of a sliding window technique where you process a long signal in small "frames" of samples. Other than this, this website has an extensive list of book references if you search for "books".

Hope this helps.


Python: Least Squares Support Vector Machine (LS-SVM)

There is a package called FukuML.

In their description (Version 0.4.1) they write:

Support Vector Machine

Primal Hard Margin Support Vector Machine Binary Classification Learning Algorithm
Dual Hard Margin Support Vector Machine Binary Classification Learning Algorithm
Polynomial Kernel Support Vector Machine Binary Classification Learning Algorithm
Gaussian Kernel Support Vector Machine Binary Classification Learning Algorithm
Soft Polynomial Kernel Support Vector Machine Binary Classification Learning Algorithm
Soft Gaussian Kernel Support Vector Machine Binary Classification Learning Algorithm
Polynomial Kernel Support Vector Machine Multi Classification Learning Algorithm
Gaussian Kernel Support Vector Machine Multi Classification Learning Algorithm
Soft Polynomial Kernel Support Vector Machine Multi Classification Learning Algorithm
Soft Gaussian Kernel Support Vector Machine Multi Classification Learning Algorithm
Probabilistic Support Vector Machine Learning Algorithm
Least Squares Support Vector Machine Binary Classification Learning Algorithm
Least Squares Support Vector Machine Multi Classification Learning Algorithm
Support Vector Regression Learning Algorithm

As you can see, they have the feature you're after.


Filter design, relationship between energy and stop-band ripple

$\delta_s$ is the stop-band (hence the ${}_s$ subscript) attenuation.

A signal in the stop-band with amplitude $1$ would have amplitude $\delta_s$ after filtering.

Since power goes quadratic with amplitude, and if you feed in white noise, $\frac{22}{25}$ of the original energy will be in stop band, and will be reduced by a factor of $\delta_s^2$.


Interesting code! This was tough to follow from the pictures, but I did run the OP's code from the linked site, I can't play the sound loud enough to really test but appears to not be very sensitive for lower volume signals. There is no indication of the frequency axis on the graphic but I assume the left side of the screen is DC and the right side of the screen is half the sampling rate.
The sampling rate in the code was 4200 (I assume Hz) which would be way too low for audio and suggest the OP is seeing the aliasing of all higher frequencies above the Nyquist frequency of 2100 Hz. Since the audio system has a likely bandwidth of 15 KHz, any signals in multiples of 2100 Hz up through that frequency range would fold in and be a more obvious explanation of multiple frequencies due to aliasing.
I did change sampleRate to 44100 and it appeared (in my minimum ability to test) to behave more as expected.
Room echoes on their own should not be non-linear meaning they would not create new frequencies but only modify the phase and amplitude of the frequencies that already exist.


What exactly is captured in a Sample of the Nyquist-Theorem variety?

It might be better to think of Nyquist sampling as shooting fish in a bucket.

The bucket limits what is in the bucket by being only so large. This is like a band limit. We can only have so much variation.

Two samples is the very least number of samples we need to specify that 1 sec of 1 Hz sine wave, although intuition tells us we are on the edge.  Three samples for a second and half is a safer proposition.   How about 120 samples for 60 seconds of 1 Hz signal. If you look at the reconstruction formula for Nyquist sampling, it says you need a need an infinite number of samples to get that 1  Hz sine back.  The summation limits are plus and minus infinity.

In practice, we usually sample a bit faster than what Nyquist specifies.

There are circumstances where one can sample less than what Nyquist tells us and there is a lot of activity in an area known as sparse sampling.

To be honest, this was a terrible analogy. The mathematics is really more straightforward. Things like mathematical limits are hard to analogize.

It also turns out that a 1Hz sine wave is also a mathematical ideal.


Combining Analog Video Signals for Diversity Receiver System

(or maximal gain combining if a simple circuit exists to variably attenuator each signal).

Nobody does maximum gain combining in analog circuitry (at least as far as I know). That would require something to estimate your channel coefficient, calculate its inverse, and convert it to a complex gain for an analog element - after having probably already converted the received signals to digital signals, to even estimate the channel coefficient.
You simply do that math in your digital signal processing. A single multiplication doesn't "hurt" much.
Equal gain combining is indeed sometimes done in analog hardware, but it's then usually called "phase shifting for analog beamforming" and mostly done in analog, because the realization in digital domain is too expensive due to a very high signal bandwidth or antenna count, typically in radar systems. You're not in one of these use cases, far as I can tell.

My question is, how is this method practically implemented?

Honest answer: very rarely... the era of analog video transmissions ended before MIMO / diversity combining became a thing.
That can be explained, and it also hints at what you need to do when implementing this!
For any combining method to be useful, you need your video signal to take different paths to different antennas (otherwise: no diversity gain).
Now, since most video transmission things back when analog video was the dominant method were television, i.e. were done in an outside environment, there was a slight problem: The transmitter would typically be a tower that can be received with a mostly straight path by antennas on rooftops.
That doesn't really lend itself to any temporal diversity: that channel typically doesn't change over time.
Now, if there's more than one independent paths for multiple antennas attached to the same system, that typically means there's some big reflector (typically: mountain behind your receiver), and these paths overlay.
That introduces multipath into every antenna, typically: Let's say your analog TV signal was 6 MHz wide: at $c_0=3\cdot 10^8\,\frac{\text{m}}{\text{s}}$, that basically means that as son as the delay between different important paths comes into to the region of 50 m, you'll have a visible problem. For larger delays, it becomes such an important problem that it got a name: Ghosting.

source
You can see the echo!
Now, your combining system relies on the existence of multipath, effectively, because NTSC runs at frequencies and over distances where the antennas get to large to have multiple independent antennas in an indoor / handheld device. But if multipath occurs, this is what it looks like... Before gaining anything from diversity combining, in a real-world system, you'd first solve the multipath issue, because that's a severe problem.
Solving that effectively requires an equalizer. So you build an equalizer, which requires your receiver to know something about the transmitted signal, and then use that knowledge to estimate what a filter would look like that "reverses" the multipath.
Proving that this can be done is an early 1990's thing. Realizing that doing that on analog transmissions is more hassle than moving to digital signalling happened exactly at that point. So, to my knowledge, not a single analog video system with diversity combining ever hit a market.

Is the combining of the signal usually done before or after the receiver?

But now you know what you'd need to do: first cancel echos, i.e. have an equalizer, for each antenna. Your combining techniques don't (or at least, not in their vanilla form) work at all if the channel is more than a complex coefficient – and a multipath channel is a complex impulse response, so you need to solve that first. After you've built an equalizer for an analog signal, you can then combine them.
Note that it's really a pain to equalize analog signals, especially things that aren't well-shaped. And NTSC really isn't a great candidate for that, there's very little that you can rely on to give you a channel estimate. That means your analog transmission needs to be slightly extended to include some kind of channel sounding sequence - and that essentially means you're adding something from a digital world (a known preamble) into your analog signal.
Let's assume you do that. (hint: TVs with such techniques did exist... you likely didn't have one. Early on, that was probably a "premium" feature, and when it became affordable, analog TV was on its way out.)

Would the method of converting to digital and scaling gain introduce too much noise to the signal in practice?

You do realize any modern transceiver, be it your digital TV, your smartphone, your wifi access point or your bluetooth thermometer, does that? The noise introduced by an 8 bit ADC will be laughably low compared to the noise you got in analog TV reception, anyways.
This isn't the 1980s...

What kind of ADC and variable amplifier would be needed in this application (NTSC composite video)?

Now, this is a very hard question. The answer to that is: figure out the bandwidth of your signal, and build a complete receiver frontend (which is a pretty complicated thing at first) that can digitize that bandwidth (Nyquist tells you how much sampling rate you need; if this is news to you, you'll have to read up on digital signal processing first). But that's just a lower limit: what ADC you pick depends on how much effort you want to spend on analog filters, what your system's digital interfacing looks like, what skills you have with FPGAs (you can't just plug an ADC into your computer)....
Effort-wise, and probably even financially, you'd just buy an SDR device instead. In these bandwidth regions, things like hackRFs, bladeRFs, even LimeSDRs become pretty affordable; I personally like (but also: am affiliated with) the Ettus product line of USRPs better. You'll need two channels that are synchronously sampled.
Now, after that you'll need to invent a viable method of adding some channel-sounding sequence to your NTSC signal (standard NTSC won't do!), then you need design an equalizer in your digital signal processing, then you can throw your combining at things.
Uff. At that point, really, moving to digital video would have saved you a lot of effort: First of all, digital terrestrial TV systems solve the multipath issues, typically through multicarrier (specifically: OFDM) methods. That directly gives you something you can combine.
Now, look at the quality of images that NTSC offers: you get interleaved video with around 525 lines... After applying a very well-intending smoothing filter to make things look less ugly, something like this:

Honestly, for the same video quality, you'd have to use far, far less bandwidth if you threw digital video encoding at the problem, you'd get a digital receiver, which is way easier to deal with if you're equalizing, and it's kind of not technology that's not been dead as a doornail for 20 years; I'd drop the NTSC requirement, and plug together something narrower in spectrum using the DVB blocks that come with GNU Radio (which I'm also affiliated with).
That gives you digital transmitter whose output you can either plug into a transmitting SDR (see above) or you can use in a simulation (which I strongly recommend before dealing with RF, because that opens up a completely new box of problems, like non-reproducibility of channel realizations).
You can use GNU Radio then to simulate your different channels, and then combine the resulting signals, to your heart's delight. GNU Radio's DVB toolbox brings both the transmitting as well as the receiving side for digital video broadcasting standards (and arbitrary modifications of these standards), and you can omit e.g. the channel coding and decoding for your experiments.


How to use Expectation maximization to estimate poisson noise matlab

Explanations of the Expectation Maximisation algorithm:

Nature.com ( Nature Biotechnology 26, 897 - 899 (2008) )

Wikipedia - Expectation–maximization algorithm

Stack Overflow - What is an intuitive explanation of Expectation Maximization technique?


If you're not interested in specific practical aspects of A/D conversion, but if you want to learn basic theory concerning sampling and digital (discrete-time) processing of analog (continuous-time) signals, I'd recommend that you read and study the chapter Sampling of Continuous-Time Signals in Oppenheim en Schafer's book Discrete-Time Signal Processing.


Plenty! In fact, since the mid-1970s, people in security have been working on automatic transmission classification.
It's a hard problem, and you logically can't detect all transmission schemes, but it's something that's being sold as service and as product, especially for government (and three-letter government agencies).

These guys, for example, focus on HF/VHF/UHF comms, and I've seen the list of classifiable transmission schemes they had for a public demo. It's long.

Now, what these recognize are just common standards (in Ham Radio, that would be called modes),  and they are specific to one greater usage scenario (security-relevant communication).

There's a lot of underlying algorithms to automatically classify the number of observed channels, bandwidths, symbol rates, constellations, spreading, scramblers, error-correction codes, frame lengths, synchronization preambles.... but even listing these would take a book.

I can point you to what a student did during Google summer of code: It's gr-inspector, a GNU Radio-based signal detection and classification program.

It's pretty generic in that it detects bands with power in them (which is no good for spread-spectrum systems, but for anything else it's a viable approach), and it gets pretty specific once you have an OFDM signal in there. Part of his job was writing a classifier just to get OFDM parameters: carriers, spacing, length of cyclic prefix. That took someone who was dangerously close to being a M.Sc. in communications engineering nearly 3 months. You can imagine how much work it is writing a set of classifiers for all useful transmission schemes!

You might be interested in the two papers he cites in his blogpost I linked to above; and especially their references.


How is a signal with SNR less than unity demodulated?

The necessary condition is that the techniques transmit at the rate less than channel capacity.

For example in AWGN channel, for complex (two dimensions of freedom) signal

$$\mathrm{C} = \log_2(1+\mathrm{SNR}) \leq_{(SNR < 1)} 1 \quad \textrm{bit per channel use}$$

A channel use is, in most cases, the generalized term of modulation symbol.

In fading channels, channel capacity (if it exists) formula may be different but the idea is the same.

Thus other than the techniques to increase SNR such as beamforming, there are two principal technique families to deal with the low SNR condition

Channel coding (the usage of pilot and equalizer, diversity MIMO, etc, can be considered as a form of coding)
Spead spectrum : this can also be considered as coding by taking advantage of multipath diversity


Difference Between monocular video and stereo video?

Monocular video is taken using a single camera.  Stereo video is taken using two cameras side-by-side, producing two video streams.  That's pretty much it.

A stereo video pair allows depth measurement, so you can find out how far away things are.  However, this is not an easy task for complex scenes.


In general you can divide the audio books into the following categories:

Audio Engineering
Audio Programming
Audio Signal Processing and Coding
Digital Signal processing

For your exclusive requirement (that fall into cat-2) I highly recommend the following book: The Audio Programming Book_BOULANGER_MIT


This may sound trivial but the sampling theorem can only be violated in the process of sampling an analog signal. In your question you would like to sample the signal
$$
x(t) = \sum_{k=0}^{N-1}\mathrm{e}^{\displaystyle -j 2\pi k\Delta f t}
$$
and you state correclty that, in order to not violate the sampling theorem, $f_\mathrm{s}>2(N-1)\Delta f$ must hold. In contrast, OFDM systems usually work with a sampling frequency $f_\mathrm{s}=N\Delta f$ and you may now ask where this discrepancy originates.

If we look at $x(t)$ we notice that it has a one-sided spectrum with discrete frequencies at $0, \Delta f, \ldots, (N-1)\Delta f$. That's an odd choice because the signal bandwidth is $2(N-1)\Delta f$ but only $N$ subcarriers are used. With other words: this system has potentially $2N$ subcarriers but the subcarriers $-N, -N+1, \ldots, -1$ are unused and the corresponding spectrum is wasted. We should rather define $x(t)$ like this:
$$
x(t) = \sum_{k=-N/2}^{N/2-1}\mathrm{e}^{\displaystyle -j 2\pi k\Delta f t}
$$
The absolute maximum frequency is now $\frac{N}{2}\Delta f$ and accordingly $f_\mathrm{s}>N\Delta f$.

For creating a discrete-time OFDM signal, however, $x(t)$ is quite irrelevant. You cannot violate the sampling theorem when generating a discrete-time signal, because it can be considered as already sampled. Furthermore, due to the cyclic property of the DFT, the limits of the sum in the definition of $x[n]$ are correct. Whether they are $0 \ldots N-1$, $-N/2\ldots N/2-1$ or $N+3\ldots 2N+2$ doesn't matter.

So I wouldn't say that the textbook is wrong. But if OFDM is motivated by sampling an analog signal, then $x(t)$ as defined in your answer is a poor choice.


Does changing the frequency of a wave cause a change of the phase of the same wave?

Let's assume your two sinusoids are at the same phase, $\phi_0$, at $t=0$ and you change the frequency of one of them at exactly that time.

The phases as functions of time, for $t \ge 0$, for those sinusoids at the two different frequencies are

$$\phi_1(t) = 2\pi f_1 t + \phi_0$$
$$\phi_2(t) = 2\pi f_2 t + \phi_0$$

And the difference between those phases is

$$\psi(t) = \phi_1(t)-\phi_2(t) = 2\pi(f_1-f_2)t$$

So the phase difference between them is a function of time and the difference between the two frequencies.


Autocorrelation of a signal comprised of only the same number

What you see are mainly numerical effects, even though the second result is actually correct, but not according to what Matlab intends to compute. But let's go slow ...

Matlab computes the autocorrelation function by first subtracting the mean of the signal and then correlating it with itself. Then the autocorrelation function is normalized such that its value for lag 0 equals 1 (this is where you get the NaNs from).

What happens in your case is that the signal equals its mean (because it is constant), so by subtracting the mean you generate a zero signal, the autocorrelation of which is zero. Then in the normalization step you divide the zero signal by its first value (which is zero), and you get NaNs. This shouldn't happen in a well written function, so you could complain with Mathworks about it.

What is even worse, is that this behavior is not consistent. You found some scaling values for which you actually get a reasonably looking autocorrelation function. This is really weird but I think I understand what's going on. Try the following. Look at the values of

signal - mean(signal)

This is what Matlab bases its computiation on. What you'll see is that for the scaling values which give you NaNs, this difference is exactly zero (which is the correct result). For the other scaling values (which give a well-behaved autocorrelation function), this difference will be some constant in the order of $10^{-17}$. So you're again left with a constant signal, albeit it much smaller than before. The autocorrelation function of such a constant signal (without subtracting the mean, because this has already been done) is triangular, simply because it is the convolution of two rectangular functions. The normalization gets rid of the $10^{-17}$ factor and you actually get a nice looking autocorrelation of your signal. This one is even correct if you were to compute the autocorrelation without subtracting the mean of the signal (but this is not what the Matlab function intends to do!).

Conclusion: the first result (NaNs) is useless but at least consistent with the definition of the autocorrelation function used by Matlab (even though a division by zero should be prevented in a good program). The second result is not consistent with Matlab's definition of the autocorrelation, but it is - due to numerical errors (!) - correct when the autocorrelation is defined without first subtracting the mean of the signal. Note that both definitions are used in practice so it's always important to check what is actually meant by 'autocorrelation'.


Biquad cookbook formula for broadcast FM de-emphasis?

You could also implement a 2nd order IIR filter which has response of 1st order filter by using some other TF but one used in RBJ's cookbook.
Here's one way to do it using MZT (C++).

double Fs = 44100; // higher sample rate results better filter accuracy
double T1 = 0.000075; // time constant
double p1 = -exp(-1.0/(Fs*T1));
double z1 = 1 + z1;

a0 = 1.0;
a1 = p1;
a2 = 0;

b0 = 1.0;
b1 = z1;
b2 = 0;

% swap between a1, b1 to select pre- or de-emphasis

Coefficients from above calculations adds +13.63dB gain so you need to normalize before use.

Same as Octave code:

% PACKAGES
pkg load control
pkg load signal

clear all;
clc;

fs = 44100;
samplingtime = 1/fs;

% analog prototype
A2 = [1];
B2 = [0.000075 1];
Ds = tf(B2, A2); % swap between A2 and B2 to select pre- or de-emphasis
Ds = Ds/dcgain(Ds);

% MZT
T1 = 0.000075; % 75us
z1 = -exp(-1.0/(fs*T1));
p1 = 1+z1;

a0 = 1.0;
a1 = p1;
a2 = 0;

b0 = 1.0;
b1 = z1;
b2 = 0;

% swap between a1, b1 to select pre- or de-emphasis

Bmzt = [b0 a1 b2]
Amzt = [a0 b1 a2]

DzMZT = tf(Amzt, Bmzt, samplingtime);
DzMZT = DzMZT/dcgain(DzMZT);

%% Plot
wmin = 2 * pi * 20.0; % 20Hz
wmax = 2 * pi * ((fs/2.0) - (fs/2 - 20000)); 20kHz

figure(1);
bode(Ds, 'b', DzMZT, 'c', {wmin, wmax});
legend('Analog prototype', 'MZT 2nd order','location', 'northwest');
grid on;

There are few other ways to calculate the coefficients as like what mentioned here.

If you don't want to code then here are (75us) IIR coefficients for most common Fs:

44100
den = [ 1 -0.2430506562509079 -0.3653116688328392 ]
num = [ 0.226236781592716 0.1524695626369773 0.01293133068655965 ]

48000
den = [ 1 -0.2625656192085276 -0.3736509425970024 ]
num = [ 0.2102915623850396 0.1415200473180778 0.01197182849135262 ]

88200
den = [ 1 -0.3702185951814965 -0.419953874650389 ]
num = [ 0.1217100727211314 0.08132561501764404 0.006791842429339071 ]

96000
den = [ 1 -0.3813262768034237 -0.424764645791489 ]
num = [ 0.1125121591800709 0.07513002079134688 0.006266897433669414 ]

176400
den = [ 1 -0.4405600251706213 -0.4505357896180226 ]
num = [ 0.06328932872397297 0.04212236025403318 0.003492496233349943 ]

192000
den = [ 1 -0.4464885169088165 -0.4531262255096845 ]
num = [ 0.05834721577148291 0.03882109525679095 0.003216946553225094 ]

352800
den = [ 1 -0.4775710488916642 -0.4667412351991018 ]
num = [ 0.03239176683609103 0.02151798489625147 0.001777964176891494 ]

384000
den = [ 1 -0.4806342452619782 -0.4680860610215655 ]
num = [ 0.02982988806358219 0.01981315821469544 0.00163664743817871 ]

Magnitude error less than 0.1 dB (20Hz-20kHz)


The phase vocoder algorithm used with overlapped FFT frames might easily catch certain kinds of discontinuities in pure sine waves.

If the audio file is compressed, note that this process itself may add certain kinds of distortion.  So you would need to find an experimental threshold.


Difference between double-sided and single-sided AWGN noise after bandpass filtering?

Both descriptions give the same result. We often use the one-sided noise power spectral density (PSD) because for real-valued processes the negative frequencies are redundant, so defining the PSD for positive frequencies is sufficient. You just have to scale the noise power spectrum such that integrating the one-sided power spectral density (PSD) over the positive frequencies gives the same result as integrating the two-sided PSD over positive and negative frequencies. I.e., for white noise with a constant PSD, defining the one-sided PSD as $N_0$ means that we would need to define the two-sided PSD as $N_0/2$ such that the noise power remains the same: $N_0W=N_0/2\cdot 2W$


The butter command expects a normalized frequency from 0 to 1, where 1 is the Nyquist frequency.  In other words, Wn is in half-cycles per sample.

So you just divide your desired frequency (1 cycle / 48 hours) by half the sample rate ((1 sample per hour)/2), giving Wn = 0.04167 half-cycles per sample.


How to filter sudden spikes and baseline changes in microphone signal?

Assuming some processing delay is acceptable (100 ms should be more than sufficient), and given there is such a substantial difference between the spike peaks and the signal level, a robust solution would be a non-linear puncture based filter by removing the spikes based on a threshold detector (that then replaces the the required number of samples before and after the threshold with a prior running mean of the filtered signal) and then pass the resulting signal through a low pass filter that has a duration in time at least ten times longer than the interval of time removed due to a spike. Use a running mean and rms of the filtered output for dynamic processing. The mean is described above to replace removed samples as a reasonable minimum error estimate of the sample otherwise hidden and the standard deviation is used to set the threshold. A threshold of 3-sigma above the mean should work well, but ultimately you could also track the statistics of the peak and place the threshold based on optimized probability of false alarm versus probability of detection trades.


How to Detect Gradients and Edges in Images?

I was working in opencv and trying to find the peak of a gradient generated by distance transform. I realised that using morphologic operations (erosion/dilatation) in grey-scsale images was very useful in this case.
If you erode dilate a grey-scale image, any pixel is going to take the value of the lower/highest neighbour. You can therefore find peaks of intensity in gradients by subtracting grey-scale image from the same dilated/eroded image. Here is my result:

And a way to do it in OpenCV/Cpp:

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"

int main( int argc, char** argv ){

    cv::Mat objects, img ,peaks,BGR;
    std::vector<std::vector<cv::Point> > contours;
    /* Reads the image*/
    BGR=cv::imread(argv[1]);
    /* Converts it to Grayscale*/
    cv::cvtColor(BGR,img,CV_BGR2GRAY);
    /* Devine where are the objects*/
    cv::threshold(img,objects,0,255,cv::THRESH_BINARY);
    /* In order to find the local maxima, "distance"
     * is subtracted from the result of the dilatation of
     * "distance". All the peaks keep the save value */
    cv::dilate(img,peaks,cv::Mat(),cv::Point(-1,-1),3);
    cv::dilate(objects,objects,cv::Mat(),cv::Point(-1,-1),3);

    /* Now all the peaks should be exactely 0*/
    peaks=peaks-img;

    /* And the non-peaks 255*/
    cv::threshold(peaks,peaks,0,255,cv::THRESH_BINARY);
    peaks.convertTo(peaks,CV_8U);

    /* Only the zero values of "peaks" that are non-zero
     * in "objects" are the real peaks*/
    cv::bitwise_xor(peaks,objects,peaks);

    /* The peaks that are distant from less than
     * 2 pixels are merged by dilatation */
    cv::dilate(peaks,peaks,cv::Mat(),cv::Point(-1,-1),1);

    /* In order to map the peaks, findContours() is used.
     * The results are stored in "contours" */
    cv::findContours(peaks, contours, CV_RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE);
    /* just draw them and save the image */
    cv::drawContours(BGR,contours,-1,cv::Scalar(255,0,0),-1);
    cv::imwrite("result.png",BGR);

    return 1;
}


How to apply a 1D filter (CSF) to an image?

I believe what you want to do is take your 1D filter column vector, call it $x$, and create a 2D filter with $xx^T$.  The result of this outer product (column vector times row vector) is a matrix that you can use to do element-wise multiplication with the 2D Fourier transform of the image.


What is the difference between phase delay and group delay?

First of all the definitions are different:

Phase delay: (the negative of) Phase divided by frequency
Group delay: (the negative of) First derivative of phase vs frequency

In words that means:

Phase delay: Phase angle at this point in frequency
Group delay: Rate of change of the phase around this point in frequency.

When to use one or the other really depends on your application. The classical application for group delay is modulated sine waves, for example AM radio. The time that it takes for the modulation signal to get through the system is given by the group delay not by the phase delay. Another audio example could be a kick drum: This is mostly a modulated sine wave so if you want to determine how much the kick drum will be delayed (and potentially smeared out in time) the group delay is the way to look at it.


How do we use signal::filter and signal::butter in R for EEG data?

I think your problem is that you want to use signal::filtfilt rather than signal::filter for the cut-offs you set to map on to the half-amplitude cut-off (otherwise, with the signal::filter function, these are the half-power cut-offs).

See my answer here for example code showing how to use the W parameter in signal::butter and double-check that your filter is working appropriately:

https://stackoverflow.com/a/49905340/9665680


What are the Elements of Filters for Highly Volatility Price Series with Jumps?

In case you're looking for an edge-preserving low pass filter, you could have a look at the bilateral filter.

Implementation is fairly straightforward.


Measuring noise floor in an impulse response

I can suggest you take a look at the supporting information for Dirac software, where they describe the process of INR (IR to Noise Ratio) calculation.
The simplest approach is to estimate the noise level based on the beginning of the IR before the main peak. For that purpose, you can also do the tail of the Impulse response.

Yet another article supporting that estimation is by Lundeby et al.: "Uncertainties of Measurements in Room Acoustics". Their algorithm is an iterative procedure which tries to estimate the background noise level and the slope:

You could also try to verify that the noise floor doesn't have any distortions, by fitting the linear regression and making sure that the slope is as close to 0 as possible.


What Is the Correct Way to Apply Richarson Lucy Deconvolution to Luminance Data?

I think the answer is in Jim Kasson's Post (The Zone System and Digital Cameras) you linked to:

If you are only moderately familiar with the Zone System, you may be asking yourself about the utility of recording light brighter than 100% reflectance. The 100% reflectance calibration assumes a matte subject with perfectly Lambertian reflectance. The higher light levels occur when the subject has some specularity or gloss to it. Water is a classic, and photographically-important example. Metals, glass, plastics, skin, stone and many paints are also often far from Lambertian. So we often need the headroom to accommodate those materials.
All of the above assumes that the camera is calibrated using UniWB, and that the UniWB hack works perfectly to predict raw file values. However, UniWB never achieved much penetration among digital camera users, and I think the number of users is waning today as cameras get more and more dynamic range, so it’s worthwhile looking at what happens with widely-used camera setups.

So the model to of Reflectance and the actual number of photons the sensors is exposed to are linear only if the camera is behaving according to UniWB Model (See also Guillermo Luijk - UniWB Tutorial) and the material used is indeed perfectly Lambertian.
Both probably doesn't hold and create some non linear effects.
CMOS / CCD Sensors are indeed linear in their response (At least before they are clipped).
So, given a RAW data with Linear Demosaicing (Pay attention that some algorithms may be non linear in this step) you can indeed apply the Lucy Richardson Convolution on the data.
Pay attention that LAB Color Space is a non linear transformation of the RGB colors (even without the Gamma Correction).
So I'd apply the deconvolution in the following order:

Demosaic RAW Data.
Convert data into YCbCr Color Space (Or similar) which decompose the colors into luminosity and color data with Linear Transformation.
Apply the deconvolution on the Y channel.
Convert back to RGB based color space.


How do I set window length and shift length for computing STFT on Scipy?

for this function:

nperseg will set the number of samples to output from your STFT.
noverlap will specify how many old samples you want to take to calculate your future STFT.

Let's say you want a stft of size 256 with an overlap of 50% you will set noverlap = 128.
128 samples will then be taken from the old segment that was used to calculate the previous STFT...
For more information about STFT and the interest of overlap I invite you to watch this discussion.
STFT: why overlapping the window?