Is it possible to encrypt an audio wave? NFL Example

Someone mentioned how occasionally, these radios will pick up interference from airplanes, and that frequently the helmets lose communication.

Is this the 1970s?

Modern headsets should be pretty robust against interference, and they have a dedicated spectrum allocation. If there's interference from airplanes, either device malfunctions.

So, no, this sounds like rumors, or long obsolete claims, especially since for superposition reasons and 6 decades of compatibility, aircraft voice comms is partly still AM; modern (as in: post-1970) headsets are either digital or at the very least FM.

  One thing I wonder about is if an opposing team's coaches could pick up what was being communicated between the offensive coordinator and the quarterback.

Sure. It's on the air. Radio waves can't tell friend from foe. And since a quarterback can be all over the field, there's certainly no high-gain antenna focusing the radio beam on the player alone.

  And then I wonder, well, why wouldn't there just be an encryption scheme between the coach and the quarterback.

Just because you can receive the mobile phone call of your neighbor, or his WiFi, or pay TV, doesn't mean you can decrypt stuff.

So, yes, there's plenty encryption systems in the technological fields of voice comms, some civil, some military... plenty to pick from.

  But then I wonder if the encrypted signal would be distorted with all this noise, and if that distortion would completely disrupt the decrypted signal.

Digitally, whether or not the signal is encrypted makes little difference to the signal quality in case of errors. You typically break one full "frame" worth of channel coding / encryption.

All wireless digital communication schemes incorporate channel coding that allows for forward error correction (FEC) so that this makes no difference, and the audio codecs (specifically: the VoCoder) for such communication channels are robust against "surviving" errors.

  Are there any methods for this kind of "loose encryption"?

I'll not go into crypto here, because a) that's not my forte and b) because this isn't the perfect platform, but:

No, that's something a cryptographer must never implement; every encryption must have zero correlation between encrypted and decrypted bits. If that wasn't the case, the cryptography would be inherently insecure, because you could, through statistics, relatively easily reconstruct (parts of) the unencrypted communication.

There are "human" encryption schemes. Like NFL coaches could change the name of their plays from game to game. This would keep it known to only the QB and the coach.

  There are "human" encryption schemes. Like NFL coaches could change the name of their plays from game to game. This would keep it known to only the QB and the coach.

That would be called a "substitution cipher", and these are extremely weak. Especially, you call a specific move "rabbit", and then the quarterback does a move, for the rest of the game it's clear what "rabbit" means.

Commentary on the situation: Obviously, the solution to this problem is organizational, not technological: It's a sport. If you cheat, you get disqualified. So, there's strong incentive not to take performance-enhancing drugs. And there's just as much incentive to not use information you mustn't get.

Just that the supervision of the former requires invasive, expensive testing, and the question what kind of knowledge a coach had can be relatively easy tested within similar margins of error as drug tests with basic statistics.


Difference between 2D DFT's and 1D DFT's of Linearized Matrices

The DFT is "Separable Operator" (Also the classic Fourier Transform) and hence can be applied on the Rows and Columns of the image separately (It can be generalized to N dimension and not only 2).
Yet, if you create  1D signal from your image (Let's say by Column Stack) and apply 1D DFT you don't get the information you would by using 2D DFT (By going on the Row and them Columns).

Remember, Fourier Transform is all about synthesizing the signal using different functions. In this case if it is 2D signal you want to build it using 2D Signals.

You need to play with the indices to first apply it on each row and later on each column.


Finding the Harmonics contained in an Image?

this questions is asked fairly often. If you don't have a grasp of fft in 1d, higher dimensions may be difficult to grasph. But logically it makes a lot of sense once you see it What does frequency domain denote in case of images?

basically sharp changes (high contast noise) usually has high frequencies. and smoother gradients have low frequencies. just as a quick crash course

edit

sorry i thought the other post had code, its just code agnostic conceptual stuff. anyway ill tell you matlab (since your tag says matlab in it) and it is more familiar to me than openCV

function fft_im = createFFTIm(grayscale_im)
    fft_im = fft2(grayscale_im);

    figure()
    subplot(2,2,1:2);imshow(grayscale_im); title('original');
    %fftshift puts the 0 frequency in the middle of the plotting area, normally
    %it would be split, part of it in the 4 corners of the image
    %we add 1 because log(0) is -inf, which we dont want
    %we scale it jsut to show the colors better
    subplot(2,2,3);imagesc(log(1+abs(fftshift(fft_im)))); colormap(gray);
    title('magnitude spectrum');

    subplot(2,2,4);imagesc(angle(fft_im));  colormap(gray);
    title('phase spectrum');
end

edit #2 i slightly modified the code

and here are two examples

notice with the chess example, there are visible peaks in the magnitude spectrum. Just like the regular fft, these peaks correspond to a sing wave of a specific frequency. Since out image is periodic there is are strong sines at the corresponding frequencies. You also notice since the image is near symmetric, the phase is also periodic.
The phase is responsible for modeling the fact that an image is not symmetric

edit #3
this decomposes the figure into individual elements. In the code I limited it to a 4x4 subblock of pixels (the upper left portion of the figure) but you can change this all the way up to the size of the image (though for any large image this would be incredibly slow)

As in my comment, I processed the image one pixel at a time, and did an FFT2 on each of those 1 pixel images. The images themselves are still the same size as the original, it just has only one pixel value populated

function fft_im = decomposedcreateFFTIm(grayscale_im)
    %i only did a small 16 pixel sub block of the image, 4x4 pixels from
    %the upper left corner of the image
    max_rows_to_process = 4;
    max_cols_to_process = 4;

    [m,n] = size(grayscale_im);

    %this is used to isolate pixels, for individual analysis
    zer_array = zeros(m,n);

    figure(1);title('magnitude')
    figure(2);title('phase')

    %i only did a small 16 pixel sub block of the image, 4x4 pixels from
    %the upper right corner of the image
    m=4;n=4;

    %for all rows
    for ii=1:1:max_rows_to_process
        %for all columns
        for jj=1:1:max_cols_to_process
            curr_idx = (ii-1)*max_cols_to_process + jj;

            %creates an image with only one pixel
            indiv_pix_pic = zer_array;
            indiv_pix_pic(ii,jj) = grayscale_im(ii,jj);

            %does fft
            fft_im = fft2(indiv_pix_pic);

            %fftshift puts the 0 frequency in the middle of the plotting area, normally
            %it would be split, part of it in the 4 corners of the image
            %we add 1 because log(0) is -inf, which we dont want
            %we scale it jsut to show the colors better
            figure(2)
            subplot(max_rows_to_process,max_cols_to_process,curr_idx);
            imagesc(angle(fft_im));  colormap(gray);

            figure(1)
            subplot(max_rows_to_process,max_cols_to_process,curr_idx);
            imagesc(log(1+abs(fftshift(fft_im)))); colormap(gray);
        end
    end
end

magnitude

phase


Method to generate binary sequences with desired cross-correlation and autocorrelation properties

Flattered though I am that a short note (D.V. Sarwate, "Bounds on crosscorrelation and autocorrelation of sequences", IEEE Trans. Information Theory, vol. IT-25, pp. 720-724, 1979) that I wrote over 40 years ago is still receiving some notice, the OP's question is way too broad for this forum; people (not me!) have written whole books on this subject.
Neither Lloyd Welch nor I ever claimed that there exist sequence sets that meet these bounds exactly or almost exactly. Rather, every set of $K$ sequences of period $N$ must be such that $(1)$ and $(2)$ in the OP's question are satisfied. In particular, the OP's belief

From this, I understand that it is possible to find two sequences whose cross-correlation is ideal ($\theta_c = 0$) and whose autocorrelation is good enough ($\theta_a \approx \sqrt{N}$).

is a false understanding. If two sequences are uncorrelated (that is, $\theta_c = 0$), then $\theta_a$ must be at least $\sqrt N$, but it is not the case that there must exist a pair of uncorrelated sequences with $\theta_a \approx \sqrt N$. My paper cited above gives examples of a set of $N$ sequences (the rows of a length-$N$ DFT matrix) for which $\theta_c = 0$ but $\theta_a = N$ which is much larger than the desired $\sqrt N$. Unfortunately, while smaller subsets do still satisfy $\theta_c = 0$, the maximum off-peak autocorrelation continues to cling to $\theta_a = N$. The paper also described the Frank-Zadoff-Chu sequences for which $\theta_c = \sqrt{N}$ and $\theta_a = 0$. Both sets are optimal with respect to these bounds.  The sequences in both sets are complex-valued, not binary.
For sets of $N$ binary sequences of period $N$, a bound due to Sidel'nikov tells us that $\theta_\max$ must be at least as large as $\approx \sqrt{2N-2}$. $\,\,$ Gold sequences of period $N = 2^{2n+1}-1$ essentially meet this bound ($\theta_\max $ equals $\sqrt{2(N+1)}+1)\approx \sqrt{2N-2}$, but for periods $N = 2^{2n}-1$, $\theta_\max$ is $\approx 2\sqrt{N}$ for Gold or Gold-like sequences. If $\theta_\max\approx \sqrt{N}$ is non-negotiable as a criterion for sequence design, then for $N=2^{2n}-1$, the small set of Kasami sequences has $\sqrt{N+1} = 2^n$ sequences with $\theta_\max= \sqrt{N+1}+1 = 2^n+1$. For more details, see D.V. Sarwate and M.B. Pursley, "Cross-correlation properties of pseudorandom and related sequences," Proc. IEEE, vol.68, pp.593-619, May 1980.
With regard to other queries in the OP's question, some work has been done towards the construction of of sequences whose autocorrelation values are small for small nonzero offsets (while larger bounds for large offsets are OK), but typically, small crosscorrelation values for small offsets while allowing larger values for larger offsets has not received much attention. The former problem is of interest because it makes fine-tuning of the synchronization process in spread-spectrum systems easier; the latter might be more of a niche problem for special systems with nearly synchronized times-of-arrival.


Multiply signal $x[k]$ with $\cos(2\pi\nu_0k)$, then given $X(\nu)$ draw resulting function in frequency domain?

Thinking in terms of convolution with shifted impulses help. Multiplication in the time domain corresponds to convolution in the frequency domain. Your example is a classic showcase of frequency/band shifting using a "carrier" such as a cosine (or sine) or a complex exponential. Now, note, that the Fourier transform of these carriers are actually shifted impulses. When you convolve a function x(t) with a shifted impulse, say, D(t-T), the result that pops out is x(t-T), i.e., basically a shifted form of x(t). Note that the amount of shift is exactly equal to the shift in the impulse (the actual math is quite easy). Likewise, in the frequency domain, convolution with a shifted impulse results in your Fourier transform of the original signal getting shifted accordingly - For an impulse shifted to the right of origin, it feels as if the impulse "pushes all the bands in the spectrum to the right" exactly by the amount it itself is shifted from the origin. Vice-versa analogy holds for left shifted impulses. This is the mathematical background of this if you needed one. And yes, those results look correct to me. Cheers!


I've tried to get the bin with greatest magnitude but that only give me right results for higher pitch signals, it doesn't matter which oversampling factor I use I still get bad data for low freq signals.

That's because the harmonics are larger than the fundamental.  Plot your spectrum and you'll see.  A better method to find the true fundamental is autocorrelation.  Then you're "sliding" the waveform past itself and finding delays at which the wave shape lines up with itself.

http://www.phon.ucl.ac.uk/courses/spsci/matlab/lect10.html

Do you really want them to sing the exact note, or is it ok if they sing an octave above or below depending on their voice register?


Fast Recursive 1D Signal Smoothing - IIR / Auto Regressive Implementation of Gaussian Smoothing

Have a look at my Fast Gaussian Blur Project at GitHub.

You will find there implementation of IIR Approximation of Gaussian Blur which implements the following papaers:

Recursive Gabor Filtering.
Recursive Implementation of the Gaussian Filter.
Boundary Conditions for Young - van Vliet Recursive Filtering.

The idea is pretty straight forward.


Any references for the trade-off between software and mechanics/optics in industrial machine vision systems?

I have found some "proverbs" like:

  Never use software to compensate for a poor lighting system.
  It is not cost effective and will result in a poor system design.

  It is cheaper to add a light-proof shroud to keep sun-light away from
  the object under inspection than to modify the software. Another
  universal truth which is often forgotten.

  Nothing exceeds the speed of light. Any processing that can be done
  optically will save a lot of computer processing later.

in the book "Intelligent Vision Systems for Industry" by Bruce G. Batchelor and Paul F. Whelan and also in B.G. Batchelor and P.F. Whelan (1994), "Machine vision systems: Proverbs, principles, prejudices and priorities", Proceedings of the SPIE - The International Society for Optical Engineering, Vol. 2347 - Machine Vision Applications, Architectures, and Systems Integration III, Boston (USA), pp 374-
383. (see here http://elm.eeng.dcu.ie/~whelanp/proverbs/proverbs.pdf).

The proverbs are also in the 2012 book "Machine Vision Handbook", Editors: Bruce G. Batchelor ISBN: 978-1-84996-168-4.


frequencies in sound: multiple possibilities?

The Discrete Fourier Transform, or DFT (the FFT is an algorithm that computes the DFT) of a length of finite duration (which any practical transform would need to be) is identical to the result of the transform of an infinitely long sequence formed by repeating the original sequence in time.

Knowing this provides the following insights related to your question:

First, anything that repeats (exact repetition) in time can only exist at discrete frequencies in the frequency domain. We see this with the Fourier Series Expansion specifically as shown in the graphic below.

The concept behind the Fourier Series Expansion is that any single valued continuous function can be represented as a sum of sinusoidal components and notably each component MUST have a frequency that is an integer of 1/T where T is the duration of the signal in time. Therefore IF those sinusoidal components were allowed to play out for all time (rather than being bound to the time interval [0,T]), the next cycle immediately after T would have to commence and proceed exactly as the waveform did at the start of the sequence (as each sinusoidal component would do the same). Thus it is often described that the Fourier Series Expansion decomposes any periodic function into a sum of sines and cosines (or equivalently and I believe mathematically simpler, complex exponentials).

The mathematical model of the time limited signal from [0,N-1] (for N samples of the DFT, which is equivalent to the analog bound of [0,T] in time), to also be a signal extending to infinity repeating with time provides for further intuition into the behavior and result of the DFT. Specifically we see that repeating a signal in time results in discrete uniform spaced impulses in the frequency domain, that can only exist at integer multiples of 1/T (including 0) where T is the length of the base waveform in time.

The second characteristic of the DFT is that it is done on a waveform that is sampled in time.  Without going into significant detail, sampling in time is associated with repetition in frequency (for those familiar with A/D and D/A conversion this will be readily apparent). So with the DFT we have both characteristics of repetition in time and sampling in time, which therefore means we will have "sampling" in frequency (discrete frequencies where only non-zero values can exist) and repetition in frequency.

The repetition is an implied construct that actually helps significantly to provide an intuitive understanding of many signal processing constructs- especially when considering both analog and digital domains. To be clear, the DFT involves only a fixed duration sequence both in the time and frequency domains, but mathematically these sequences can be repeated for infinite duration with the same result.

For example, to your question of what happens where there is a partial cycle? If we realize this equally represents a waveform repeating for all time, we see that in this case there would be an abrupt transition in the waveform.

Such a waveform cannot be created or represented with a single sinusoidal tone. Going back to the Fourier Series Expansion, we know that it can be represented by multiple tones as long as they are spaced at integer multiples of the repetition rate. So where we thought one tone exists, similarly in the DFT there will be multiple tones as required to create such an abrupt transition.  According to the DFT these frequencies really exist as what we are solving for in that process is the frequency components that are needed to create the time limited time domain waveform.

Below shows the case above (lower plot) compared to what we would get if there was a complete integer number of cycles over the time duration used (upper plot).  This is one explanation of "Spectral Leakage" with the DFT but given as example insight into the relationship between the time duration of the DFT chosen and the frequencies that would result.

If the waveform was changing with time (either in frequency, or amplitude beyond the sinusoidal component itself, such as an envelope) this would require many frequency components to represent it. This is no different than a modulation view of the time domain waveform: to transmit signals over the air we modulate the amplitude or frequency of a carrier frequency (that is better suited to go over the air) which results in several frequencies being present around that carrier. If the waveform is not actually repeating with time (which is likely the case), then instead of discrete tones we will get a continuous band of frequencies around the carrier.


IQ signal FM descrimination using phase derivative

I simply use unwrapped atan2(IQ(i)) - atan2(IQ(i-1)) to estimate a discrete derivative, then low pass filter to below 15 kHz.  Although with a shallow slope, the 1st order approximation to atan() given by Boschen will work just as well.

Your noise might be due to not unwrapping the phase delta, or to not low pass filtering after doing the phase discrimination.  Not sure whether your angle() function is the same as atan2().


My 3x3 Weighted Median Filter doesn't seem to be improving my image quality

Apparently, you intend to weight the patch pixel per pixel in:
a1_R= W*patch_selected_R

and then apply a median on the piece-wise product, and divide by 15. By the way, beware of the * product and the input type of the array, you don't want a matrix product.
This is not the weighted median I know of. To me, the elements in  patch_selected_R should be duplicated, triplicated, $n$-plicated with respect to the corresponding integer weight in the mask. On a smallest example: if $w=[1,3,1]$ and $m=[1 4 3]$, the intermediate (longer, of length $1+3+1=5$) buffer is (before sorting):
$$b=[1,4,4,4,3]$$
thus $[1,3,4,4,4]$ after sorting, and the weighted median is $4$ (instead of $3$ with the classical median).
The resulting array here is indeed 15 times larger, and you then compute the median, and replace the central value, with no need to further divide by 15. I am thinking of an interpretation of your result. Meanwhile, the theoretical origin of the weighted median is describe here:
What is the advantage of weighted median filter over median filter?
Note: there are combinations of linear and nonlinear filters, like the mean-medians filters.


It depends on which DCT (there are more than four of them), the idea is to take your original data and append to it a mirror-reflection copy of that data.  Now your data is twice as long, and you perform an FFT on twice the length.  but you have some symmetry to it that allows you to toss half of the data in the result of the DFT.

exactly how, depends on if it's DCT I, DCT II, DCT III, DCT IV, or the "MDCT" or whatever is the DCT flavor of the month.

looks like something useful here.  it may need further exploration.


Distort voice to hide talker identity as much as possible

To preserve the talker's identity for privacy reasons, I would make available for listening a distorted version of the original audio file, while using the original one for real recognition to improve the system.

This won't work.

The point is that the whole job of your speech recognition algorithm is to figure out which "components" in your audio signal are of significance, and what they mean, and to ignore things that are irrelevant to the textual content, namely identity of the speaker.

By wanting an algorithm that removes a dominant aspect of the signal but preserves the original speech content, you'd first need to have an algorithm able to tell the content-wise important components from the unimportant. That would be the major part of your speech recognition algorithm, which is exactly what you wanted to build in the first place.

Generally, this sounds like you want to do something that other companies literally spent several millions on with a small team and little experience. I'm afraid that's not really a promising approach.


How would I fit a distribution to this image noise?

Since you say your distribution is discrete, if your samples are independent, you could try a multinomial distribution as a model
https://en.m.wikipedia.org/wiki/Multinomial_distribution

You can use either maximum likelihood or a Bayesian estimator.

You could also take the mean and standard deviation and use a Gaussian model with rounding.  A goodness of fit test

https://en.m.wikipedia.org/wiki/Goodness_of_fit

would tell you how good the model is.

If your samples are independent noise generation is straightforward .  A Bayesian model is different because you use another (prior) noise generator as an input to the fixed parameter noise generator.

If your samples are correlated you can use a Gibbs sampler.

You should perform an independence test on your data, of which a Google search will provide a number of candidates

You also might consider boot-strapping

https://en.wikipedia.org/wiki/Bootstrapping_(statistics)

It is also very likely that a physics based noise model exists for your device.


Given a discrete time signal, what is the sequence of possible frequencies I can get from DTFT?

The discrete-time Fourier transform (DTFT) is defined by

$$F(\omega)=\sum_{n=-\infty}^{\infty}f[n]e^{-jn\omega}\tag{1}$$

Note that this is the common definition, and it is different from the one in your question. First of all, the sequence $f[n]$ doesn't need to be of finite length. Second, the frequency variable $\omega$ is a continuous variable. Your formula can be derived from $(1)$ by assuming that $f[n]$ has a finite length, and by computing the DTFT on a finite set of discrete (angular) frequencies. (Note that the frequency variable $\omega$ in $(1)$ is normalized by the sampling frequency, unlike the discrete frequencies $\omega_n$ in your question.)

Another related transform is the discrete Fourier transform (DFT):

$$\tilde{F}[k]=\sum_{n=0}^{N-1}f[n]e^{-j2\pi kn/N},\qquad 0\le k<N\tag{2}$$

which, for finite length sequences of length $N$, is just a sampled version of the DTFT:

$$\tilde{F}[k]=F\left(\omega_k\right),\qquad \omega_k=\frac{2\pi k}{N}$$

Unlike the DTFT, the DFT can be computed efficiently. Efficient algorithms for computing the DFT are referred to as Fast Fourier transform (FFT) algorithms.

If this doesn't clear up your confusion, leave a comment below.


How does parallel structure implementation of IIR speed up the process ? What happens in case of FIR?

As already pointed out in the comments, it's a misunderstanding that the parallel structure of IIR filters "speeds up the process". The reason for breaking up the transfer function into second-order (and first-order) blocks is to increase numerical stability, especially for fixed-point implementations. The low-order filter sections can be arranged in a cascade structure - i.e., effectively factoring the transfer function -, or in a parallel structure, which corresponds to a partial fraction expansion of the transfer function.

Note that if we don't take quantization effects into account, the direct implementation of the total transfer function, the cascade structure, and the parallel structure give exactly the same output signal. There is no difference in delay, as you seem to believe. Even though in the direct implementation you have a longer delay line, you don't need to wait till the input signal has traveled through the delay line, because the delay line is tapped. So the first non-zero input sample can directly produce a non-zero output sample.

Let's now come to your question about FIR filters. The transfer function of a causal $N^{th}$-order FIR filter with impulse response $h[n]$, $n=0,1,\ldots,N$ is given by

$$\begin{align}H(z)&=\frac{h[0]z^N+h[1]z^{N-1}+\ldots +h[N]}{z^N}\\&=h[0]+h[1]z^{-1}+\ldots +h[N]z^{-N}\end{align}\tag{1}$$

where the second line of Eq. $(1)$ is (trivially) a partial fraction expansion. That means that the "parallel structure" of an FIR filter simply equals a direct form implementation. In other words, there is no distinct parallel structure implementation of an FIR filter.


Difference in frequency spectrum of continuous time signal and discrete time signal?

This is the classical sampling theorem. When you sample a continuous signal $x(t)$, you are basically multiplying it by a sample train $s(t) = \sum_{-\infty}^{+\infty}\delta(t-kT_s)$, the value at consecutive $T_s$ being your samples. In frequency domain, it effect is to convolve $X(\omega)$ and $S(\omega)$

$$S(\omega) = \sum_{-\infty}^{+\infty}\delta(\omega-2\pi k/T_s)$$

So the resulting spectrum would be
$$
X_s(\omega) = \sum_{-\infty}^{+\infty}X(\omega-2\pi k/T_s)
$$
There are infinite shifted copies of $X(\omega)$. To begin with your original signal was a single-frequency signal which traced the circumference with frequency $f_m$. When you sampled it at $f_s$ rate, you got infinite copies at $f_m \pm kf_s$, where $k$ is any integer in $(-\infty,+\infty)$ hence it referred to as infinite number of lines after sampling (discrete time).


How to extract the vein pattern from an IR image of finger?

You might want to also try running a ridge filter on the dataset; something like the frangi filter might work.

A sample result is

NB
the white border is due to the border in the image I pulled from this site


It is a very broad question. You need to specify the type of categorization of this recommendation system. Then you can think about which features might be relevant to this system. Generally you need to apply a classification approach on music signals.


How about:

import cv2
import numpy as np

img = cv2.imread('inputs/Lenna-220px.png', 0)

out = []

for k in range(0, 7):
    # create an image for each k bit plane
    plane = np.full((img.shape[0], img.shape[1]), 2 ** k, np.uint8)
    # execute bitwise and operation
    res = cv2.bitwise_and(plane, img)
    # multiply ones (bit plane sliced) with 255 just for better visualization
    x = res * 255
    # append to the output list
    out.append(x)

cv2.imshow("bit plane", np.hstack(out))
cv2.waitKey()


As you have correctly derived, the difference of carrier and instantaneous frequency is

$$
g(s(t)) = \omega(t) - \omega_\mathrm{carrier} = 2\pi k s(t) + 2\pi k s(t)'t
$$

I will omit $t$ for convenience. If $g(s)$ is a linear function of $s$ it must fulfill the following condition:

$$
g(c(a + b)) = cg(a) + cg(b)
$$
So here:
$$
g(c(a+b)) = 2\pi k c(a + b) + 2\pi k \left[\frac{d}{dt}(c(a + b))\right]t\\
= 2\pi k c a + 2\pi k c b + 2\pi k c a' t + 2\pi k c b' t\\
= c(2\pi k a + 2\pi k a' t) + c(2\pi k b + 2\pi k b' t)\\
= cg(a) + cg(b)
$$
Thus, $g(s(t))$ is indeed a linear function of $s(t)$.


Having problems with understanding the meaning of PSD for colored Gaussian noise

Think in terms of noise realizations. When your head is full to the brim of covariances and multivariate Gaussian variables, generate a Gaussian white noise realization of length N, and create, using this sequence, a sequence $y(n) = x(n) - x(n-1)$. Let us call the latter a blue Gaussian noise realization. I plot histograms of these sequencies which show how many values fall into each bin, those bins uniformly partitioning a range of $y, x$ values.

In compliance with your convention, a blue colored graph represents the white noise data, a red colored, the blue noise data. This histograms visually agree with your observation that

$w_n \sim N(0,2σ_n^2)$, 2x higher variance than $x_n$

Each histogram plotted here is an approximate representation of the distribution of noise signal data, id est, the Probability Density Function of noise signal. The realization sequence's length N is 220 samples. The longer sequence length gets and the greater bin size gets, the closer these numerical PDF graphs resemble theoretical curves.
Go on to Power Spectral Density. Here are discrete Fourier transforms (actually, FFTs) of our white-noise and blue-noise realizations, squared and plotted. Color scheme according to your convention (white noise in blue, blue noise in red):

A white-noise PSD (blue graph) is near-constant, a blue-noise PSD graph (red graph) is two symmetrical parabolic segments sewn at 1/4 and 3/4 of the entire range (the shapes are revealing $\omega$ squared function and FFT mirroring). The high-frequency end is in the middle of horizontal axis.
No relationship between variance and PSD is revealed in our calculations and pictures. It appears that these concepts are things that can only be related indirectly, using formulas that have yet to be derived at this stage of learning.
Variance 'measures how far a set of numbers is spread out from their average value' (Wikipedia). Variance is defined in the time/space domain, and this definition never mentions frequencies or wavenumbers or Fourier transform. Although variance can be measured for any "set of numbers", in many signal processing scenarios said set of numbers is almost synonymous with sampled noise values.
A PSD thing, on the contrary, belongs in the frequency domain. The article in Wikipedia states from the beginning

According to Fourier analysis, any physical signal can be decomposed
into a number of discrete frequencies

and continues

PSD of the signal describes the power present in the signal as a
function of frequency. PSD of noises is a specific case of PSD of
signals in general.

So variance and PSD are sort of "orthogonal" concepts; whatever, variance cannot be a function of frequency in any meaningful way, while PSD is.
Compare also the histograms of the white-noise realization PSD (top) and the blue-noise realization PSD (bottom). The white-noise PSD shape replicates, as it should, a normal distribution of source signal values. A peak of the blue-noise PSD shape shows greater share of the high-frequency components in the blue noise signal.

PS There is a number of Guassian Process tutorials on github and elsewhere, mostly Jupyter notebooks. See, for example, https://peterroelants.github.io/posts/gaussian-process-kernel-fitting/.


The power spectral density is NOT the Fourier transform of the variance. For wide-sense-stationary (WSS) signals, the power spectral density is the Fourier transform of the Autocorrelation function, and the variance is the autocorrelation at time offset $\tau = 0$ (when the WSS process is zero-mean).
Consider for example a white noise process (which is done if the OP's randomly generated samples were completely uncorrelated from sample to sample), to a tightly filtered process, both can have the same variance (which is the power if WSS and zero-mean), but VERY different spectrums. A white noise spectrum is constant over all frequencies while the filtered process would be limited to a narrow bandwidth. The relationship between the two is given by Parseval's theorem; the sum of the squares of each frequency component will be the same in both cases.
See this answer from Matt L further detailing variance with relation to the signal power: variance in the time domain versus variance in frequency domain
Update:
Dilip stated this all much more elegantly and concisely in the comments, paraphrased here: The area under the power spectral density function of a WSS process equals the value of the autocorrelation function at offset $\tau=0$. For a zero-mean WSS process, this value is the variance. Therefore deducing the power spectral density from knowledge of the variance is equivalent to asking for a curve given that we only know the area, which we know is not possible.


simulating directional interference to evaluate model performance

You do need to be careful how exactly you define the SNR. I think you are actually talking about SIR (signal to interference ratio). The way you describe it, it sounds like you plan to just compute the power (time averaged energy) emitted by the signal and interference and use that to get after the SIR. This would be fine but, as you mentioned, sometimes the interference is not present. I might think about using a different definition such as instantaneous SIR, since it sounds like sometimes you have signal + interference and sometimes it is just signal. For this reason, I'm not convinced that computing the power as you described is the right thing to do.

I would suggest doing this: take the source audio file and interference audio file, figure out which is longer, zero pad the shorter one so they are of equal length. Now compute the metric, "sample-by-sample SIR", by doing $\frac{s[n]^2}{i[n]^2}$. Just be careful when $i[n]=0$. This is where knowledge of the noise (real noise, not interference) would be helpful so that you could ideally report SINR (signal to interference plus noise ratio) and avoid the denominator being $0$ problem.


Wavelet transform of a spatial convolution

I cannot say I have a clear understanding of this at this time.  However, a few pointers. I'd love to see somebody provide a detailed account. Others bits at: Multiplication in the wavelet domain, what does it look like in real space?

The nonexistence of a wavelet function admitting a wavelet convoluton theorem of the Fourier type, 1994, A. R. Lindsey, unpublished report
Convolution using the undecimated discrete wavelet transform, 1996

  Convolution is one of the most widely used digital signal processing
  operations. It can be implemented using the fast Fourier transform
  (FFT), with a computationalcomplexity of $O(N \log N)$. The
  undecimated discrete wavelet transform (UDWT) is linear and shift
  invariant, so it can also be used to implement convolution. In this
  paper, we propose a scheme to implement the convolution using the
  UDWT, and study its advantages and limitations.

Convolution theorems for linear transforms, 1998

  This correspondence explores the existence of convolution theorem for
  linear transformations under a variety of different assumptions. There
  are eight convolution theorems, all Fourier-related with only N
  operations in the transform domain and no ordering constraints on the
  convolution components in the result. They include circular
  convolutions and correlations.

The convolution theorem for the continuous wavelet tranform, 2004

  We study the application of the continuous wavelet transform to
  perform signal filtering processes. We first show that the convolution
  and correlation of two wavelet functions satisfy the required
  admissibility and regularity conditions. By using these new wavelet
  functions to analyze both convolutions and correlations, respectively,
  we derive convolution and correlation theorems for the continuous
  wavelet transform and show them to be similar to that of other joint
  spatial/spatial–frequency or time/frequency representations. We then
  investigate the effect of multiplying the continuous wavelet transform
  of a given signal by a related transfer function and show how to
  perform spatially variant filtering operations in the wavelet domain.
  Finally, we present numerical examples showing the usefulness of
  applying the convolution theorem for the continuous wavelet transform
  to perform signal restoration in the presence of additive noise.

On the uniqueness of the convolution theorem for the Fourier transform, 2008

  This paper shows that members of the fourier transform family are the
  only linear transforms that have a convolution theorem, that is, that
  can replace $O(N^2)$ operations of a convolution in a time domain by
  $O(N)$ operations in a transform domain. Generally, there is an
  additional cost to compute the transform itself. Our observation is
  motivated by recent activity in wavelet and subband decompositions and
  related spectral analyses, which are attractive alternatives for
  signal compression applications. A natural question when using such
  techniques is to determine if convolutions of $N$-point signals can be
  calculated with fewer operations in a compressed transform domain than
  in an uncompressed time domain. The answer is negative for a broad set
  of assumptions. This paper indicates what assumptions must be relaxed
  in seeking a linear transform that has a convolution theorem
  comparable to the convolution theorem for Fourier transforms.


You have to define numerator and denominator polynomials. In your case you have

b = 1;
a = [1,-0.2, -0.1];


Essential Matrix from line correspondences

When I think of lines in correspondence, my mind immediately goes for Pluecker representations, which are very suitable for such estimation tasks. In fact this work exists for camera pose estimation from lines and might be a good starting point:

  Přibyl, Bronislav, Pavel Zemčík, and Martin Čadík. "Camera Pose Estimation from Lines using Pl\" ucker Coordinates." arXiv preprint arXiv:1608.02824 (2016).
  (link to the project page)

As you mentioned, there are various resources, especially when your camera is calibrated (I assume this because you ask for the essential matrix). If you really like to have a quick and dirty solution, you could either use the endpoints of the line segments (if you have them) or the intersection points of the lines.

In fact this just reminded me of one nice work, which generates virtual points on cameras along the lines and uses them for pose computation:

  Line-Sweep: Cross-Ratio for Wide-Baseline Matching and 3D Reconstruction
  S. Ramalingam, M. Antunes, D. Snow, G. H. Lee, S. Pillai
  Computer Vision and Pattern Recognition (CVPR), 2015
  (link to the project page)

Let's keep here posted with updates.


How to find SNR of received DSB:AM signal with white noise after demodulation & LPF'ed?

Since you don't want a full answer, I'll just provide a couple of hints. Your calculation for the noise power is correct. Regarding your questions:

The power of a signal is the integral of its PSD.
You're looking for the Wiener-Khinchin theorem, and especially its application to LTI systems.


Limiter, compressor, AGC, dynamics processor?

one category has a protection property and has a very fast response to large amplitude transients.

another category attempts to maximize sensitivity.

another group tries to put the receiver in a predictable “sweet spot” or operating point.

some can be “slow” and others are “fast”.

none of these are mutually exclusive.

advances in converter bit resolution (dynamic range) moves these techniques along.

One might consider schemes that adapt to be related.


How to plot the poles and zeros of this Bandpass Filter in MATLAB?

You've designed a second-order band pass filter, so the number of terms in the denominator is just fine (it's a second-order polynomial). Your design is OK, apart from the scaling, i.e., the maximum of the magnitude of the transfer function is not equal to $1$, which might be desirable.

Your transfer function has the form

$$H(z)=\frac{z^2-1}{(z-p)(z-p^*)}=\frac{z^2-1}{z^2-2zr\cos(\omega_p)+r^2}\tag{1}$$

where $p$ is a complex-valued pole:

$$p=re^{j\omega_p}\tag{2}$$

with radius $r$ and pole angle $\omega_p$. Note that in general the pole angle $\omega_p$ and the peak frequency $\omega_0$ (where the magnitude of the frequency response $H(e^{j\omega})$ achieves its maximum) are not equal. However, in your case they are almost identical because $\omega_p$ is relatively close to $\pi/2$. See this answer for more details on this matter.

If you want the frequency response peak value to be equal to $1$, you need to scale your transfer function by the factor $(1-r^2)/2$ (also explained in the answer quoted above).

The pole angle $\omega_p$ depends on the desired center frequency, so it can't always be equal to $\pi/2$ as you seem to have thought (even though in your calculations you chose it according to the desired center frequency).

The Matlab function freqz computes the complex frequency response of a discrete-time filter.
For plotting the poles and zeros you can use the function zplane.


What are the parameters of a calibrated spherical camera?

...if we have a pinhole camera model several parameters describe the specific camera (such as aspect-ratio, focal length, principal point, distortion parameters etc).

"Distortion parameters" does not sound like a typical pinhole camera model. A pinhole camera does not have a "finite apperture" or way of focusing light other than a tiny little opening which all rays of light (that will go on to produce the image) go through. You can take lens distortions into account by adding them to the pinhole model as additional operations but that then is not the standard pinhole camera model.

The pinhole camera model is a very simple transformation of a point from a global 3D coordinate system to a local 2D coordinate system. It does nothing more than that. It maps the angle at which a point appears within its field of view to a point on its focusing plane. Another useful way to be thinking about lenses is in terms of mapping points in front of them to points behind them.

The pinhole camera model is the simplest of those transformations and requires:

The position of the camera in global coordinates (Where is it?).
The orientation of the camera (Where does it point to?)
The "Focal Length"
The sensor size
The principal point

The focal length in particular is exactly the "thing" that creates this simple map. The shorter the focal length, the wider is the field of view within which points get mapped to a plane behind the pinhole and vice versa.

But notice here, rays in front of the pinhole enter it at some angle and leave exactly behind the pihnole at the same angle.

  I was wondering however if I replaced the pinhole camera with a spherical camera (such as 360 camera for example) what parameters describe such camera?

A pinhole 360 camera is a pinhole camera with 360 degrees field of view. Exactly the same set of parameters but now practically points can enter the pinhole from anywhere around it and depart at the same angle behind it. Notice here that in this view the "focal length" is still valid as a "ratio" of lengths.

But in reality, there are no 360 degree "pinholes". Instead, the wide field of view is acquired in some other way. That might be a mirror for example or multiple cameras with various lenses.

In this case, it is even more useful to be thinking in terms of rays entering the lens and rays departing from the  lens (or optical system in general). The general model in this case is established around a transfer matrix that maps "entry points" to "exit points" and basically can model any sort of configuration you like.

In terms of references, Geyer and Daniilidis:A unifying theory for central panoramic
systems and practical implications is probably still relevant and since then, Daniilidis has also put together a volume of relevant work, the contents of which you might find useful, even if it is just for more relevant terminology to get you closer to what you are dealing with.

Hope this helps.


Can I use frequency hopping spread spectrum technique for audio steganography?

Yes, it is possible to hide data inside audio files using FSSS technique. Have a look at this publication and this one.

Since HAS (Human Auditory System) is more sensitive than HVS (Human Visual System), you are usually limited with audio steganography comparing to image steganography. There are not many detalied publications with "Audio Steganography" tag. Interests about hiding digital data inside audio files have tended into watermarking field. So I strongly suggest you to research whatever you are looking for as "Audio Watermarking". You will find more contents related to watermarking.

Few more words about Audio Steganography:

LSB technique provides you more capacity to hide inside audio files, yet it will be very fragile to resist against any attacks such as mp3 compression. According to my experiences, DSSS and Echo Hiding methods are more robust comparing to LSB Coding and Parity Coding. DSSS provides a bit more capacity than Echo Hiding, but the noise is audible comparing to Echo Hiding.


Normalized autocorrelation of a sum of two signals

That would be a scaled autocorrelation such that its maximum value, hopefully attained at $\tau=0$, has value $1$.


Investigating correlation for unequal signal lengths

If you are using MATLAB, there is a built in object type, called

  timeseries

, that is very useful for this sort of problem.

>    Methods for class timeseries:
>
> addevent              getabstime            gettsbeforeatevent    mean
> pvset                 times                  addsample
> getdatasamples        gettsbeforeevent      median
> rdivide               timeseries             append
> getdatasamplesize     gettsbetweenevents    min
> resample              tsprops                createTstoolNode
> getinterpmethod       hasduplicatetimes     minus                 set
> utArithCommonData      delevent              getprop
> idealfilter           mldivide              setabstime
> utArithCommonOutput    delsample             getqualitydesc
> init                  mode                  setinterpmethod
> utArithCommonTime      detrend               getsamples            iqr
> mrdivide              setprop               utChkforSlashInName
> display               getsampleusingtime    isequal
> mtimes                setuniformtime        utGetEventTime         eq
> gettsafteratevent     isequalwithequalnans  plot                  std
> utStatCalculation      filter                gettsafterevent
> ldivide               plus                  sum                   var
> get                   gettsatevent          max
> pvget                 synchronize
>
> Static methods:
>
> loadobj               tsChkTime             tsgetrelativetime

The setuniformtime method will do the interpolation for you. The other methods are often useful.


Confusion regarding wav to MFCC conversion

There are many variations in the ways MFCCs can be extracted. The most common ones are:

STFT parameters like frame rate, FFT size, FFT window.
Filter bank parameters like number of bands, whether filter normalization is used, and whether the filters are applied to the magnitude or power spectra.
Subtle variants in the definition of the discrete cosine transform.

And it is actually not correct to say that there are 13 MFCC - it is indeed very common to keep the 13 first coefficients, but you could very well use a larger cosine transform matrix and get more coefficients. The farther you go, the less information the coefficients will capture (since they will correspond to dimensions of decreasing variance - assuming the discrete cosine transform is a good approximation of a PCA on the mel coefficients data, which it is).

Check this to learn more about all the subtle differences between several "famous" MFCC extraction toolboxes.

In the end, such or such combination of MFCC extraction parameters might give you a difference of at most one or two percentage points at whatever task you are using the MFCC for (for example speech recognition), but there are more interesting things to try to tweak and improve than these parameters.


The tricky part tends to be time shifting the state variables. You try to avoid spending extra cycles on explicitly doing this but ideally this can be all kept in registers or just happens during the natural flow of the instructions. The best way to do this will highly depend on your instruction set and/or specific capabilities of the hardware.

In general Direct Form II is a poor choice. The transfer function between the input and the state variables is vine solely through the poles which leads frequently to stability and noise problems. Transposed Form II is almost always the better choice. Only exception would be a fixed point processor with a double precision accumulator in which Direct Form I would be the preferred choice.


Backward finite difference differentiation filter frequency response

In discrete-time signal processing terminology a backward (one sample) difference filter is
$$y[n] = x[n] - x[n-1]$$
which is a causal, LTI system with an FIR impulse response of
$$h[n] = \delta[n] -\delta[n-1]$$
The associated frequency response is then the DTFT of the impulse response as
$$\begin{align}
H(e^{j\omega}) = \mathscr{F}\{ h[n] \} &= 1 - e^{-j\omega} \\
&= e^{-j\omega/2}(e^{-j\omega/2} - e^{-j\omega/2}) \\
&= 2j e^{-j\omega/2} \ \tfrac1{2j}(e^{-j\omega/2} - e^{-j\omega/2}) \\
&= e^{-j\omega/2} 2j  \ \sin(\omega/2) \\
&\approx e^{-j\omega/2} j  \ \omega \qquad \text{for } |\omega| \ll \pi \\
\end{align}$$
So you have a delay of $\frac12$ sample, an additional phase shift of +90° (due to the $j$ factor), and a gain proportional to frequency (at least for low frequencies, relative to Nyquist).  Other than the half-sample delay, this is what you expect from a differentiator.
If you happen to use an $N$-sample backward difference such as
$$y[n] = x[n] - x[n-N]$$
then the associated frequency response will be
$$H_N(e^{j\omega}) = 1 - e^{-jN\omega} $$
This is similar to the first orde case, but the delay is $\frac{N}2$ samples, and the gain is $N \omega / 2$ .


Improve contrast when image has a flood of white balancing

Histogram equalization is a global operation which can result in some areas being adversely effected at the cost of the rest of the image looking better. Consider using the Retinex algorithm. It works very well at improving both local and global contrast enhancement. There's a flurry of papers on this topic and the algorithm is somewhat straightforward to implement. The most challenging part is the selection of parameters.

Here's a list of papers and here's one that's very useful.


Which is the most suitable programming language for image processing in scientific researches?

Your question as is has a very simple answer: "the best language is the one I use". Pros and Cons (P&C) would give some more value to it, as commented by Peter and Marcus. So, back to the topic, but concerning with P&C, here is my $\epsilon$...

Your application does not seem to be restrict only to simulation, but also to web interface; moreover, you want a roadmap to it, probably in a asynchronous way dealing with branches, and therefore with some version control.

Matlab can be a choice if your audience has all toolboxes you have for this work (web + image + signal + ... ). Does it?

Forget control versioning if you intend to use Simulink.

I would make thumbs up for Python as well. You have support of a huge community. The scientific packages are numerous. It is a complete language, from web applications to visual to scripting to ... And it is free!


Matlab applying an FIR filter to a signal, generated in Matlab

In general your result match expectations, I don't see anything in there that's unusual or unexpected.

Filter design is a complicated trade off between complexity/cost, stopband attenuation, pass band ripple, transition steepness, phase distortion, time domain ringing, non-causality, latency, etc. The best way to go about is to clearly articulate the requirements for your specific application and then evaluate any design against these requirements line-by-line.

I'd also be careful with spectral analysis: direct FFT works sometimes, but something like pwelch() is typically much more robust.


Understanding where the constant $2/N$ comes from in Fourier transformation

The two comes from the fact that a real valued sinusoid is really the average of two complex one.

$$ \cos( \theta ) = \frac{ e^{i \theta} + e^{-i \theta} }{2}  $$

There is the "2" in the denominator.  In the DFT, the other half is located at bin $N-k$, if $k$ is the bin number.

That definition of cosine comes straight from Euler's equation:

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

The $N$ comes straight from the definition of the DFT.  The relationship beteen the amplitude, N, and the magnitude of the DFT bin is why I like to use the $1/N$ normalized form of the the DFT over the more conventional unnormalized form.

I recommend that you read my first blog article The Exponential Nature of the Complex Unit Circle for an understanding of Euler's equation and my blog article DFT Bin Value Formulas for Pure Real Tones for understanding the bin values of a pure tone and what "leakage" really means.  Equation (19) is the answer to your "big question".

The $\Omega$ comes from the continuous FT.  The DFT is the Discrete Fourier Transform.  Too often, they are conflated leaving confusion like yours.  It takes a lot of heavy math to understand the relationship between the two.  The premise of my blog articles is that the DFT can be learned and understood straight from the summation definition without any reference to the continuous case.

Hope this helps.

Ced


How does 3D Time-of-Flight (TOF) camera work?

The (unwrapped) phase between two constant same frequency signals is linearly related to any delay between the two signals.

Measuring phase of a fixed frequency periodic signal over many periods is sometimes more accurate than trying to detect delay locally (by zero crossing, etc.)


Any experiences for plotting a stationary wavelet transform?

I don't do Python, I'm an old person sticking to his old Matlab (codes and) habits. However, up to extension/wavelet/border issues,  SWT is a discrete equivalent to CWT. And in most versions, the number of samples in approximations or details is the same (which is not the case for the DWT). Hence, you can concatenate 1D rows of details (or approximations) into 2D images, akin to traditional scalograms. The following images and code show the process. I have generated a random piecewise polynomial signal with increasing degrees. Then rows of details, and an image of the rows concatenated. The motivation is to address the impact of vanishing moments on piece-wise polynomials, with border effects.
Here are two different realizations.

The Matlab code, that may be reproduced in Python:

dataLength = 512;
data = zeros(dataLength,1);
nChunk = 4;
lChunk = dataLength/nChunk;
for iChunk = 0:nChunk-1
    idxChunk = iChunk*lChunk+(1:lChunk)';
    polyChunk = rand(iChunk+1,1)-0.5;
    dataChunk = polyval(polyChunk,linspace(-0.5,0.5,lChunk));
    data(idxChunk) = dataChunk;
end

nLevel = 4;
[swa,swd] = swt(data,nLevel,'db2');

figure(1);clf
subplot(1,3,1)
plot(data,'x');axis tight
for iPlot = 1:nLevel
    subplot(nLevel,3,3*(iPlot-1)+2)
    plot(swd(iPlot,:));;axis tight
end
subplot(1,3,3)
imagesc(swd)

Other related questions/answers:

STFT and DWT (Wavelets)
Does SWT/ISWT require intermediate approximation coefficients to represent/reconstruct the original signal?
Extract approximation and detail coefficients
Undecimated DWT vs. CWT. In what cases is one preferable over another?
Additive white Gaussian noise and undecimated DWT
Reproducing paper results about a wavelet transformation using python


Kalman Filter Estimate vs ACF Least Squares Estimate

As noted on p. 11 of this site: http://www.pages.drexel.edu/~pyo22/mem640/lecture04-Estimation/mem640Lecture-Estimation.pdf  least squares is a "special" case of Kalman Filtering. Additional details can also be found on other web sites, such as: http://www.insidegnss.com/node/3445

Also, as noted in both links above - although least squares is a special case of Kalman Filtering - the two algorithms differ from each other.  Note that the 'weighting matrix' W reduces to the identity matrix in the least squares case.  So you might expect that the results would also differ when you process the data.

As to which one will be more accurate - it all depends on how closely the 'signal processing model' of your algorithm matches the process that generates your input data.  I've programmed Kalman Filters whose results were completely wrong because: 1) it started with a bad first guess, or 2) the state model being used was wrong, or 3) the input data was far too noisy.  I've also programmed 'hypothesis testing' arrays of Kalman Filters, where all filters operated on the same data, but each filter was set up to match a particular target motion (ie: zig right, zig left, go straight, etc.) - in that case, only the filter that closely predicted the actual input data was the correct one (ie: residuals were very small), and the residuals of all the others were way out of range.


First, on distinguishing low contrast features from noise, then on a possible alternate course of action to solve what I think you are trying to do.

From my experience, distinguishing between noise and low contrast features is a hard thing to do. The problem arises quite often when you try to do edge detection on an image. Using the most common algorithms such as Canny, Sobel, LoG, etc. generally results in a good deal of noise, which is dependant of the parameters you set (e.g.: a 'threshold' needs to be set (chosen) in order to define at which intensity value an edge-candidate pixel is classified as part of the edge or not, so what if the noise is more 'intense' than the correct edge?)

So, what's the point with stereo correspondence/disparity? Some scholars tackled on the problem of noisy edge detection and there is right now a pretty good edge detector, namely the Berkeley edge detector, which use machine learning techniques to 'improve' the edge detection.

Here are 3 images, left is original, centre is with Canny edge and right is with Berkeley. (Images credit: ETHZ-Shape-Class Dataset)

You can see that (while maybe some tuning could have been done on Canny) Berkeley is providing a 'better' (less noisy) edgemap than Canny.

So, you could see that there are techniques that you could possibly try to adapt in order to tell the salient features from the noise. I can't guaranty it will work since I haven't done it myself. Also, please note that these machine learning algorithms are pretty slow which prevent their use in realtime processing of images.

This said, what you really want to do is measure the dissimilarity between the to images by mean of correlation techniques, which was not successful apparently. Let me tell you about what is known as keypoint-descriptors. You could read on this topic there. The broad idea is to be able to evaluate the correspondence of some (a lot of) features between two images. It is often used in object detection (where you want to find the presence and location of a certain object in an input image) or in image stitching (such as panoramas maker, where you want to 'align' two or more images). Now, how can you tell that you have two images corresponding to the same object? These methods using keypoint-descriptors are, in the end, giving a measure of similarity between both images. Actually, lots of them are measuring the dissimilarity between the features to tell if they are similar or not...

I hope you will be able to better understand the situation with this. This is not a direct pre made solution, but it could help you in your work!


I don't know how WCDMA does it, but CDMA2000 encodes the data twice or three times (it's been a while): once with a "long code" (42 or 43 bits, so it takes days to cycle through even at chip rates in the MHz), once with a "short code" (16 bits, with zero-stuffing so that the cycle is exactly 65536 chips long, and the spectrum is flat), and IIRC once with a 64-bit long Walsh code.  But There's a channel that's only encoded with the long code; that one is used by the base station to assign a specific set of codes and timings to the phone.

So, basically, it's done dynamically.

The CDMA2000 specs are out there.  When I looked (over ten years ago, which is why this answer is a bit fuzzy) the "real" spec was behind a paywall, but there were draft specifications that were pretty close.


Sound card == audio interface

The only difference is that the name "audio interface" is usually used for pro audio equipment. On the other hand, "sound card" usually refers to consumer audio.

Currently there are no cheap multi-input audio interfaces. You may want to check out second hand market, but even there USB and 1394 interfaces are expensive.

The only cheap solution that comes to my mind is to use multiple stereo-input interfaces and synchronize theirs clocks using SPDIF. Then you can link them on Linux (in asoundrc) or Mac OS. I don't know if it's possible on Windows.

However, I have never tried synchronizing clocks using SPDIF myself.

Without synchronizing clocks you need to use resampling which is computationally expensive and can degrade quality. alsa_in and alsa_out tools for Linux/JACK can do this. But they're usable only if you don't need accurate timing as they don't do latency compensation.


Intuition behind commutativity of convolution in LTI systems

In a discrete-time system such as the one that you have, the number $y[n_0]$ (here $n_0$ is a fixed integer) is a sum of the form $$\sum_{k=-\infty}^\infty h[k]x[n_0-k]$$ which can be re-arranged via a change of variables (replace $k$ by $n_0-\ell$) to $$\sum_{\ell=-\infty}^\infty h[n_0-\ell]x[\ell].$$ So, the commutativity of the convolution is trivial.  The issue is the interpretation that you put on it.  As Laurent Duval/s answer points out, the systems A and B are not equivalent in any sense of the term.  If the signal $x$ were replaced by a different signal $\hat{x}$, then system A would have output $\hat{y} = h \star \hat{x}$, but you wouldn't get the same output $\hat{y}$ if system B were excited by $h$; the impulse response of system B continues to be $x$, and system B thus has output $x \star \hat{x} =  \hat{x}\star x \neq h \star \hat{x}$.


Is a leaky integrator the same thing as a low pass filter?

A so-called leaky integrator is a first-order filter with feedback. Let's find its transfer function, assuming that the input is $x(t)$ and the output $y(t)$:

$$
\frac{dy(t)}{dt} + Ay(t) = x(t)
$$

$$
\mathcal{L}\left\{\frac{dy(t)}{dt} + Ay(t)\right\} = \mathcal{L}\left\{x(t)\right\}
$$

where $\mathcal{L}$ denotes application of the Laplace transform. Moving forward:

$$
sY(s) + AY(s) = X(s)
$$

$$
H(s) = \frac{Y(s)}{X(s)} = \frac{1}{s + A}
$$

(taking advantage of the Laplace transform's property that $\frac{dy(t)}{dt} \Leftrightarrow sY(s)$, assuming that $y(0) = 0$).

This system, with transfer function $H(s)$, has a single pole at $s = -A$. Remember that its frequency response at frequency $\omega$ can be found by letting $s=j\omega$:

$$
H(j\omega) = \frac{1}{j\omega + A}
$$

To get a rough view of this response, first let $\omega \to 0$:

$$
\lim_{\omega \to 0} H(\omega) = \frac{1}{A}
$$

So the system's DC gain is inversely proportional to the feedback factor $A$. Next, let $w \to \infty$:

$$
\lim_{\omega \to \infty} H(\omega) = 0
$$

The system's frequency response therefore goes to zero for high frequencies. This follows the rough prototype of a lowpass filter. To answer your other question with respect to its time constant, it's worth checking out the system's time-domain response. Its impulse response can be found by inverse-transforming the transfer function:

$$
H(s) = \frac{1}{s+A} \Leftrightarrow e^{-At}u(t) = h(t)
$$

where $u(t)$ is the Heaviside step function. This is a very common transform that can often be found in tables of Laplace transforms. This impulse response is an exponential decay function, which is usually written in the following format:

$$
h(t) = e^{-\frac{t}{\tau}}u(t)
$$

where $\tau$ is defined to be the function's time constant. So, in your example, the system's time constant is $\tau = \frac{1}{A}$.


Why does Phase-locking value show spurious effects at the beginning & end of a signal?

What was previously posted as an answer should have been an extended comment -- sorry for that. There simply was not enough space for what I thought I'd write about how the Octave code seemed to behave. Now that the code is cleared, what should have been the answer in the first place (given the rather obvious signs), is now the answer: spectral leakage.
The Hilbert transform is calculated by zeroing half the FFT of the signal and then recomputing the analytical signal through the IFFT. Brute-force zeroing will always lead to some leakage due to the nature of the limited length of the data. This can be easily tested with your code. s1 is a simple sin(), so its Hilbert transform will be a -cos(). Simply plot the difference between -cos()-hx.imag, and you'll see the leakage. The same goes for the second signal.
This is the changed code in Python:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
f = 10
t = np.linspace(0, 100/f, 1000)
N = 50f = 10
nt = 1000
t = np.linspace(0, 100/f, nt)
N = 50
hor = np.ones(nt)
ver = np.ones(N)
r = np.array( [ np.random.uniform(low=0, high=2*np.pi) for i in range(N) ] )
s1 = np.array( [ np.sin( 2*np.pi*f*t ) for i in range(N) ] )
s2 = np.array( np.sin( 2*np.pi*( f*ver[:,None]*t + hor*r[:,None] ) ) )
c1 = -np.array( [ np.cos(2*np.pi*f*t) for i in range(N) ] )
c2 = -np.array( np.cos( 2*np.pi*( f*ver[:,None]*t + hor*r[:,None] ) ) )
hx = signal.hilbert(s1)
hy = signal.hilbert(s2)
nx = hx/np.absolute(hx)
ny = hy/np.absolute(hy)
p = np.absolute( np.mean( nx*ny.conj(), axis=0 ) )
plt.plot(p)
plt.show()

and this is the test (c2 will show non-overlapping plots due to the randomness):
plt.plot(c1.transpose() - hx.imag.transpose())
plt.show()

The Octave code, for checking:
f = 10
t = linspace(0, 100/f, 1000);
N = 50;
c = ones(N, 1);
ft = f*t.*c;
s1 = sin(2*pi*ft);
c1 = -cos(2*pi*ft);
r = rand(N,1).*c;
s2 = sin(2*pi*(ft + r));
c2 = -cos(2*pi*(ft + r));
hx = hilbert(s1');
hy = hilbert(s2');
nx = hx./abs(hx);
ny = hy./abs(hy);
p = abs(mean((nx.*conj(ny))'));
grid on
subplot(2, 1, 1)
plot(p)
subplot(2, 1, 2)
plot(c1' - imag(hx))


What does the "smoothing optical flow" mean?

Optical ﬂow is an approximation of the local image motion based upon local derivatives in a given sequence of images. That is, it speciﬁes how much each image pixel moves between adjacent images. (The same concept extends on 3D as well).

The concept of Optical flow derives from the fact that the moving patterns cause temporal varieties of the image brightness. It is assumed that all temporal intensity changes are due to motion only. This is usually true but there are many exceptions (see below). These variations are measured using derivatives. The optical flow methods try to calculate the motion between two image frames which are taken at times $t$ and $\delta t$ at every pixel position. The following constraint equation can be used to derive the optical flow - which is called as Image constraint equation or optical flow equation:

$$
{{\partial I} \over {\partial x} }\Delta x + {{\partial I} \over {\partial y}} \Delta y + {{ \partial I} \over {\partial t}} \Delta t = 0
$$

Refer to wiki page, or this lecture note for more details.

Thus the computation of diﬀerential optical ﬂow is, essentially, a two-step procedure:

measure the spatio-temporal intensity derivatives (which is equivalent to measuring the velocities normal to the local intensity structures) and
integrate normal velocities into full velocities, for example, either locally via a least squares calculation or globally via a regularization.

However, the above equation alone is not sufficient for computation since it involves two variables with only one equation/constraint. BKP Horn, introduced another constraint in their seminal paper.

BKP Horn, Brian Schunck "Determining  Optical  Flow" Artificial Intelligence, Vol. 16, No. 1–3, August 1981, pp. 185–203.

Accordingly, in a rigid body motion (with non deformation assumption), neighboring pixels have similar velocity. Based on this they derived another constraint based on velocity:

$$
\nabla ^2 u = {{\partial ^2 u} \over {\partial x^2}} + { {\partial ^2 u} \over {\partial y^2}} = 0
$$

and
$$
\nabla ^2 v = {{\partial ^2 v} \over {\partial x^2}} + { {\partial ^2 v} \over {\partial y^2}} = 0
$$
where
$$
u = {{\delta x} \over {\delta t}} \text{ and  } v = {{\delta y} \over {\delta t}}
$$

The above two equations are considered as Smoothing Constraints of optical flow.
Read this: cv-online page and another reference here.


For linearity it is irrelevant if the output depends on future values. This only affects the system's causality. For linearity you just need to check whether the response to the signal

$$ax_1[n] + bx_2[n]$$

equals

$$ay_1[n] + by_2[n]$$

where $y_1[n]$ and $y_2[n]$ are the individual responses to inputs $x_1[n]$ and $x_2[n]$, respectively.

And of course you can write the input/output relation as

$$y[n+1] = -\frac12 x[n+2] - y[n]\tag{1}$$

I will now characterize all sequences that satisfy the recursion (1) for the specific input signal $x[n]=\delta[n]$. Applying the $\mathcal{Z}$-transform to (1) we get (with $X(z)=1$)

$$Y(z)=-\frac12\frac{z^2}{1+z}\tag{2}$$

The function $Y(z)$ corresponds to one left-sided sequence and to one right-sided sequence. Note that left-sided and right-sided doesn't necessarily mean anti-causal and causal, but it only means that $y(n)=0$ for $n>N$ or for $n<N$, respectively, for some positive or negative integer $N$. These two sequences are given by

$$y_{p1}[n]=\frac12 (-1)^{n}u(n+1)\tag{3}$$
and
$$y_{p2}[n]=-\frac12 (-1)^n u(-n-2)\tag{4}$$

Note that neither (3) nor (4) is causal. The sequences $y_{p1}[n]$ and $y_{p2}[n]$ are particular solutions of the difference equation (1) for $x[n]=\delta[n]$. They are in fact impulse responses of two different linear time-invariant (LTI) systems described by the difference equation (1). The impulse response $y_{p1}[n]$ was mentioned by Dilip Sarwate in a comment below, and something similar to $y_{p2}[n]$ was mentioned by Robert B.-J. I completely agree with them that the difference equation (1) describes an LTI system, if it is used to describe a system with either (3) or (4) as impulse response.

However, the sequences (3) and (4) are not the only sequences that satisfy (1). Looking at the homogeneous difference equation

$$y[n+1]=-y[n]\tag{5}$$

we see that

$$y_h[n]=c(-1)^n,\quad c\in\mathbf{R}\tag{6}$$

is a solution of (5). Consequently, the general solution of (1) (with $x[n]=\delta[n]$) can be written as

$$y[n]=y_{p1}[n]+y_h[n]=\frac12 (-1)^{n}u(n+1) + c(-1)^n\tag{7}$$

Of course we could as well have used $y_{p2}[n]$ to define the general solution. Note that $y_{p2}[n]$ is obtained from (7) by choosing $c=-1/2$.

Now assume that we scale the input signal and use $x[n]=2\delta[n]$ instead. Obviously, the solutions $y_{p1}[n]$ and $y_{p2}[n]$ will change accordingly, but the general solution (7) will not (unless $c=0$ or $c=-1/2$), because the homogeneous solution (6) does not depend on $x[n]$. This is why I claim that the difference equation (1) does not necessarily describe an LTI system. In our case, it describes exactly two LTI systems and an infinite number of systems which are neither linear nor time-invariant.


finding phase and group delay from frequency response

Define a frequency response

$$\tilde{H}(e^{j\omega})=\frac{e^{-j\omega}-|\alpha|}{1-|\alpha|e^{-j\omega}}\tag{1}$$

and note that the group delay $\tilde{\tau}(\omega)$ corresponding to $(1)$ is related to the group delay $\tau(\omega)$ of the original frequency response $H(e^{j\omega})$ by

$$\tau(\omega)=\tilde{\tau}(\omega-\theta)\tag{2}$$

because $(1)$ is just a scaled and frequency-shifted version of the original $H(e^{j\omega})$.

$\tilde{H}(e^{j\omega})$ is the frequency response of a first-order all-pass filter with a real-valued pole at $z_{\infty}=|\alpha|<1$, and a real-valued zero at $z_0=1/|\alpha|$. The original frequency response $H(e^{j\omega})$ describes a first-order all-pass filter with a complex-valued pole at $z_{\infty}=\alpha=|\alpha|e^{j\theta}$ ($|\alpha|<1$), and a zero at $z_0=1/\alpha^*$.

Now write

$$\tilde{H}(e^{j\omega})=\frac{B(e^{j\omega})}{A(e^{j\omega})}\tag{3}$$

and note that

$$\tilde{\phi}(\omega)=\phi_B(\omega)-\phi_A(\omega)\tag{4}$$

where $\tilde{\phi}(\omega)$ is the phase of $\tilde{H}(e^{j\omega})$, and $\phi_A(\omega)$ and $\phi_B(\omega)$ are the phases of $A(e^{j\omega})$ and $B(e^{j\omega})$, respectively. Consequently,

$$\tilde{\tau}(\omega)=\tau_B(\omega)-\tau_A(\omega)\tag{5}$$

where $\tau_A(\omega)$ and $\tau_B(\omega)$ are the group delays of $A(e^{j\omega})$ and $B(e^{j\omega})$, respectively.

We have

$$\phi_B(\omega)=-\arctan\left(\frac{\sin(\omega)}{\cos(\omega)-|\alpha|}\right)\tag{6}$$

and

$$\phi_A(\omega)=\arctan\left(\frac{|\alpha|\sin(\omega)}{1-|\alpha|\cos(\omega)}\right)\tag{7}$$

where I've ignored any ambiguities by $\pm\pi$ due to the principal value of $\arctan()$ because they are irrelevant for computing the group delay.

Taking the negative derivative of $(6)$ and $(7)$ with respect to $\omega$ gives

$$\tau_B(\omega)=\frac{1-|\alpha|\cos(\omega)}{1-2|\alpha|\cos(\omega)+|\alpha|^2}\tag{8}$$

and

$$\tau_A(\omega)=\frac{|\alpha|^2-|\alpha|\cos(\omega)}{1-2|\alpha|\cos(\omega)+|\alpha|^2}\tag{9}$$

From $(5)$ we get

$$\tilde{\tau}(\omega)=\frac{1-|\alpha|^2}{1-2|\alpha|\cos(\omega)+|\alpha|^2}=\frac{1-|\alpha|^2}{\big|1-|\alpha|e^{-j\omega}\big|^2}\tag{10}$$

From $(2)$ the group delay of the complex first-order all-pass is given by

$$\tau(\omega)=\frac{1-|\alpha|^2}{1-2|\alpha|\cos(\omega-\theta)+|\alpha|^2}=\frac{1-|\alpha|^2}{\big|1-|\alpha|e^{-j(\omega-\theta)}\big|^2}\tag{11}$$

The difference with the given solution in Schaum's outline is the fact that in $(11)$ you have $|\alpha|$ instead of the complex-valued $\alpha$. Note that $(11)$ is correct and Schaum's formula is wrong because a complex-valued group delay doesn't make much sense.


Is there a recommended amount of attenuation for rectifying edge artefacts?

My guess is that for Spotify, it is not just about avoiding clicks, but also to improve the DJing kind of experience. The fade length seems to be a user setting in Spotify. Based on my testing with 20 Hz and 440 Hz test tones and Virtual Audio Cable, YouTube doesn't fade in or out or do any other kind of click removal when pausing, here with 20 Hz test tone:

The worst-case test signal to be faded in or out would be a full-scale constant-valued signal (for 16-bit audio that would be 32767, 32767, 32767, 32767, 32767, 32767, ...), because it has no audible sound of its own, and because it will still allow to hear the fade as loud as possible. Then we'd need the best possible fade shape, but I don't know what it is, so here is a sub-optimal fade in:
$$f(x) = \begin{cases}0&\text{if } x < 0,\\
\frac{1}{2}-\frac{1}{2}\cos(\pi x)&\text{if } 0 \le x \le 1,\\
1&\text{if } x > 1.\end{cases}\tag{1}$$

The following "floatbeat" program plays this fade twice a second (once in and once out per second):
f = function(x){return 0.5-0.5*cos(x*PI)}, fadeLength=0.03, loopLength=1, sampleRate=44100, p=(t/sampleRate)%loopLength, q=p-loopLength/2, (p < fadeLength? f(p/fadeLength): (q < 0)? 1: (q < fadeLength? 1-f(q/fadeLength): 0))

Try running it and adjust fadeLength. At least for me, using good headphones at normal music listening volume, a fade length 0.03 s is only barely audible as a thump, and increasing it to 0.04 s makes the thump disappear. I think 0.03 s is pretty safe.
But is 0.03 s already so long that it sounds like a fade? Let's apply the fades to a 400 Hz sine tone, in floatbeat:
f = function(x){return 0.5-0.5*cos(x*PI)}, fadeLength=0.03, loopLength=1, freq=400, sampleRate=44100, p=(t/sampleRate)%loopLength, q=p-loopLength/2, (p < fadeLength? f(p/fadeLength): (q < 0)? 1: (q < fadeLength? 1-f(q/fadeLength): 0))*sin(2*PI*t*freq/sampleRate)

Running this, I can't hear a click in fade in or fade out. I know there is a fade from the way it sounds and because I do a lot of music editing, but it does not register as a gradual process, more like a smooth event.
Let's try some other fade shapes of length 0.03 s too, with my comments comparing each to Eq. 1:

Repeating Eq. 1 here for easy clicking.
A linear fade $f(x) = x$ sounds much clickier to me.
Smoothstep $f(x) = 3x^2-2x^3$ sounds maybe tiny bit worse to me.
Asymmetrical polynpomial $f(x) = 6x^2 - 8x^3  + 3x^4$ sounds significantly thumpier.
Smoothstep $f(x) = 10x^3-15x^4+6x^5$ sounds about as good, maybe even a bit better. Less mid-freq but more bassy thump.

Calculating the Fourier transform $F(f)$, where $f$ denotes frequency, of each fade of length 0.03 s appears to confirm what I heard: (For this, I calculated the Fourier transform of the derivative of each fade as function of time $t = \frac{x}{0.03s}$, and compensated for derivation by multiplying the result by a factor $-\frac{i}{f}$.)

I also included the pure tone hearing threshold curve from ISO 226:2003. Some caveats: I do not know the listening volume level so I just shifted the curve in dB scale until it matched my subjective hearing experience with the fade curves peaking above the threshold at roughly those frequencies at which I heard them. Also, the threshold curve has been measured for pure sinusoidal tones whereas here we are dealing with transient sounds. The level comparison between the two cannot be well defined. It would be better to do the analysis with something like a compressive gamma-chirp filter bank which better models hearing. I used closed-back headphones, but still my listening environment has some fan noise and such, so not perfect silence, so there could be some auditory masking going on. Despite all the caveats, I think the threshold curve is useful to show here as it seems to match what I heard. It also hints to that there may be a better fade shape that moves some of the energy from the first spectral side lobe to the main lobe, using window function terminology. Perhaps a better fade shape, one that barely touches the hearing threshold at multiple early lobes, could be obtained by integrating a suitably parameterized Ultraspherical window. Or perhaps there are ways to trace the hearing threshold for all lobes.
I'd choose between 5th order smoothstep and Eq. 1, or 3rd order smoothstep if you want to save some processing power. The choice is not too critical.


Formula to calculate Cepstral coefficients (not MFCC)

Linear Prediction Cepstral Coefficients (LPCC) can easily be computed from LPC (Linear Prediction Coefficients) and I think that a LPC function is implemented in the same package as the MFCC.

All you need to know for extracting $p + 1$ first LPC is here. The formula that link LPC and LPCC is:

$c_0 = ln(p)$

$c_1 = a_1$

$c_i = -a_i + \sum_{n=1}^{i-1}\frac{n}{i}a_{i-n}c_n$ for $1 < i \leq p$

where $c_i$ are the LPCC and $a_i$ the LPC and $p$ the order (i.e. the number of coefficients - 1). I'm not sure about the negative sign before $a_1$. Reference here.

I'm not sure what you mean with 'windowing' but yes it may be necessary depending on your signal and application.

Hope that it will help.


FFT of sine wave not coming as expected i.e single point

Acually Matt's answer already gives one view on the problem here: the DFT is implicitly periodic in both time and frequency domain (see this question). From your parameters we can calculate that your observation period is 1 s. That means you observe 50 periods of a 50 Hz tone. Periodically extending that observation interval will always result in a seemless sine wave. If you take the 50.1 Hz tone, you're transforming 50.1 periods of an oscillation. Periodically extending that signal will result in phase jumps that cause additional spectral tributaries.

Another view is the frequency resolution that is defined as $f_\mathrm{s}/N_ \mathrm{DFT} = 1024 \text{Hz}/1024 = 1 \text{Hz}$. As 50 Hz is an integer multiple of 1 Hz, this tone can be reresented by a single DFT bin. As the result of the DFT is discrete values that are space 1 Hz apart, a 50.1 Hz cannot be represented by a single bin and the energy is smeared over the frequency domain. In your example this is the reason for the slightly asymmetric magenta spectrum and for the lower magnitude of the 50 Hz bin compared to the 0.7 magnitude of the 50 Hz spectrum.

Both effects described above contribute to the spectrum you're observing.


Anti-Aliasing and the Fourier Transform, Gonzalez Digital Image Processing

Note that in order to obtain the Fourier transform of the windowed time domain signal $f(t)h(t)$, you need to convolve the Fourier transforms of $f(t)$ and $h(t)$. We know that the Fourier transform of $f(t)$ is zero outside some interval (because $f(t)$ is band-limited). However, since $h(t)$ is time-limited, i.e., it is zero outside some time interval, we know that its Fourier transform cannot be band-limited. For this specific example of a rectangular window $h(t)$, we know that the corresponding Fourier transform is a sinc function. Convolving any function with a sinc function results in a function that extends from $-\infty$ to $+\infty$. Consequently, the Fourier transform of the windowed signal $f(t)h(t)$ also extends from $-\infty$ to $+\infty$. Consequently, the signal $f(t)h(t)$ cannot be band-limited.

In sum, a band-limited signal cannot be time-limited, and a time-limited signal cannot be band-limited. Note, however, that you cannot conclude that a signal is band-limited if it is not time-limited and vice versa. There are signal that are neither time-limited nor band-limited.


My guess is that it is an alpha filter, as defined in the context of alpha/beta filtering.


To get a nice visualisation lets consider a colour image with zero mean noise & standard deviation 20! Which means the values of each RGB component can randomly vary from +20 to -20. So pixels will be added with this noise.

Here is the code for generating noise to an Image using opencv c++

Mat mSource_Bgr,mNoise_Bgr;
mSource_Bgr= imread(FileName_S,1);

// We need to work with signed images (as noise can be
// negative as well as positive). We use 16 bit signed
// images as otherwise we would lose precision.
Mat mGaussNoise_Bgr(mSource_Bgr.size(), CV_16SC3);

double StdDev=10;
randn(mGaussNoise_Bgr, Scalar::all(0), Scalar::all(StdDev));
Mat mTemp;
mSource_Bgr.convertTo(mTemp,CV_16SC3);
addWeighted(mTemp, 1.0, mGaussNoise_Bgr, 1.0, 0.0, mTemp);
mTemp.convertTo(mNoise_Bgr,mSource_Bgr.type());

mGaussNoise_Bgr.convertTo(mGaussNoise_Bgr,CV_8UC3);
imshow("Noise Image",mGaussNoise_Bgr);
imshow("Original Image",mSource_Bgr);
imshow("Original Image + Noise",mNoise_Bgr);

Here is the Original Image

Noise Image with Mean : 0 & Std Deviation : 20

Original Image + Noise

Same Image with Noise Mean:0 Std Deviation:10


Your example seems to imply you are talking about mutually exclusive, non-overlapping regions.

In this case, given N such regions, use an image where each pixel has enough values (depth), to represent N regions. Just set the pixel value to the index of the region in that position.

You could think of the result image as a sum of N binary images, each multiplied by the $i$ subscript of the corresponding $\Omega_i$.

Alternatively, if you are asking about how to segment such an image into 2 (or more) regions, the you should use a connected component algorithm, which is common is most image processing libraries.


Hexagonal Modulator Comparison with Square Modulator

For a rectangular QAM constellation, calculation of the symbol error rate is straightforward.
Exact calculation of the bit error rate is somewhat more complicated, but approximations
are more easily calculated based on the following assumptions.

Each signal point has four nearest neighbors (interior point)
or three nearest neighbors (edge points), or two nearest neighbors (corner points),
Note that 16-QAM has 4 interior points, 8 edge points, and 4 corner points.
The probability that a transmitted signal point is demodulated into
a nearest neighbor is orders of magnitude than the probability that it is
demodulated into a more distant point
A syymbol error results in exactly one bit error (that is, data bits
are assigned to signal points in Gray code order) so that the data bits assigned
to a signal point differ from the data bits assigned to a nearest neighbor
in just one bit position.

For the 16-ary hexagonal constellation with regular equiangular hexagons
(similar to but not quite the same as the constellation shown in the question),
the 4 interior points have 6 nearest
neighbors instead of 4 while the edge points have 4 nearest neighbors on average (some
have 3 and some have 5, others have 4 exactly), and two corner points have
3 nearest neighbors.  This has two effects. First, for a given spacing
between signal points, a symbol error is more likely since there are more
neighbors. Second, a symbol error can result in more than one bit error.
Thus, it is necessary to adjust the spacing between nearest neighbors
to compensate for this if one wants to have the same BER, and this will
typically increase the signal power needed.
In short, it is not immediately obvious that

  It is supposed to give better BER measurement than Square as it has
  wider decision boundary

is a believable claim.  In any case, the width of the decision
boundary is not quite as relevant to the BER calculation as the
distance of the constellation points from the decision boundaries,
and the number of nearest neighbors that each signal point has.

I don't decipher MATLAB very well, but if the picture with the
16 blobs was created by the MATLAB program included, then I
don't think that the constellation that has been created is
a (regular) hexagonal constellation at all. It looks for all
the world like you began with a rectangular 16-QAM signal constellation
(which has points at $(\pm m, \pm n)$ for $m, n \in \{1, 3\}$)
and moved the second and fourth rows by $1$ to the right, so
that the $x$ coordinate is $-2, 0, 2, 4$ instead of $-3, -1, 1, 3$.
Note that this increases the maximum signal amplitude in
one branch (I in this case) from $3$ to $4$ so that
the amplifiers have to have greater dynamic range and
the energy per bit also has increased. Any comparison
of the BER performance must take into account the fact
that the hexagonal constellation shown is operating at
a larger SNR than the standard rectangular 16-QAM constellation.
In a regular hexagonal constellation (which is what my remarks
above apply to), the points are at the
vertices (and centers) of regular (equiangular) hexagons and
are equidistant from their nearest neighbors, while
your distances are $2$ and $\sqrt{5}$ respectively.


Why do we assume zero mean noise in sensor data?

The brutally honest answer here is: The noise is considered zero-mean because that's what the author decided to do. Without looking deeper into the signal model employed, it's impossible to answer.

However, for many systems this makes a lot of sense physically, since the processes leading to a noise realization are very often zero-mean in nature. For example, the thermal noise over any resistor should be zero-mean, simply because, well, that's how these random fluctuations work.


GCC-PHAT (Generalized cross correlation) always peak at delay=0 on real audio signal

I think that's because of the DC offset of the signal. A time domain or frequency domain plot provided by you could have been helpful here for the confirmation of your doubt. Following are suggestions

1.try to filter the signals before FFT depending upon the signal band you are interested in

2.Subtract the signal by the mean of the entire signal (which is a reasonable trick to find DC offset) or by the DC offset directly if you know that.

The signal of your's must have less swing amplitude in comparison to DC offset making the '0' point correlation of greater amplitude than the actual delay.


Multiplying the imaginary part of DFT with a linear ramp to get a derivative

I've changed my mind during the answer. Instead of using frequency domain processing to get the derivative signal, I will provide a more direct and clean approach to get the derivative of a signal by using what's called a digital-differentiator.  An LTI filter, which is applied into the samples $x[n]$ of a signal $x(t)$ and computes the samples $y[n]$ of the derivative $y(t) = x(t)'$ of $x(t)$.
Below is the matlab code you can use to compute the samples of the derivative of a bandlimited signal $x(t)$:
N = 512;                   % length of x[n]; samples of signal x(t).

% Assume Gaussian pulses added:
t = linspace(-1,1-2/N,N);  % time is sampled at T = 2/N period.

sigma = 0.001;             % Gaussian pulse shape parameter (adjust accordingly)
x1 = 2*exp(-(t.^2)/(2*sigma) );         % Gaussian pulse one
x2 = 0.5*exp(-((t-0.5).^2)/(3*sigma) ); % Gaussian pulse two
x = x1 + x2;               % total signal

% Generate differentiator impulse response hc[n]:
T = 2/N;                  % sampling period (implicit, defined above)
n = -N/2:N/2;
h = ((-1).^n)./(n*T);
h(N/2+1) = 0;
w = hamming(2*N/2+1)';
hc = w.*h;                % causal,stable impulse response hc[n]

% Filter x[n] to compute samples of derivative signal
y = conv(x,hc);
yd = y(N/2+1:N/2+N);     % extract relevant portion

% display the result:
figure,plot(t,x);
hold on
plot(t,0.1*yd);
legend('x(t)','dx/dt');
grid on
title('x(t) vs dx(t)/dt  (scale adjusted to fit both');

The result is :

For getting the same result using DFT techniques,
the following code can be used :
M = N+N;
H = fft(hc,M);          % generate the RAMP mask to apply on X(w)
Y = fft(x,M).*H;        % obtain ramp multiplication in freq domain
y = real(ifft(Y,M));    % invert back to time domain
y = y(N/2+1:N/2+N);     % extract relevant portions

figure,subplot(2,1,1)
plot(linspace(0,2,M),abs(H)); title('RAMP mask that is applied on X(w) in freq-domain')

subplot(2,1,2)
plot(t,x);
hold on
plot(t,0.1*y);
legend('x(t)','dx/dt');
grid on
title('x(t) vs dx(t)/dt  (amplitude scale adjusted)');

with exactly the same result on the time-domain sequence as shown below:


Relaxation time in a stationary stochastic process

Considering the complexity of the question, and the rating of your reputation, i assume you are still alive.

In the following, i am assuming $k$ is a constant, since it is not a variable -fixed or random- expressed in the question.

If you are still interested, which is a numerically strong assumption, the system can be rewritten trivially as:

$$\frac{d}{dt}
\left(\begin{array}{cc}
x_1  \\
x_2
\end{array}\right)=
\left(\begin{array}{cc}
-(t_1+d) & t_0 \\
t_1 & -(t_0+d)
\end{array}\right)
\left(\begin{array}{cc}
x_1  \\
x_2
\end{array}\right)
+
\left(\begin{array}{cc}
k  \\
0
\end{array}\right)
=
\left(-d
\left(\begin{array}{cc}
1 & 0  \\
0 & 1
\end{array}\right)
+
\left(\begin{array}{cc}
-1 & 1  \\
1 & -1
\end{array}\right)
\left(\begin{array}{cc}
t_1  \\
t_0
\end{array}\right)\right)
\left(\begin{array}{cc}
x_1  \\
x_2
\end{array}\right)
+
\left(\begin{array}{cc}
k  \\
0
\end{array}\right)$$

Which has the standard form:
$$Dx=Ax+x_0$$

The constant vector $x_0=[k;0]$ or equivalently the constant vector $A^{-1}x_0$ is a state reference.

Through the variable change:
$$x+x_0 \to x$$
you can safely discard the $x_0$ vector.

This is an autonomous system which, unless you forget to include perturbations or inputs, is a very boring system. I am sure you did.


The impulse response and the transfer function of an LTI system are different, sometimes extremely different (except for a specific shape of Gaussian).

The DC component of any impulse response is it's average value.  Your example time-domain squiggle seems to have a non-zero average value (integrate to determine the area under the curve), thus is has a non-zero DC component in the frequency domain.


Here is the code for adding the two signals(the delayed and not-delayed):

    x = getaudiodata(recObj); n1 = 1:size(x,1);%audiodata of original signal
    y = time_delay(x , 50000 ); n2 = 1:size(y,1);%audiodata of delayed signal
    mixed = sigadd(x,y,n1,n2); %audiodata of mixed signal
    mixrecObj = audioplayer(mixed,44100);
    play(mixrecObj);

This is the sigadd function:

function [ y,n ] = sigadd( x1,x2,n1,n2 )
%implements y[n] = x1[n] + x2[n]
%y = sum sequence over n wich includes n1 and n2
%x1 = first sequence over n1
%x2 = second sequence over n2
%

n = min(min(n1),min(n2)):max(max(n1),max(n2)); %duration of y(n)
y1 = zeros(1,length(n)); y2 = y1;
y1(find((n >= min(n1))&(n <= max(n1)) == 1)) = x1;
y2(find((n >= min(n2))&(n <= max(n2)) == 1)) = x2;
y = y1 + y2;
end

This works perfectly, but is not enough for getting a real echo effect. In order to archieve this we have to
create more delayed signals(with diferent delays of course), and then add them all,like this:

x = getaudiodata(recObj); n1 = 1:size(x,1);
y = tim_delay(x , 5000 ); n2 = 1:size(y,1);
s = time_delay(x , 4000 ); n3 = 1:size(s,1);
d = time_delay(x , 3000 ); n4 = 1:size(d,1);
mixed1 = sigadd(s,d,n3,n4);n5 = 1:size(mixed1,2);
mixed2 = sigadd(x,y,n1,n2);n6 = 1:size(mixed2,2);
totmixed = sigadd(mixed1,mixed2,n5,n6);
mixrecObj = audioplayer(totmixed,44100);
play(mixrecObj);

I have tested this one with the three delays plus the original and you can hear a pretty nice echo.

Note: In the third piece of code it would have be better if there were a function that returned the signal with the echo by passing an echo effect index  as a parameter to it, but I left it that way for the sake of simplicity.


Using only positive frequencies in fourier domain, How will it affect the ifft?

You can do things that way, but unless you're manipulating things properly, you won't get the right answer.  A procedeure for doing a N=32 point example is:

1) generate 32 point time domain data,
2) forward transform it,
3) divide results by 32,
4) multiply positive frequencies by 2 (frequencies f= 1 to 15),
5) zero negative frequencies (frequencies f=17 to 31),
6) inverse transform it,
7) ignore the imaginary outputs.

The real results of the above procedure will match the original input waveform.

Note that in the above, you don't multiply the f=0 or f=N/2 outputs of the FFT by 2.  For this case (N=32 - an even number), they are unique - they don't have negative frequency counterparts.  Now if you're using an appropriate FFT routine that allows for odd N, then you have to account for the fact that you'll have an f=0 point, but no f=N/2 point.

Unless you are doing very large FFT's, you're probably not going to save much in the way of processing time by using just the positive frequencies.


Find horizontally and vertically aligned sets of points

Using alignment, it's mostly a matter of heuristic (sweeping horizontally and vertically). You can also work out something similar to Hough transform here.

Using inter and intra-dot distances is also known as linkage. What you want to do can be obtained through hierarchical clustering using a single linkage distance computation. It works as follows:

you start with all your dots as (unconnected) leaves in a tree;
you test all the leaves pairs, compute their spatial distance, and connect the two leaves that yield the smallest distance;
replace these 2 leaves by only 1 leave that represents the cluster;
go back to 2 and repeat until you have only 1 leave left.

The last leave is the root of a tree that represents the hierarchical clustering result, which is also called a dendrogram.
By using single linkage, you connect each dot to the closest dot, then recursively all closest pairs of dots, etc. You can finally discard the hierarchical tree part that corresponds to the smallest distance between two different lines.