The two main measures of scatter in the Gaussian case are the mean absolute deviation:

$$d_n = 1/n\sum |x_i-\overline{x}|$$

and the mean square deviation:

$$s_n = \sqrt{1/n\sum (x_i-\overline{x})^2}$$

If the random Gaussian noise has deviation $\sigma$, $s_n \to \sigma$ and   $d_n \to \sqrt{\frac{2}{\pi}}\sigma$. I don't have the exact references right now, you can find them in Robust Statistical Procedures, 1996, P. Huber, or texts by Hampel, or Rousseeuw, on robusts statistics.

The derivation is OK, let us use the zero-mean case: the expected value of the absolute value of the random variable:

$$E ( | x | ) = \int_{- \infty}^{ \infty} | x | p ( x ) d x$$

then you derive:

$$E(|x|) =
2\int_0^{\infty} \frac x {\sigma\sqrt{2\pi}}
e^{-\,\frac {x^2} {2\sigma^2} } dx$$

then

$$E(|x|) =
\left[-\,\sigma \sqrt{\frac 2 {\pi}}  e^{-\,\frac {x^2} {2\sigma^2} }\right]_0^{\infty}
=  \sigma \sqrt{\frac 2 {\pi}}$$

But this is for the mean absolute deviation, I am not sure of the equivalence with  the Median Absolute Deviation.


How can I estimate the spatial resolution of consumer grade electronics?

Consider the following figure:

The cameras are separated by baseline 2b - say 10cm, so b=5cm.  Your focal length of 4" sounds looong (about 100mm) for a webcam, but hey, let's go with it.  f=100mm (or also 10cm!).

From $xr-xl = \frac {2bf}{z}$ an object at 10m away will have a pixel difference between the two cameras ($xr-xl$) of 0.1*0.1/10 = 0.001m.  So, if your pixels are 5 micron square $xr-xl$ is about 200 pixels.

Which is a lot - large disparities are expensive to search for.

With a more normal focal length of say 5mm, your disparity is more like 10 pixels.  11m is about 9 pixels, so your distance resolution is about 1m @ 10m.

If you have a specific camera in mind you can run the numbers for your particular geometry.


Subtracting sensors output that have mutual inductnace between them

This sounds like blind source separation.  In general, you can't unmix things after they've been mixed. If you have two different recordings of two sources with some of each source in each recording, you can sometimes use independent component analysis to separate them.

I have a Python example here.  There's also FastICA for MATLAB.  If they are audio signals, being picked up magnetically, there's probably no appreciable delay between them.  ICA works well in this case.


What's the frequency transform of a random matrix?

That is the DC offset (at zero frequency).

img = 255*randn(512);  % thanks @Jason's comment
img1=img-mean(img(:));
fftimg = abs(fft2(img1));
imagesc((fftshift(fftimg)))

gives you the spectrum of the normal distribution image. Note that white noise is flat in the frequency domain, but the Gaussian noise in the time domain is still Gaussian in the frequency domain. You may need to average a large number of FFTs of white noise to approach the average power spectral density.


...take the FFT of the image, and discard the DC (constant) components.

That would result in removing the mean brightness from an image. Due to the fact that grayscale images have positive numbers between 0-255, what would be returned from that "high pass" filter (that only "cuts" the DC) would be the image with the same grayscale variance but now centered around 0 (instead of its original mean level).

  I want to analyze how much these fringes are moving around. Is there an easy way to do this?

You could obtain the absolute difference between two successive frames and then the sum of all values of that quantitty as a simple metric of how much movement was there.

If you would also like to be able to estimate how much displacement was there between two successive frames then the simplest thing would be to use cross correlation (on successive frames) and track the position of its maximum.


GnuRadio CC Decoder or How to use FEC Extended Decoder

Inside CCSDS channel coding standard, there is a synchronization word called the ASM.
It is a pattern sequence of 32, 64 or 256 bits which is used for synchronization inside the stream. It is usefull to solve the two points I describe below.
Inside your GNU block diagram I see two things that might be an issue:

Phase ambiguity: There is no mechanism to solve the phase ambiguity of the BPSK in demodulation. I am not really sure if Viterbi decoder works fine if not solved (the ambiguity).

Two bits ambiguity of CC: This type of code with a rate of 1/2 produces two encoded bits for one noncoded. At the input of decoder, you must be synchronized with the first encoded bit of a couple.

This could explain that sometimes your system works after switching on/off. On chance over two or maybe four...
Other remark, using an alternating pattern of 0 and 1 could make you miss something debuging. Try a 32 bits pattern.
Tell us if you found something


What does following statement mean: " Simulation was done at the rate of 12 samples per cycle of 50Hz a.c."

One cycle duration of 50Hz is 0.02seconds. So this means you need to take 12 samples in a duration of 0.02s. In one second, there will be 600 samples (sample rate of 600 samples per second).


Why is bit error rate = symbol error rate / number of bits per symbol in QPSK

In general, the bit error probability depends on the symbol error probability and the coder, i.e., the mapping between bits and symbols. The formula you refer to

$$P[\textrm{bit error}]\approx\frac{1}{M}P[\textrm{symbol error}]\tag{1}$$

with $M$ being the number of bits per symbol ($M=2$ for QPSK) is an approximation for a certain type of coder under the assumption that the SNR is relatively high.

In $(1)$ it is assumed that the system implements a Gray code, where nearest neighbor symbols only differ by one bit. So if the SNR is high, most symbol errors occur due to confusing a symbol with one of its nearest neighbors, which causes only one out of $M$ bits to be wrong.


Morphological transformations for images with discrete values

Morphological operators are primarily defined over binary images where values 0,1 could correspond to areas identified as foreground and background respectively.

In your case, you can create a binary image of the region you want to apply the operator on by assigning the "islands" of noisy classified regions to the "foreground" (e.g 4  now becomes the value 1) and the broad area around them that is consistently classified in a stable way to the "foreground" (e.g value 3). Apply erosion to this image and then substitute the pixels back from the new binary image to the segmentation mask.

Alternatively, if you want to eliminate spurious classifications so that you end up with consistent areas, you could apply nonlinear filters such as median, max and min over some sliding m by n window.

For more information, please see this link and this link

As far as Python is concerned please see this or this and this.  Another excellent python module for this job is this one.

(Also, please note that you could apply any kind of filter that could even produce a Real type output (e.g 2.43, 1.76, etc) but either round or truncate its final output to make it to an integer. Not sure these operators would be very helpful though in this specific case.)

Hope this helps.


Some examples would be the null to null bandwidth(s) as well as the sidelobe levels of the first/second sidelobes.

Example vector of 5 numbers:
[mainlobe gain, first null location, 1st sidelobe gain, second null, 2nd sidelobe gain]

From this you could produce an estimate for the 3 dB points of each lobe using some "reasonable" formula, i.e. halfway between the peak and the null.


Orthonormal Dictionaries for Band Limited Signals

Why would you add the constraint of being Orthonormal Dictionary?
It doesn't make sense in the context of what you ask.

First we need to define resolution.
If you mean the grid to be denser than indeed what you need is to create the DFT matrix of zero padded signal and take the subset which you're interested in as you wrote in your question.

If you define resolution as the ability to discern between 2 close (In frequency relative to the observation window) harmonic signals then there is a subtle thing to keep in mind.
While indeed by Uncertainty Principle you're limited by the observation time what's important is under what model.
If you have no knowledge about the signal but it is band limited than this holds and probably there is nothing to do.

Yet if you have some prior on it you can do something and indeed gain some ability to "Bypass" the limitation.

For instance, if you assume the signal is sparse relative to the DFT matrix than solving the problem with $ {L}_{1} $ regularization (Sparse Model / MAP Model with Laplacian Prior) will yield ability to have "Super Resolution".

I can't, currently, think on a prior to work on range of resolutions, but if you believe your signal is sparse in that range you certainly can build an optimization problem to gain extra resolution.


What is the relationship between the PSD of a continuous signal and the PSD of its periodically sampled one?

If you sample a finite-power continuous-time random process $x(t)$ you get a discrete-time random process $y_k$. If $x(t)$ is wide-sense stationary (WSS) you get for the autocorrelation function of $y_k$

$$R_y(k,l)=E\{X(kT)X^*(lT)\}=R_x((k-l)T)=R_x(mT)$$

where $T$ is the sampling period. Obviously, $y_k$ is also WSS (it only depends on the difference $m=k-l$), and its autocorrelation function is a sampled version of the autocorrelation $R_x(\tau)$ of the original continuous-time process. Since the power spectrum is the Fourier transform of the autocorrelation function, the relation between the two power spectra is the same as the relation between the spectrum of a continuous signal and its sampled version:

$$S_y(e^{j\omega T})=\frac{1}{T}\sum_{k=-\infty}^{\infty}S_x(\omega-2\pi k/T)\tag{1}$$

Since the anti-aliasing filter has a cut-off frequency of $\pi/T$ it is sufficient to consider the interval $-\pi/T\le\omega<\pi/T$ (i.e. just consider the term $k=0$ in (1)):

$$S_y(e^{j\omega T})=\frac{1}{T}S_x(\omega)\tag{2}$$

This equation is equivalent to Eq. (10.65) in your question. Just note that I've used $\omega$ in a different way (analog frequency in radians), whereas Oppenheim uses $\omega$ as the normalized angular frequency. So with $\omega'=\omega T$ you finally get

$$S_y(e^{j\omega'})=\frac{1}{T}S_x\left(\frac{\omega'}{T}\right)$$


Characteristic function of a random Gaussian variable

In general, the characteristic function of a random variable is related to the fourier transform of the distribution as follows:

$$\varphi_Z(-\omega) = \mathscr F \big\{ f_Z(z) \big\}$$

Why?

Because, fourier transform of a PDF $f_Z(z)$ is:

$$\mathscr F \big\{ f_Z(z) \big\} = \int^{\infty}_{-\infty} f_Z(z)\cdot e^{-j\omega z} \ \mathrm dz$$

And Characteristic function of a PDF is given by:

$$\varphi_Z(\omega) = \mathbf E_Z \big\{e^{j\omega z} \big\} = \int^{\infty}_{-\infty} f_Z(z) \cdot e^{j\omega z} \ \mathrm dz$$

As you can see the difference is just in the argument $\omega$. I think now you can get the characteristic function of a Gaussian Random Variable. Just calculate the fourier transform of the PDF and replace $\omega$ by $-\omega$.

Fourier Transform of Gaussian PDF :

$$\mathscr F \big\{ f_Z(z) \big\} = \int^{\infty}_{-\infty} \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{1}{2\sigma^2}(z-\mu)^2}\cdot e^{-i\omega z} \ \mathrm dz$$

Now, just consider the power of exponential and simplify as follows:
$$(\frac{-(z-\mu)^2}{2\sigma^2} - i\omega z )= -\frac{1}{2\sigma^2}(z^2+\mu^2-2\mu z + i2\sigma^2 \omega z)$$ $$= -\frac{1}{2\sigma^2}(z^2 -2(\mu - i\sigma^2 \omega)z + \mu^2 + (\mu - i\sigma^2 \omega)^2 - (\mu - i\sigma^2 \omega)^2)$$ $$ = -\frac{1}{2\sigma^2}((z - (\mu - i\sigma^2 \omega))^2 - (\mu - i\sigma^2 \omega)^2 + \mu^2)$$

Now, basically, the integral becomes the following:
$$\mathscr F \big\{ f_Z(z) \big\} = \int^{\infty}_{-\infty} \frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{1}{2\sigma^2}(z-(\mu - i\sigma^2 \omega))^2}\cdot e^{\frac{1}{2\sigma^2}((\mu - i\sigma^2 \omega)^2 - \mu^2)} \ \mathrm dz$$ The second exponential can come out of the integral since it does not contain $z$, then what you are left with is just area under a Gaussian PDF with mean $(\mu - i\sigma^2 \omega)$ and variance $\sigma^2$. Area under any PDF is equal to $1$. So, basically the fourier transform of the Gaussian PDF becomes:
$$\mathscr F \big\{ f_Z(z) \big\} = e^{\frac{1}{2\sigma^2}((\mu - i\sigma^2 \omega)^2 - \mu^2)} = e^{-\frac{\sigma^2 \omega^2}{2} - i\mu \omega}$$
Once you have fourier transform, you can get the characteristic function by substituting $\omega = -\omega$ as follows: $$ \varphi_Z(\omega) = e^{-\frac{\sigma^2 \omega^2}{2} + i\mu \omega}$$


Yes, the 30 kHz are correct, but only assuming your receiver has filters that filter out any noise outside these bands. That might or might not be the case – if your receiver in fact observes 0 – 45 kHz, for example, then the noise energy it'll see is defined by that 45 kHz bandwidth.

I'd also recommend you not just peruse the formula from Wikipedia (which I'd argue is in terrible shape), but try to understand how this noise floor calculation comes to be, and then, after that, apply the conversion to dBm.

Namely, you get Johnson-Nyquist thermal noise. This noise has a power density of $k_B T$, that is, Boltzmann's constant times the temperature of your system. This constant represents Energy per temperature above absolute zero (eg. Joules per Kelvin), so $k_B T$ is Energy.

In signal-processing systems we commonly care about Power more than energy (because, for example, if we know how much power gets dissipated by a resistor, we know the voltage across it), so we divide the Energy by time, or, multiply with bandwidth/frequency (which is the inverse of time):

$$P_\text{N, thermal}=k_B T \Delta f$$

Real-world active systems, especially receiver's amplifiers, add more than just the noise present due to thermal noise – the factor between the noise going in and coming out is the Noise Figure $\text{NF}$. Therefore:

$$P_\text{N, thermal}=k_B T \Delta f\, \text{NF}$$

Because we like our $\text{dB}$-based units, let's cancel out all units (I'm doing this in-line here just for clarity) and then do the dB-conversion. Then we get

$$P_\text{N, thermal}[\text{dBW}]=10 \log_{10}\left((k_B \cdot \frac KJ)( T \cdot \frac 1K)(\Delta f \cdot s) \, \text{NF}\right)$$

You wanted $\text{dBm}$, but a milliwatt is simply $\frac1{1000}$ of a Watt, so

$$\begin{align}
P_\text{N, thermal}[\text{dBm}]=10 \log_{10}\left((k_B \cdot \frac KJ)( T \cdot \frac 1K)(\Delta f \cdot s) \, \text{NF}\,\cdot 1000\right)
\end{align}$$

Let's apply a few logarithm rules here:

$$\begin{align}
P_\text{N, thermal}[\text{dBm}]&=10 \log_{10}\left((k_B \cdot \frac KJ)( T \cdot \frac 1K)\right)&+&10 \log_{10}\left(\Delta f \cdot s \right)&+&10 \log_{10}\left( \text{NF}\right) &+& 10\log_{10}\left(1000\right)\\
&=10 \log_{10}\left((k_B \cdot \frac KJ)( T \cdot \frac 1K)\right)&+&10 \log_{10}\left(\Delta f \cdot s \right)&+&\text{NF}_\text{dB} &+& 30
\end{align}$$

If you drop in the values of Boltzmann's constant and room temperature, you'll end up with the classical -174dBm/Hz; simply add your bandwidth (in dBHz) and your Noise Figure (in dB), and you get your observed noise power.


Is RMS of a sound wave good measure of loudness?

No

The sensitivity of human auditory is very dependent on frequency and without any type of frequency weighting, a raw RMS estimate will be pretty much useless.

There are some standard frequency weighting curves, most popular being the "A" weighting.  See https://en.wikipedia.org/wiki/A-weighting. The sensitivity is also dependent on the absolute level and hence there are multiple weighting curves.

Sp frequency weighted RMS is a good starting points, but they typical ignore any sensitivity to time domain behavior like impulse, transient, masking, etc. There are good models to predict actual perceived loudness but they tend to be rather complicated.


How to store the signal in a variable and plot the signal in MATLAB

I suggest you to do something like (one plot with two curves):

figure
plot(sound_L(:,1), 'b'); % blue
hold on; grid on; % prevent from redrawing and add grid\
plot(sound_R(:,2), 'r'); % red
legend({'left', 'right'});

Or if you want two subplots:

figure
subplot(211);
plot(sound_L(:,1), 'b'); grid on;
subplot(212);
plot(sound_R(:,2), 'r'); grid on;

Please mind that in your case sound_L and sound_R are the same. You are not making any use of beep variable.


Why is state-space representation more numerically stable than zeros-poles representation?

This is because polynomial root finding algorithms are terrible without any implication that people code it wrongly or else. It is an ill-posed problem whereas eigenvalue computations and hence generalized eigenvalue problem solutions are much more stable.

Any mistake in any of the pole/zero data permeates to all entries of the polynomial entries and contaminates them which in turn causes big numerical inaccuracies in the polynomial root finding.

Plus, MIMO system representations in both are simply asking for trouble since common terms are very difficult to retain and after a few manipulations you can easily create dummy roots -which could have been considered to be numerical noise- becoming actual modes of the representations.

Similarly, canonical representations suffer from the same problem in state space. You can test it via creating 2 identical A matrices and perturb the entries of the row/col that holds the entries. Then create a balanced realization of the same system and perturb the A matrix entry of that representation and you'll see that the difference is usually very big.

In your comment you mention the nonuniqueness of the SS representation and it is true. But the dynamics are identical, hence they form an equivalence class up to transformations with nonzero transformation matrix $T$. That's actually a feature not a bug.


Fourier Transform: interpretation of continuous spectrum at specific frequencies

You can see the non-zero magnitude of the spectrum, but the interpretation that the signal contains a sinusoidal component at a specific frequency where the magnitude is non-zero is wrong. If there is a sinusoidal component present, then there is a Dirac impulse at the respective frequency. A non-zero Fourier transform at a certain frequency is not sufficient.

The squared magnitude of the Fourier transform is the energy spectral density, and the energy inside a certain frequency band can be obtained by integrating the energy density over the respective interval.

An analogy from statistics is the difference between probability and probability density function (PDF). The values of the probability density function are not probabilities, but probabilities can be obtained from the PDF by integrating it over certain intervals.


Continuous phase for phase delay calculus in FIR filters

The problem is that magnitudes close to zero in the stop band change their sign almost arbitrarily / randomly.  And that's why you're seeing arbitrary jumps of $\pi$ radians --- it's a sign change.

Because of that, I'd be tempted to just ignore the phase of any part of the magnitude response below some low value. What value to use will depend on your application.

Zero padding will not generally help because that just gives you more "random" sign changes around the frequency axis.

I use "random" in quotes as it's not stochastic, but it is hard to predict.

If I interpolate the phase on 100,000 points for an FIR filter of length 126 taps and unwrap the phase I get the picture below.

The reason the unwrapping algorithm fails is because the phase in the stop band is changing very quickly compared with that in the passband. Comparison below of unwrapped phases.  The red lines are $\pm \pi/2$.

Changing the complex_magnitude expression to mimic the fftshift idea from @hotpaw, to:

    complex_magnitude <- polyval(f1[c(63:125,1:62)], exp(-1i*freq))

I get:

which is slightly better, but not ideal.

R Code Below

#29899

library("signal")

f1 <- remez(125, c(0, 0.3, 0.4, 1), c(1, 1, 0, 0))
plot(f1,type="l")
dev.copy(png, 'Q29899-ImpulseResponse.png')
dev.off()

freq <- seq(0,pi,pi/100000)

complex_magnitude <- polyval(f1, exp(-1i*freq))

plot(freq, unwrap(Arg(complex_magnitude)), type="l");
dev.copy(png, 'Q29899-UnwrappedPhase.png')
dev.off()

plot(freq[60000:100000],Arg(complex_magnitude[60000:100000]), type="l")
lines(freq[c(60000,100000)],c(pi/2,pi/2), col="red")
lines(freq[c(60000,100000)],c(-pi/2,-pi/2), col="red")
title('Stop band phase detail')
dev.copy(png, 'Q29899-PhaseDetailWrapped.png')
dev.off()

plot(freq[1:10000],Arg(complex_magnitude[1:10000]), type="l")
lines(freq[c(1,10000)],c(pi/2,pi/2), col="red")
lines(freq[c(1,10000)],c(-pi/2,-pi/2), col="red")
title('Pass band phase detail')
dev.copy(png, 'Q29899-PhaseDetailWrappedPassband.png')
dev.off()


Checking if a lecturer can be heard clearly from the back of a room

There's possibly an easy way out: Add two mics, one near the speaker and one in the back of the room. Run both signals through speech recognition software. If there is textual output from both and they agree, then the voice must be loud and clear. The software must be able to handle the language spoken though.

Other than that, there is research and standards on speech intelligibility. Perhaps it is better to start with a simplified model of intelligibility. Two important factors are speech level and noise level. By carrying out preliminary listening tests and measurements in the space, you can determine what the baseline noise level ($N_0$ dB) in the back is and how loud the speech must be ($S_0$ dB) in the minimum at the front mic for it to be intelligible in the back. Then during a lecture, you would detect voiced segments of the audio with the front mic and measure the speech level ($S$ dB) at the front during voiced segments and noise level ($N$ dB) in the back during non-voiced segments. If $S-N \ge S_0-N_0$, the speech should be intelligible. When measuring levels, you should weight the spectrum (by running the signal through a filter) so that level best reflects speech intelligibility or the detrimental quality of the noise. You would need to look into voiced/non-voiced detection, or simply have a threshold level for $S$ above which you classify the input as voiced.


How to implement Fourier Descriptor of an image?

Yes, these k are the same.

In the basic sense, the rotational invariance is obtained by only taking the magnitude of the descriptors. However, in practice this might lead to information loss and decreases discrimination power. For this reason it is typically prefered to augment the phase information with additional terms. There are many ways to do this but for a simple approach follow:

Invariant Fourier Descriptors Representation of Medieval Byzantine Neume Notation

Dimo Dimov, Lasko Laskov


If the cyclic autocorrelation output has the properties you need in your application (e.g., it's peak-to-average ratio is high enough for picking out the peaks), then it's probably okay to substitute the cyclical correlation for the linear one. For example, consider a length 7 maximal-length sequence:
$$
\{+1, -1, -1, +1, +1, +1, -1\}
$$
The cyclical and acyclical autocorrelation functions are plotted below.

Notice that both functions have a peak value of 7. The acyclical autocorrelation function has a largest off-peak magnitude of 3, whereas the cyclical autocorrelation function has a largest off-peak magnitude of 1. The peak-to-off-peak amplitude ratios are $7$ and $7/3$ for cyclical vs. acyclical, respectively. This means that the cyclical autocorrelation function is basically three times more peaky, making it much easier to detect. If the extra peak-to-off-peak amplitude ratio is not necessary for your application, it would be okay to correlate using linear correlation instead.

If instead, you have something whose acyclical autocorrelation properties are optimal (such as a Barker code), and you would like to calculate that acyclical correlation but don't want to use the computationally inefficient direct form of convolution, there are algorithms that calculate the acyclical autocorrelation using the computational savings of FFTs. Two of these I know of are called the overlap-add and overlap-save methods. They work on blocks of the input data at a time, and using FFTs of the blocks of data and multiplication in the transformed domain, calculate the acyclical correlations/convolution of the data against your given sequence.


Convert sine wave signal to square wave signal

What you can do is process this signal with a windowed RMS function (I'm writing in matlab code, but it can easily be translated in python) :
windowlength = 5; %you can change that, but keep it odd
hw = (windowlength-1)/2; %hw stands for "half window"
%pre-assigning the output RMS signal
sig_RMS = zeros(1,length(sig));
%I'm adding zeros the the signal's left and right for boundary calculation
sig_zero_padded = [zeros(hw,1), sig, zeros(hw,1)];

sig_zero_padded = sig_zero_padded.^2; %squaring the samples
for i=(hw+1):(length(sig_zero_padded)-hw)
   % calculating square-mean inside the window and assigning the result
   %to the appropriate sample in sig_RMS
  sig_RMS(i-hw) = sqrt(mean(sig_zero_padded(i-hw:i+hw)));
end

Then you set a minimum value for which superior values of the RMS signal will be assigned to 1, and the other to zero :
trigger = 0.5; %I chose arbitrarily
square_sig = (sig_RMS > trigger);

And you should get your square signal.
What we are doing here, in simple language, is taking little bits of the signal and calculating the RMS value of each bit. If this calculated value is significant enough to be considered as "not noise", we assign '1' to it.
EDIT :
I managed to compute it and it worked. Note that I had to take Fs = 10^5 because my computer does not have enough capacity for Fs = 10^6.
Here is the signal (blue), and the calculated RMS signal (red) for windowlength = 8191 :

Then I set trigger = 0.4 and obtained two crystal clear zones identifying the useful signal :

So in conclusion :

Set windowlength to get the RMS signal mean value far enough of the noise while keeping accuracy. In my case it was approximately 1 % of the total number of samples.
Set trigger between the RMS crests and the noise

However, I'm having that feeling that I'm reinventing something so feel free to look up for some literature about this method, I'm sure someone did the math to obtain optimal values of windowslength and trigger. ;)


Calculating $E_b/N_0$ from given SNR ratio

HINT:
$$
\frac{S}{N}=\frac{E_b \cdot R_b}{N_0 \cdot B }
$$
With all units linear:
$$\implies \frac{E_b}{N_0} = \frac SN\cdot \frac BR_b\tag{linear}
$$
Now with $S/N$ in $[\rm dB]$ and $R_b$ and $B$ linear, you get $E_b/N_0$ in $[\rm dB]$ as follows:
$$\implies \frac{E_b}{N_0} = \frac SN +10\log_{10}\left(\frac BR_b\right)\tag{in [dB]}$$


You must use notch-stop or band-stop filters.

There are different methods to implement them, But you could implement them easily using fdatool in MATLAB. After you designed the filter, filter your signal using filter function.


To have the desired effect of image denoising, the Wiener filter for can only be implemented with a point-spread function (blurring filter) that filters out the noise by averaging in some neighborhood of image pixels. As you do not disclose your program text, one cannot say for sure, but most probably your disappointing result comes from failure to apply the averaging in the vicinity of each pixel.
The typical derivation of the Wiener filter in Fourier domain arrives at the expression
$$
H(u,v) = {\frac {B^*(u,v)S_{ff}(u,v)} {|B(u,v)|^2 S_{ff}(u,v) + S_{nn}(u,v)}}
$$
where $H(u,v)$ is the Wiener filter in agreement with your designation, $B$ is the point-spread function (blurring filter), and I rename your $P$'s into $S$'s in order to reserve the character $P$ for a common PSD. In keeping with your use of $f$ subscripts for the original image and $g$ for this image with WGN added, $S_{ff}(u,v)$ is the (supposedly unknown) power spectrum of the original image before the noise was added.
Using the periodogram technique to estimate the power spectrum from the DFT of the observation $X(u,v)$, one can write $S^{periodogram} = |X(u,v)|^2/N^2$ and then arrive at the expression similar to your $H(u,v) = (P_g - P_n)/P_g$:
$$
H(u,v) = {1 \over B}{\frac {S^{periodogram}(u,v) - S_{nn}}{S^{periodogram}(u,v)}}
$$
with the remarkable appearance of the ${1 \over B}$ factor.
Now, there can arise problems with the inverse of the point-spread function $B(u,v)$. The solution is to use matrix pseudoinverse.
Back to your question: you failed to reconstruct a denoised image because you omitted the averaging in the pixel neighborhoods. Experimentation with MATLAB function wiener2(I,[m n]) for the pixel-wise adaptive low-pass Wiener filtering with different sizes of the pixel neighborhood would help you understand the working of the Wiener filter on noisy images. You would see how the choice of optimal m,n values depends on the local image variance.


How can I plot the spectrum of a signal in MATLAB?

One way is simply to plot the absolute value:

X_total = fft(x_total);
figure;
plot(10*log10(abs(X_total).^2));

Note that if your signal is real, you can plot only the first part of the vector.
Another option is to use the Welch Periodogram

figure;
pwelch(x_total);


Analogous Library to OpenCV for Audio Processing / Analysis

Consider the following:

clam-project.org:

  CLAM (C++ Library for Audio and Music) is a full-fledged software framework for research and application development in the Audio and Music Domain. It offers a conceptual model as well as tools for the analysis, synthesis and processing of audio signals.

MARF:

  MARF is an open-source research platform and a collection of voice/sound/speech/text and natural language processing (NLP) algorithms written in Java and arranged into a modular and extensible framework facilitating addition of new algorithms. MARF can run distributedly over the network and may act as a library in applications or be used as a source for learning and extension.

aubio:

  aubio is a tool designed for the extraction of annotations from audio signals. Its features include segmenting a sound file before each of its attacks, performing pitch detection, tapping the beat and producing midi streams from live audio.


Equivalence of Maximum Likelihood (ML) and Discrete Fourier Transfrom (DFT) Peak Finding for Single Tone Estimation

Maximum Likelihood under the assumption of Additive White Gaussian Noise (AWGN) is always equivalent to finding the hypothesis with the minimum distance to given data.

Since minimizing distance is equivalent (In the euclidean Space) of maximizing the correlation you can always build the idea of Match Filter for parameter estimation in the settings of ML under AWGN.

Matched filter means that you create a set of hypothesis parameterized on a grid and goes one by one until you find the one which maximizes correlation (Minimizes distance) to data.

This is exactly what's the DFT is doing.
The grid is based on the different tones (Frequencies) and the correlation level is given by the bin value.


High Level Process, Calculating Cepstrum When You Only Have Spectrum

The cepstrum is defined as:

$C = | \mathcal{F} \log( P(f) ) |^2$

$\mathcal{F}$ is a Fourier transform and $P(f)$ is the power spectrum.  If you've got the power spectrum for the signal there is no need to go back to the time domain.  Just take the logarithm, compute an FFT, take the magnitude squared of the result.

Also, the data you list is strictly positive and real so I'm assuming it's the spectrum.  At the end of your question you put spectrum (FFT).  If what you had was an FFT of a real-valued time domain signal, and not a spectrum, I would expect it to have an imaginary part.


extending windowing function to higher dimensions

Without being familiar with the library they use, this pretty much looks to me like a numerical normalizing condition.

Why you need this

This is needed in order to not change the magnitude of the numbers in your data set after filtering. Consider for example a large (multidimensional) region in your image data that exhibits the same value everywhere, e.g. 529. If you filter/smooth this dataset with a filter that has a FWHM that is considerably smaller than the size of the region, you would expect that the value in this region does not change. If you wouldn't include this normalizing condition, this behavior is not guaranteed and you could end up with totally different values... e.g. 1E6 or 0.02 instead of 529. This would seriously complicate further computation, if you had a threshold or so in later processing steps. Or if you converted your data later on... The values would always depend on the filter settings.

A somewhat analytical example

Consider a 1D Gaussian curve. It has a normalization factor $\displaystyle\frac{1}{\sqrt{2\pi}\sigma}$ that is needed to satisfy $\displaystyle\int G(x)\,\textrm{d}x = 1$. For higher dimensions this becomes a little bit more complex: $\displaystyle\frac{1}{\sqrt{(2\pi)^k |\Sigma|}}$, with $|\Sigma|$ being the determinant of the covariance matrix.

A somewhat numerical example

If you do not know the analytical formula, you could also just construct your multidimensional Gaussian without the prefactor, numerically integrate it, and use the resulting factor as the normalization constant.


The reason why you are not able to figure out whether these are two peaks is because your frequency plot has no "time element/information". It is the entire frequency content of the signal accross all time.
Take a short time Fourier transform of the signal. The time granularity of the STFT should be that of one heartbeat or close to it (any approx value of a typical human). If then you see two peaks in different time instances then these are two peaks at those times. maybe they are close or spread apart. That will help you determine whether it is a legitimate beat or not
A short time Fourier transform gives you a time contained frequency view of the signal.
How's it done is that you take your signal and divide it into small pieces, each piece length is what determines your time resolution. Then you take the FFT of each peice and plot it as a time-frequency plot.
You would need to know the sampling rate of the acquisition device. For ex: let's say it's 1Mhz, then every digital sample has 1us of information. Now if you want a granularity of 1ms then you need to take 1000 digital samples at one time and plot their FFTs. You can set ur time resolution based on a typical human heart beat or the granularity you want that will help you resolve the time occurrences of the two peaks.


Relationship between Discrete Deconvolution and Toeplitz Matrices

I think you are constructing the convolution matrix incorrectly. Let's do an easy example with some numbers to see how it is constructed

Let's imagine that youhave the following vectors of dimensions $M$ and $N$, And let's for the sake of it, say that $M = 5$ and $N = 4$

$\mathbf{a} = [a_{1},a_{2},a_{3},a_{4}]^{T}$

$\mathbf{b} = [b_{1},b_{2},b_{3},b_{4},b_{5}]^{T}$

$
\mathbf{c} = \mathbf{a}*\mathbf{b}
$

In matrix form, if we assume $\mathbf{a}$ is the system response, the convolution is given by

\begin{equation}
    c = \mathbf{A}\mathbf{b}
\end{equation}

with

\begin{equation}
    \mathbf{A} = \left[
    \begin{array}{ccccc}
        a_{1} & 0 & 0 & 0 & 0 \\
        a_{2} & a_{1} & 0 &  0 & 0 \\
        a_{3} & a_{2} & a_{1} &  0 & 0 \\
        a_{4} & a_{3} & a_{2} & a_{1} &  0 \\
        0 & a_{4} & a_{3} & a_{2} & a_{1} \\
        0 & 0 & a_{4} & a_{3} & a_{2} \\
        0 & 0 & 0 & a_{4} & a_{3} \\
        0 & 0 & 0 & 0 & a_{4} \\
    \end{array}\right]
\end{equation}

By doing the convolution with this matrix, you get the full operation including the transients. So the correct matlab code would be

A = toeplitz([a,zeros(1,M-1)],[a(1),zeros(1,M-1)])

(Check help toepliz to verify how it constructs the matrix)

Use this matrix and the result of the convolution is exactly as with the matrix multiplication.

Now, a word of wisdom, as you can clearly see, the full convolution matrix with the transients is not square, so in order to "invert" it and do deconvolution you will need to use the Moore-Penrose Pseudoinverse. And you will get some differences in the beginning and end of the resulting vector, but it is the best you can do, as using the pseudoinverse would be the same as performing least squares, since the deconvolution is given by

$\mathbf{b} = \left(\mathbf{A}^{T}\mathbf{A}\right)^{-1}\mathbf{A}^{T}\mathbf{c}$

This result will always hold, because the dimensions of $A$ are always $(M+N-1)\times M$, so the matrix is always "skinny" and the equation holds.

Try this, it should give you the result you were looking for.


Variance of periodogram estimate of the power spectrum

There are two sources of error in estimating a quantity: variance and bias.  In your question you have only considered variance.

The improvement in the spectrum estimate using the periodogram comes because there is less smearing of frequency components as you include more samples in your estimate of the power spectrum.  Using a longer data window causes less power from one frequency to be leaked into other frequencies.  However, this only reduces the bias of the estimate.

So, as you use more samples the bias is decreased, which means you do have a better estimate.  However, the variance remains unchanged, as stated in Numerical Recipes.  A discussion of this is given on this web page:

http://matlab.izmiran.ru/help/toolbox/signal/spectra7.html


One application of the Hilbert Transform is to obtain a so-called Analytic Signal. For signal $s(t)$, its Hilbert Transform $\hat{s}(t)$ is defined as a composition:

$$s_A(t)=s(t)+j\hat{s}(t) $$

The Analytic Signal that we obtain is complex valued, therefore we can express it in exponential notation:

$$s_A(t)=A(t)e^{j\psi(t)}$$

where:

$A(t)$ is the instantaneous amplitude (envelope)

$\psi(t)$ is the instantaneous phase.

So how are these helpful?

The instantaneous amplitude can be useful in many cases (it is widely used for finding the envelope of simple harmonic signals). Here is an example for an impulse response:

Secondly, based on the phase, we can calculate the instantaneous frequency:

$$f(t)=\dfrac{1}{2\pi}\dfrac{d\psi}{dt}(t)$$

Which is again helpful in many applications, such as frequency detection of a sweeping tone, rotating engines, etc.

Other examples of usage include:

Sampling of narrowband signals in telecommunications (mostly using Hilbert filters).
Medical imaging.
Array processing for Direction of Arrival.
System response analysis.


Convert fixed frequency bandwidth samples to time domain

As I suspect, it was an issue with insufficient number of points in freq domain data. I convert time domain data from the green trace to freq domain and got exactly the same values in freq domain as on the first picture with original freq data from network analyzer.


well, it depends on how similar the audio signal from the two tracks are.  i dunno if the two tracks are identical, to within some degree of error or noise, and to within some degree of time alignment, but  somehow, in a manner that is salient to your model or to your problem, you have to have a measure of something in the two signals to align.

so let

$$ a_x(t) = f\left\{ x(t) \right\} $$

$$ a_y(t) = f\left\{ y(t) \right\} $$

be that something that defines what it is you're trying to synchronize.   so you want to find a delay offset $\tau$, that maximizes the cross-correlation between $a_x(t)$ and $a_y(t)$

$$ R_{xy}(\tau) = \int_{-\infty}^{+\infty} a_x\left(t+\frac{\tau}{2}\right) a_y\left(t-\frac{\tau}{2}\right) w(t) \ dt $$

with whatever window function $w(t) \ge 0$ you like best.  find the $\tau$ that maximizes $R_{xy}(\tau)$

so i'm not telling you how to do this using discrete-time code.  you get to figure that out.


differences between DIT & DIF algorithms

The DFT is (let's say that $N$ is even):

$$ X[k] \ \triangleq  \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} $$

The radix-2 FFT works by splitting a size-$N$ DFT into two size-$\frac{N}{2}$ DFTs.  (Because the cost of a naive DFT is proportional to $N^2$, cutting the problem in half will cut this cost, maybe, in half.  Two size-$\frac{N}{2}$ DFTs appear to cost less than one size-$N$ DFT.

The Decimation-in-Time FFT splits the two DFTs into even and odd-indexed input samples:

$$\begin{align}
 X[k] \ &\triangleq  \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n\text{ even}} x[n] \, e^{-j 2 \pi \frac{nk}{N}} + \sum\limits_{n\text{ odd}} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m] \, e^{-j 2 \pi \frac{(2m)k}{N}} + \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m+1] \, e^{-j 2 \pi \frac{(2m+1)k}{N}} \\
\\
&=  \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m] \, e^{-j 2 \pi \frac{mk}{N/2}} + e^{-j 2 \pi \frac{k}{N}} \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m+1] \, e^{-j 2 \pi \frac{mk}{N/2}} \\
\\
\end{align}$$

This shows that if the evens are stuffed into one DFT and the odds are stuffed into another DFT, that the results of the two DFTs can be combined by multiplying the result of the right-hand (odd) DFT with a single factor $ e^{-j 2 \pi \frac{k}{N}} $ and added to the left-hand (even) DFT.  The result for the two frequencies $X[k]$ and $X[k+\tfrac{N}{2}]$ are the same except the right hand term is subtracted than added.

$$\begin{align}
X[k] &= \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m] \, e^{-j 2 \pi \frac{mk}{N/2}} + e^{-j 2 \pi \frac{k}{N}} \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m+1] \, e^{-j 2 \pi \frac{mk}{N/2}} \\
\\
X[k+\tfrac{N}{2}] &= \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m] \, e^{-j 2 \pi \frac{mk}{N/2}} - e^{-j 2 \pi \frac{k}{N}} \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m+1] \, e^{-j 2 \pi \frac{mk}{N/2}} \\
\end{align}$$

This shows how the butterflies are defined on the last pass of the FFT.

The Decimation-in-Frequency FFT splits the two DFTs into the first half and last half of the input samples:

$$\begin{align}
 X[k] \ &\triangleq  \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} + \sum\limits_{n=\frac{N}{2}}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} + \sum\limits_{n=0}^{\frac{N}{2}-1} x[n+\tfrac{N}{2}] \, e^{-j 2 \pi \frac{(n+N/2)k}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} +  e^{-j 2 \pi \frac{(N/2)k}{N}} \sum\limits_{n=0}^{\frac{N}{2}-1} x[n+\tfrac{N}{2}] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} +  e^{-j \pi k} \sum\limits_{n=0}^{\frac{N}{2}-1} x[n+\tfrac{N}{2}] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} +  (-1)^k \sum\limits_{n=0}^{\frac{N}{2}-1} x[n+\tfrac{N}{2}] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} \Big( x[n]  +  (-1)^k x[n+\tfrac{N}{2}] \Big) \, e^{-j 2 \pi \frac{nk}{N}} \\
\end{align}$$

For even $k$ we add the $x[n+\tfrac{N}{2}]$ to $x[n]$.  For odd $k$, we subtract.

$$\begin{align}
X[2m] &= \sum\limits_{n=0}^{\frac{N}{2}-1} \Big( x[n] + x[n+\tfrac{N}{2}] \Big) \, e^{-j 2 \pi \frac{n(2m)}{N}} \\
\\
 &= \sum\limits_{n=0}^{\frac{N}{2}-1} \Big( x[n] + x[n+\tfrac{N}{2}] \Big) \, e^{-j 2 \pi \frac{nm}{N/2}} \\
\\
X[2m+1] &= \sum\limits_{n=0}^{\frac{N}{2}-1} \Big( x[n] - x[n+\tfrac{N}{2}] \Big) \, e^{-j 2 \pi \frac{n(2m+1)}{N}} \\
\\
 &= \sum\limits_{n=0}^{\frac{N}{2}-1} \Big( x[n] - x[n+\tfrac{N}{2}] \Big) \, e^{-j 2 \pi \frac{n}{N}} e^{-j 2 \pi \frac{nm}{N/2}} \\
\end{align}$$

This shows how the butterflies are defined on the first pass of the FFT.


From the (limited) description the uHz rotator algorithm sounds like one of the phase-weighted averages from this site, but it's not an algorithm I am familiar with.

The Cramér–Rao lower bound$^1$ for estimating the frequency of sinusoid with amplitude $A$ in white noise with variance $\sigma^2$ is given by:

$$
\mathrm{var}(\hat{f}) \ge \frac{12}{(2\pi)^2\eta N(N^2-1)}
$$

where $N$ is the number of samples and $\eta = A^2/(2\sigma^2)$ (taken from Kay's Fundamentals of Statistical Signal Processing, Estimation Theory equation 3.41).

So, assuming $A =1$, $\sigma^2 = 0.1$, and assuming that frequency accuracy can exceed 8 significant digits means that $\mathrm{var}(\hat{f}) < (10^{-8})^2$ we get:
$$
\begin{array}
\ (10^{-8})^2 &\ge& \frac{12}{(2\pi)^2 5 N(N^2-1)}\\
\frac{20}{12} \pi^2 \times 10^{-16} &\ge& \frac{1}{N(N^2-1)}\\
N(N^2-1) &\ge& 6.079271 \times 10^{14}\\
N &\ge&  84713
\end{array}
$$

This means that your data length would need to be 84713 samples to achieve this accuracy.

For different values of $\sigma^2$ see the following table.

$$
\begin{array}
\ \sigma^2 & \eta & 10\log_{10}(\eta) & N\\
\\
0.1&    5.\ \ &    6.9897\ \ &    84714.\\
\\
    0.01&    50.\ \ &    16.9897\ \ &    39321.\\
\\
    0.001&    500.\ \ &    26.9897\ \ &    18252.\\
\\
    0.0001&    5000.\ \ &    36.9897\ \ &    8472.\\
\end{array}
$$

PS: I have emailed the author of that web-site and asked for their comment. Let's see if they respond.

$^1$ The Cramér–Rao lower bound is the lowest achievable variance of an unbiased estimator.


This is a homework type problem, so I'll give you a few hints to help you solve it yourself (and learn something while doing so).

From what is given, we can assume that the given difference equation describes a causal discrete-time system. Let $h[n]$ denote the impulse response. It must satisfy

$$h[n]=h[n-1]+\frac{1}{N}\big(\delta[n]-\delta[n-N]\big),\qquad n\ge 0\tag{1}$$

with $h[-1]=0$. Note that the second term on the right-hand side of $(1)$ only contributes to $h[n]$ for $n=0$ and $n=N$. Now start at $n=0$ and continue up to $n=N$ to see what happens. It will turn out that the impulse response has a finite length, i.e., the given system is an FIR filter.

The step response is the convolution of the impulse response with the unit step sequence, which simplifies to

$$a[n]=\sum_{k=0}^{\infty}h[k]u[n-k]=u[n]\sum_{k=0}^nh[k]\tag{2}$$

So the step response $a[n]$ is just the cumulative sum of the impulse response $h[n]$.


Where does homomorphic filtering stand in regards to DSP applications?

The cepstrum, or cepstral analysis, which is one form of homomorphic processing, has been used in music and speech analysis to estimate the deconvolution of source (exciter) pitch and vocal tract (and/or other) resonances; sometimes just used as a music/speech pitch detection/estimation method, but also useful for formant tracking as well.


what is the theory behind Butter function?

The trigonometric functions enter the calculation because of the use the bilinear transform for transforming an analog filter to a discrete-time filter. The bilinear transform warps the frequencies of the analog filter, that's why we have to pre-warp the frequencies of the analog filter in order to obtain the desired cut-off frequencies of the discrete-time filter.

The frequency warping of the bilinear transform is described by the following equation:

$$\Omega=\tan\left(\frac{\omega}{2}\right)\tag{1}$$

where $\Omega$ is the frequency (in radians) in the analog domain, and $\omega$ is the normalized frequency (in radians) in the discrete-time domain. Eq. $(1)$ is used to pre-warp the specified frequencies to obtain the edge frequencies for designing the analog prototype filter. After applying the bilinear transform to that analog prototype, the resulting discrete-time filter will have the correct edge frequencies.

The design procedure is as follows:

Pre-warp the given edge frequencies according to $(1)$.
Design an analog band pass filter with the pre-warped edge frequencies by applying a frequency transformation to a low pass prototype filter.
Apply the bilinear transform to the analog band pass filter to obtain the desired discrete-time filter.


just something to add, there are two kinds of canonical comb filters.  there are those made with only feedforward paths and those that have feedback.  in fact, you can have additional long delay elements of length $\frac{1}{f_0}$ where $f_0$ is the frequency of the spacing between the teeth of the comb.

in fact, a comb filter can be designed as generally as an FIR or IIR filter, but with the unit delay elements replaced by delay elements with delay $\frac{1}{f_0}$.  whatever frequency response you design for a low-pass filter (either FIR or IIR) can be used to design the shape of the teeth of the comb filter.

if you use a delay of half that above, a delay of $\frac{1}{2 f_0}$, your comb filter will be tuned to only the even-numbered harmonics.  in addition, you can instead design for a prototype high-pass filter (again with delay elements of half length $\frac{1}{2 f_0}$) and your comb filter will pass only the odd-numbered harmonics.

in continuous-time, using the Direct Form II, it's

$$ v(t) = x(t) - \sum\limits_{n=1}^N a_n v(t - n/f_0) $$

$$ y(t) = \sum\limits_{n=0}^N b_n v(t - n/f_0) $$

transfer function is

$$ H(s) = \frac{Y(s)}{X(s)} = \frac{\sum\limits_{n=0}^N b_n e^{-n s/f_0}}{1+ \sum\limits_{n=1}^N a_n e^{-n s/f_0}} $$

or frequency response

$$ H(j2 \pi f) = \frac{Y(j2 \pi f)}{X(j2 \pi f)} = \frac{\sum\limits_{n=0}^N b_n e^{-j n 2 \pi f/f_0}}{1+ \sum\limits_{n=1}^N a_n e^{-j n 2 \pi f/f_0}} $$

i gotta go now, but later we can discuss how you might design your LPF prototype that depicts the shape of the teeth.  and also, maybe we should discuss how to do a precision delay that has sub-sample precision so that $f_0$ can be precisely tuned.

in the other answer, the precision delay $D$ would be integer multiples of $\frac{1}{f_0}$ in seconds. and the normalized (to the sample period) precision delay (with integer and fractional part) is

$$ \frac{D}{T} = f_\text{s} D = i_D + f_D = \frac{n f_\text{s}}{f_0} \qquad n=1,2,... N $$

(don't conflate or confuse the symbols $f_\text{s}$, the sample rate, with $f_0$, the fundamental frequency or "pitch", or the fractional delay $f_D$.)


Pixels outside the image borders must be extrapolated.
Now, you need to chose the model of your extrapolation.

For instance, if you're working within the Discrete Fourier model a periodic extrapolation is the native choice.
Usually the most "Eye Friendly" solution is replication based on "Nearest Neighbor".

Your choice is as valid as any choice, you just need to figure which model it satisfies and if that model suits you.

Yet it has 1 noticeable caveat, it makes your kernel Spatially Variant, which means, without a caring treatment it will hurt performance.
The right thing to do is isolate the pixels which should be carefully treated because of the boundary condition and work very efficiently on the rest.

Tim Holy talks about it in the video Image Representation and Analysis (Go to 42:30).


You are right that the system $h[n] = 2^n.u[n]$ is unstable in absolute summability sense. And absolute summability of impulse response of an LTI system is the required condition for BIBO stability. Meaning Bounded Input bounded output stability. Here the input is bounded and it is only corresponding to 4 frequencies in digital frequency domain $\omega = -\pi /3, \pi /3, -\pi/2, \pi/2$. The magnitude of input cannot be greater than 3 at any value of $n$. But, the output will not be bounded because the system is not BIBO stable.

Even if you start the system at 0, meaning $y[n] = 0, \forall n < 0$, $y[n]$ will blow up as $n$ grows.
$$y[n] = \sum^{\infty}_{m=0}2^m.x[n-m], \forall n \ge 0$$
Putting $x[n] = cos[\frac{\pi n}{3}]$ as input:
$$y[n] = \sum^{\infty}_{m=0}2^m. \frac{e^{j\pi (n-m)/3} + e^{-j\pi (n-m)/3}}{2}$$
$$y[n] = \frac{e^{j\pi n/3}}{2}.\sum^{\infty}_{m=0}(2.e^{-j\pi/3})^m + \frac{e^{-j\pi n/3}}{2}.\sum^{\infty}_{m=0}(2.e^{j\pi/3})^m$$
Since both the sums are going to blow up as $n$ grows, we can say that $y[n]$ is not bounded.

Another way of saying this is that $H(e^{j\omega})$ does not exist for this system ,because DTFT sum does not converge. And so we cannot get $|H(e^{j\omega})|$ and $\angle{H(e^{j\omega})}$ for any $\omega \in [-\pi, \pi]$.

If we considered the system to be anti-causal, then the unit circle will fall inside ROC because ROC will become $|z|<2$. And the anti-causal but stable system's impulse response will be given by : $h[n] = 2^n.u[-n]$.

But this will be wrong thing to do, because the difference equation represents a Causal System, since the current output sample depends only on current input sample and a past output sample.


Both.  For baseband sampling (thus using only low-pass filtering), a sampling rate over over 2x the highest frequency is necessary.  For band-pass under-sampling, sampling at a bit over twice the bandwidth will work, as long as you remember the original band and don't fold (alias) the given band onto itself.  1000 Hz sampling will fold and invert a (slightly narrower than) 500 to 1000 Hz signal down to 500 to 0 Hz, and thus the signal won't be aliased with itself, but with an empty frequency band (after band-pass filtering stuff out, if needed).  Thus, no information will be lost (aliased), as long as you remember the original frequency band (e.g. it was not baseband).  Reconstruction may require a band-pass anti-aliasing filter, rather than a low-pass anti-aliasing filter.  Using knowledge of the original band, you can de-invert and upconvert the samples back to the original lower side-band of a 1000 Hz carrier, for example.


What do you mean by radio sensing a channel?

The station here is continuously monitoring the energy in the medium. Let us say you want to transmit in channel 1 of 2.4GHz band. Your RF hardware is always in receive mode if you don't want to transmit anything. You are continuously sampling the incoming signal in this particular channel. You use ADC to sample and compute the energy $E = \sum |x[n]|^2 $ of these samples corresponding to window of say, 100us. This is a simplified concept of sensing. You don't need to transmit anything to sense the channel. The value $E$ is computed in terms of $dBm$ unit.
After the above step, the $dBm$ value is compared against a threshold specified in standard. For example, let us say the minimum sensitivity needed for receiving QPSK modeulated data is -60dBm. You need to make sure that your transmission is above this. So you sense the channel. If the RSSI (that is $E$ computed above) is more than -60dBm, then you cannot transmit as the channel is BUSY. If it is less than -60dBm, you can transmit. But to be on safer side you can give a margin of say 5dBm. For higher modulation order (like 256QAM) the threshold will be lesser as you need to have very less ambient noise while transmitting 256QAM signals.


How to measure a detector's performance under extremely complex real-world environment?

on an ROC curve your plotting FP(x-axis) vs TP(y-axis) They are calculated by:

True-positives = # of correctly detected positives/# of actual positives

False-Positives = # of positives that are not true-positives/ number of positives

Now on your ROC curve, because you did not mention incuding this into your question above, you most likely did not include the "coin-flip" line. That is the line y = x (from x = 0-1)

This line represents how well your algorithm performs by just guessing randomly (Yes, No)

If your algorithm is a decent detector, the ROC curve will it will be above the coin-flip line, if it is below this line it is better at being an anti-detector, i.e., when there is no event and would be a poor detector for your application.

I've worked extensively on this for my master's thesis on onset detection. Below is an example of ROC curves plotted against the coin-flip line.


2-PAM or 4-PAM if I want to send data as fast as possible?

4-PAM transfers 2 bits ($log_2(4) = 2$) per symbol, while 2-PAM transfers 1 bit per symbol ($log_2(2) = 1$), so 4-PAM definitely transfers data faster.  The catch is that you need higher SNR to get the same bit-error rate (BER) that you would get with 2-PAM.


Effect of redundant training data in HMM-based speech recognizer/synthesizer?

I guess you confuse GMM and HMM trainings. Although in both cases EM algorithm is employed, Baum-Welch is used for HMM training.


A biquad is simply a two-pole IIR filter, where both the numerator and denominator of the transfer function are quadratic equations.  Thus, you can factor the coefficient equations, and then use any reference regarding a filter's behavior related to pole/zero placements in the Z-plane to determine the behavior.

One tutorial method I recall was a software package where you could drag the poles and zeros around and see a plot of the resulting response.

Q is a term from classic analog filter design which relates the inverse of bandwidth relative to center frequency.  So, given a known center frequency, the use of Q or bandwidth (depending on your exact interpretation or measurement method) is roughly interchangeable (after reciprocating one of them).


Your approach is putting the I and Q samples on a "sub-carrier" at frequency $f$. If that is really what you wanted to do, the typical demodulation approach is to multiply by the sub-carrier again and low pass filter to get I or Q. Dividing by cosine or sine creates time samples that go to infinity since you will divide by zero with those functions. The multiplication operation produces a baseband component and a component at twice the carrier frequency; this is made clearer with a quick refresh of the following trigonometric identities:
\begin{align}
\cos(\alpha)\cos(\beta) &= \frac{1}{2}\cos(\alpha-\beta)+\frac{1}{2}\cos(\alpha+\beta)\\
\sin(\alpha)\sin(\beta) &= \frac{1}{2}\cos(\alpha-\beta)-\frac{1}{2}\cos(\alpha+\beta)\\
\sin(\alpha)\cos(\beta) &= \frac{1}{2}\sin(\alpha-\beta)+\frac{1}{2}\sin(\alpha+\beta)
\end{align}
So in your case:
$$
\cos(2\pi ft)\cos(2\pi ft) = \frac{1}{2}\cos(0)+\frac{1}{2}\cos(2\pi ft) = \frac{1}{2} +\frac{1}{2}\cos(4\pi ft)
$$
The $\frac{1}{2}$ represents the amplitude of I and the $\frac{1}{2}\cos(4\pi ft)$ is the doubled frequency component that is removed with a low pass filter.
However I am not sure that you intended or need to have such a "sub-carrier"? If you were going to simulate QPSK modulation for example, you can modulate directly to your final carrier as in:
\begin{align}
I_t &= I \cos(2\pi f_c t), \quad\text{where}\quad I = {1,-1}\\
Q_t &= Q \sin(2\pi f_c t), \quad\text{where}\quad Q = {1,-1}
\end{align}
And further the I and Q samples are sent simultaneously by adding the modulated outputs $I_t$ and $Q_t$ above (hence the increased data rate; we send two bits in each "symbol).
To demodulate, assuming we have done the carrier recovery work to have a good estimate of $f_c$ with all phase offsets removed, we can multiply the received sequence by $2\cos(2\pi f_c t)$ and low pass filter to recover I, and multiply the received sequence by $2\sin(2\pi f_c t)$ and low pass filter to recover Q:
Noting the demodulation steps below where LPF() denotes removing the $2f_c$ frequency component with a low pass filter:
\begin{align}
I &= LPF\big\{2 I_t\cos(2\pi f_c t)\big\}\\
Q &= LPF\big\{2 Q_t\sin(2\pi f_c t)\big\}
\end{align}
The Q component goes to zero on the I output since the result is in the $\sin(\alpha)\cos(\beta)$ form resulting in $\sin(0) = 0$:
$$
0 = LPF\big\{2 Q_t\cos(2\pi f_c t)\big\}
$$
And similarly the I component goes to zero on the Q output:
$$
0 = LPF\big\{2 I_t\sin(2\pi f_c t)\big\}
$$
Therefore we can "multiplex" I and Q together given this ability to separate the two that operate simultaneously and independently in time.
Note I highly prefer to think and model using exponential frequency terms instead of sines and cosines - it is so much simpler! In fact this best explains implementations with I and Q paths as well, since to implement a complex signal you will need two real sequences $(e^{j\omega t} = \cos(\omega t) + j\sin(\omega t))$. Please see Frequency shifting of a quadrature mixed signal for more details on that.


IQ modulators, demodulators, and test equipment (oscilloscopes, etc.) can have channel imbalance in either or both gain and offset.  It looks like the Q channel of either the source or instrumentation changed or adjusted its offset during your measurement.


As Olli has said, Isotropic means the same in all directions. So the diffusion behaves the same regardless of which direction it is diffusing in.

Heat for example diffuses in an isotropic manner. (Well ... unless you have a strange environment with varying conductivity constants or something ...  but that's another story)

The heat equation is

$\frac{\partial u} {\partial t} = k \nabla ^2 u $

Where $U(x,y)$ is heat at point $(x,y)$ and $t$ is time

$k$ is the thermal diffusion constant.

I think it helps to see visually what this looks like

(From wikipedia)

This is actually just like doing Gaussian Blur on an image.

(Gaussian convolution is a solution to the heat equation)

Gaussian blur before and after

Note that the whole image gets the same blur

Anisotropic diffusion on the other hand may diffuse differently in different directions.

For example, perona and malik suggested doing diffusion based on the Laplacian (that's the 'edginess' $\nabla^2$) of that part of the image.

This can be used for performing blur only in certain parts of an image.

See here for an image which has had anisotropic diffusion. This has had the effect of removing noise, while preserving the main features of the image.

The equation proposed by Perona and Mailk is

$\frac{\partial U}{\partial t} = c(x,y,t) \nabla ^2 U + \nabla c . \nabla U$

where $c $ is a function such as

$C(||\nabla U||) = e^{-(||\nabla U||/K)^2}$ for example


How is the ROC of a transform function determined?

The ROC of $Y(z)=X(z)H(z)$ is the intersection of the ROCs of $X(z)$ and $H(z)$, unless there are pole-zero cancellations, which change the ROC of the product. In the given example we have to take into account the pole-zero cancellation at $z_0=2$. We can define a new function $X'(z)=X(z)(1-2z^{-1})$ with the pole at $z_0=2$ removed. We do not need to define a new $H'(z)$ with the corresponding zero removed, because that zero doesn't influence the ROC anymore (since the corresponding pole has already been removed). The ROC of $Y(z)$ is now the intersection of the ROCs of $X'(z)$ ($|z|>\frac13$) and $H(z)$. Consequently, the only possible ROC for $H(z)$ is $|z|>\frac23$. This means that the corresponding system is causal and stable, because the ROC is outside the pole with the largest magnitude and the ROC contains the unit circle.


High-pass filtering pre-processing before computing audio features

I agree that a large-order (12 taps+) high-pass Butterworth filter will be sufficient.

Why?  Well, Butterworth filters are commonly used in professional audio-editing applications to do filtering, because their frequency response is designed to be as flat as possible in the passband. This means that the filtered audio is of higher quality/fidelity, as there are less frequency "ripples" or "artifacts" in the passband than other comparable filter types (Chebyshev, elliptic). In other words, Butterworth has a more linear phase response in the pass-band.

The downside is that the Butterworth filter has a slower roll-off, and will require a higher-order to implement a particular stopband specification. But this shouldn't be an issue for you.

To answer your second question, this is why the above solution is preferable to using an FFT and zeroing out the bins: Why is it a bad idea to filter by zeroing out FFT bins?


Stabilizing an unstable system and removing oscillations from step response

The OP has a mistake in the numerator vector as there is an additional space between the minus sign and the 1.027e14 coefficient which will cause that to operate as an expression (-2.58e13 - 1.027e14) instead of negate the sign.
Had the system been correct, the resulting closed loop system with the compensator and gain chosen is underdamped. This can be resolved in this case by lowering the gain such that the closed loop poles are equidistant from the real and imaginary axis. Further details below.
The  resulting ringing is because the closed loop poles are too close to the $j\omega$ axis relative the the distance to the real axis, an underdamped system. The root locus shown by the OP shows all the possible closed loop pole locations depending on what loop gain is used, and thus as is with the current controller chosen the OP can lower the gain to select a condition with closed loop poles equidistant from both axis (which is an ideal trade of settling time and ringing). The general considerations for closed loop pole locations are summarized in the chart below.


If the length of the window should be $23$, $M$ must equal $(23-1)/2=11$:

M = 11;
m = -M:M;
w = .54 + .46*cos(pi*m/M);
plot(m,w)


Natural frequencies of a transfer function

The natural frequency is a real number, and Matlab computes it correctly by taking the magnitude of the (complex-valued) pole.

As an example, the denominator of the transfer function of a second-order continuous-time system is given by

$$D(s)=s^2+2\omega_n\zeta s+\omega_n^2\tag{1}$$

where $\zeta$ is the damping ratio, and $\omega_n$ is the natural frequency. The poles are given by the roots of $(1)$. Assuming an underdamped system ($\zeta<1)$, the poles are given by

$$s_{\infty}=-\omega_n\left(\zeta\pm j\sqrt{1-\zeta^2}\right)\tag{2}$$

from which we get

$$|s_{\infty}|^2=\omega_n^2\left[\zeta^2+(1-\zeta^2)\right]=\omega_n^2\tag{3}$$

So the natural frequency $\omega_n$ is given by the magnitude of the pole $s_{\infty}$.

Note that the natural frequency $\omega_n$ is not the actual oscillation frequency unless there is no damping (i.e., $\zeta=0$).


How to choose between taking the real part or the absolute value of an inverse discrete Fourier transform?

If your processed 2D FFT is conjugate symmetric, then any imaginary components in the complex IFFT result are just numerical noise or rounding errors.  Thus using them to compute an absolute value is mostly a waste of CPU cycles.


If the magnitude is $1$ and the phase were $0$, you would get

$$H(e^{j\omega})=1$$

which corresponds to $h[n]=\delta[n]$. So this is obviously wrong.

The fact that the phase of a negative real-valued number equals $\pm \pi$ can be easily seen by drawing the number as a vector in the complex plane and measuring its angle with the positive real axis. Analytically, the phase angle of a complex number must be computed using the atan2 function. The formula with the $\arctan$ function that you used only gives the correct result if the real part of the complex number is positive. Note that by dividing the imaginary part by the real part, sign information is lost: if the result is negative you don't know if the real part or the imaginary part was negative, and if the result is positive, real and imaginary part must have the same sign, but we don't know which. This is why the $\arctan$ function of the ratio of the imaginary part and the real part is generally not sufficient for determining the argument (apart from the fact that you can't compute the ratio if the real part is zero).

In sum, a negative real-valued number $a<0$ can be written as

$$a=|a|e^{\pm j\pi}$$

Also take a look at this related question.


the relationship between bandwidth and the datarate

Ideal symbol pulses are sequences of rectangular pulses.  Rectangular pulses, unfortunately, have horrible frequency profiles, as shown below.

The main lobe is where most of the pulse power is, and the other lobes are called side lobes.  Though they do diminish in power the farther out they go, they don't diminish very quickly.  The main purpose of pulse shaping is to dramatically reduce the power of those side lobes.

The width of the primary lobe is $\frac{2}{T}$, where $T$ is the rectangular pulse width, which is equivalent to the symbol period.  Thus, once pulse shaping is employed, the signal bandwidth is $\frac{2}{T}$.

Since BPSK transmits one bit per symbol, and the symbol rate is $\frac{1}{T}$, BPSK can transmit $\frac{1}{T}$ bits per second.  Therefore, the bit rate divided by the bandwidth is $\frac{\frac{1}{T}}{\frac{2}{T}}$, which reduces to $\frac{1}{2}$.  Thus, the bit rate is $\frac{1}{2}$ the bandwidth.

When the signal is at baseband half of the mainlobe is in the positive frequencies and half is in the negative frequencies.  Given that it is a real signal, the two halves mirror each other so some people consider it the bandwidth to be halved, or $\frac{1}{T}$.  Thus, at baseband the bit rate would be equal to the bandwidth.

QPSK transmits two bits per symbol, so the bit rate for QPSK is $\frac{2}{T}$.  It follows that QPSK can transmit $2$ bits per Hz of bandwidth at baseband, and $1$ bit per Hz at passband.


Scale FFT frequency range for a bars graph

To achieve that, use frequencies rather than bins. To convert from bin number to the frequency it represents use:

$$f = \frac{k}{N_{FFT}} \cdot Fs$$

Where:

$k$ is a bin number whose actual frequency you want to recover
$Fs$ is the sampling frequency
$N_{FFT}$ is the number of points of the FFT

You can work out your $N_{FFT}$ to ensure that your scanned frequency range will be divided into 16 segments. For example, to cover the frequency range $f_l=300 f_h=3000$ Hz with $N_{bars}=16$ bars and $Fs=44100 Hz$, you would have to use:

$$N_{FFT} = \frac{Fs \cdot N_{bars}}{|f_h-f_l|}$$

Where basically, you work out the spacing between the two frequencies $f_l, f_h$ for 16 segments and then "translate" it to the $Fs$ (but the ratio has been simplified there and that's why $Fs \cdot N_{bars}$)

This works out to be:

$$N_{FFT} = \frac{44100 \cdot 16}{|3000-300|} \approx 261.3333 $$

Which will have to be rounded to $262$ points.

So, now you know that your FFT has 262 points with a spacing of $\frac{Fs}{N_{FFT}}=\frac{44100}{262}\approx 168.3206$ Hz.

Your first bin is at $\frac{300}{168.3206} \approx 1.7823$ which will have to be rounded to $2$ and 3000 Hz should be at bin $2 + 16 = 18$.

Let's see, $\frac{18}{262}\cdot 44100 \approx 3029.7709 $ Hz.

Hope this helps.


How to measure algorithmically side lobes in real data

Zero pad your FFT (append zeros) such that you get more samples of the Discrete Time Fourier Transform (DTFT). The DTFT would give you the sidelobes and roll-off that you want to see (specifically for your diagram you would zero pad an input that is all ones, or an input that is any constant value).


The primary reason that root raised cosine filters are used is because they can be expressed in closed form and they are simple to understand. There are far better root Nyquist filters out there.

Usually they are optimized through least squares or remez exchange to improve stopband attenuation or some other channel specific criterion.

The ISI-free response (obtained through the cascade of a root Nyquist, the channel, and another root Nyquist) is only accurate in a flat channel (although it is also a good approximation in an almost flat channel). The shortcomings of root Nyquist filters have largely been addressed by DMT/OFDM, SC-FDM, and other proposed (and still experimental) technologies. These technologies result in ISI-free responses even for channels that the root Nyquist filters do not.

However, all of this is very dependent on channel models. If you really want to study the advantages and disadvantages then you need to provide a channel model, because the limitations and drawbacks of any method will be dependent on the channel model you are considering.


Time delay estimate and sampling frequency

A colleague of mine has investigated time delay estimation using compressed sensing in this paper: http://vbn.aau.dk/en/publications/compressive-parameter-estimation-for-sparse-translationinvariant-signals-using-polar-interpolation(56d35506-d4ae-4296-9a36-7ec95b1160c7).html
You can find the accompanying Matlab code here: http://www.sparsesampling.com/cpe/


What are some methods for detecting a sinusoid in noise?

Depending on the frequency of sampling, a FFT (Fast Fourier Transform) will work. For example, if your sampling rate is only once every 20 ms, you'll only get a few samples of the sinusoid, but if you're sampling every 0.5 ms then you'd get a lot more samples. FFTs usually work best with a large number of samples to average over. In that case even if your signal is the sum of a few sinusoids, they can be accurately determined.

Alternatively, you could have a look at the MUSIC algorithm. I'm not too sure about the details of how it is implemented, but it has been implemented in several real-time detection scenarios. An alternative to MUSIC is the Esprit algorithm.

Still, if your sample size is large enough and the noise doesn't swamp the signal entirely, then an FFT will be (on average) the fastest option.


Sequential Form of the Least Squares Estimator for Linear Least Squares Model

Slope from all samples obtained

To summarize the question's problem, you want to calculate the slope based on all samples obtained thus far, and as new samples are obtained, update the slope without going through all the samples again.

On the page you cite is the equation for calculation of the slope $m_n$ that together with $b_n$ minimizes the sum of squares of errors of the first $n$ of the $y$-values approximated by the linear function $f_n(x) = b_n + m_n x$. Converting the equation to zero-based array indexing for $x$ and $y$, we get:

$$m_n = \frac{\displaystyle\sum_{i=0}^{n-1}(x_i-\overline X)(y_i-\overline Y)}{\displaystyle\sum_{i=0}^{n-1}(x_i-\overline X)^2},\quad n \ge 2,\tag{1}$$

where $\overline X$ is the mean of the $x$ values and $\overline Y$ is the mean of the $y$ values. The problem with Eq. 1 is that the sums are actually nested sums:

$$m_n = \frac{\displaystyle\sum_{i=0}^{n-1}\left(x_i-\frac{\sum_{j=0}^{n-1} x_j}{n}\right)\left(y_i-\frac{\sum_{j=0}^{n-1} y_j}{n}\right)}{\displaystyle\sum_{i=0}^{n-1}\left(x_i-\frac{\sum_{j=0}^{n-1} x_j}{n}\right)^2}, \quad n \ge 2.\tag{2}$$

As $n$ is incremented, the inner sums can be updated recursively (iteratively), but by direct inspection the outer sum would need to be recalculated from the beginning. Least squares is an old and well-studied problem, so we don't try to bang our heads but look elsewhere. Converted to zero-based array indexing, Eqs. 16–21 of the Wolfram page on least squares fitting state:

$$\sum_{i=0}^{n-1}(x_i-\overline X)(y_i-\overline Y) = \left(\sum_{i=0}^{n-1} x_iy_i\right) - n\overline X\overline Y,\tag{3}$$

$$\sum_{i=0}^{n-1}(x_i-\overline X)^2 = \left(\sum_{i=0}^{n-1} x_i^2\right) - n\overline X^2.\tag{4}$$

These allow to rewrite Eq. 1 as (also given as Eq. 15 on the Wolfram page):

$$m_n = \frac{\left(\displaystyle\sum_{i=0}^{n-1} x_iy_i\right) - n\overline X\overline Y}{\left(\displaystyle\sum_{i=0}^{n-1} x_i^2\right) - n\overline X^2}, \quad n \ge 2.\tag{5}$$

It is easy to update the sums as $n$ is increased, and also updating $\overline X$ and $\overline Y$ is easy, by keeping and updating an accumulator containing the sum of all $x_i$ and another containing the sum of all $y_i$. The accumulated sum is then divided by the current $n$. The factor $n$ in the numerator and the denominator of Eq. 5 cancels together with one such division each, giving a formula for you to use:

$$m_n = \frac{\left(\displaystyle\sum_{i=0}^{n-1} x_iy_i\right) - \left(\displaystyle\sum_{i=0}^{n-1} x_i\right)\left(\displaystyle\sum_{i=0}^{n-1} y_i\right)/n}{\left(\displaystyle\sum_{i=0}^{n-1} x_i^2\right) - \left(\displaystyle\sum_{i=0}^{n-1} x_i\right)^2/n}, \quad n \ge 2.\tag{6}$$

Uniform sampling

If the samples are uniformly distributed, then further optimizations are possible. If the $x$ values are successive ascending integers, for example:

$$x_i = i,\tag{7}$$

then we get:

$$m_n = \frac{12\left(\displaystyle\sum_{i=0}^{n-1} i y_i\right) - 6(n - 1)\displaystyle\sum_{i=0}^{n-1} y_i}{n^3 - n}, \quad n \ge 2.\tag{8}$$

Avoiding large sums

Based on your comment to the question and to @M529's answer, perhaps this is not enough and you want to only ever store normalized sums that do not grow indefinitely with increasing $n$, similarly to what was done here. We can rewrite Eq. 8 as:

$$\begin{gather}m_n = 12 A_n - 6 B_n,\\
A_n = \frac{\displaystyle\sum_{i=0}^{n-1} i y_i}{n^3 - n},\quad
B_n = \frac{\displaystyle\sum_{i=0}^{n-1} y_i}{n^2 - n},\quad
n \ge 2,\end{gather}\tag{9}$$

which has recursions for $n \ge 3$:

$$\begin{gather}A_n = A_{n-1} + \frac{y_i}{n^2 + n} - \frac{3 A_{n-1}}{n + 1},\\
B_n = B_{n-1} + \frac{y_i}{n^2 - n} - \frac{2 B_{n-1}}{n},\end{gather}\tag{10}$$

that were found by simplifying the boxed parts of:

$$\begin{gather}A_n = A_{n-1} + \boxed{\frac{A_{n-1}\left((n - 1)^3 - (n - 1)\right) + (n - 1)y_{n-1}}{n^3 - n} - A_{n-1}},\\
B_n = B_{n-1} + \boxed{\frac{B_{n-1}\left((n - 1)^2 - (n - 1)\right) + y_{n-1}}{n^2 - n} - B_{n-1}},\end{gather}\tag{11}$$

where:

$$\begin{gather}A_{n-1}\left((n - 1)^3 - (n - 1)\right) = \sum_{i=0}^{n-2} i y_i,\\
B_{n-1}\left((n - 1)^2 - (n - 1)\right) = \sum_{i=0}^{n-2} y_i.\end{gather}\tag{12}$$

The recursive (iterative) algorithm can be tested in Python against the direct algorithm by:

from random import random
N = 10000  # total number of samples
n = 1
y = 2*random()-1  # Input y_0
sum1 = 0  # Direct (also needed to init recursive)
sum2 = y  # Direct (also needed to init recursive)
print("n="+str(n))
n = 2
y = 2*random()-1  # Input y_1
sum1 += y  # Direct (also needed to init recursive)
sum2 += y  # Direct (also needed to init recursive)
A = sum1/(n**3 - n)  # Recursive A_2
B = sum2/(n**2 - n)  # Recursive B_2
m = m_ref = 12*A - 6*B  # Recursive and direct m_2
print("n="+str(n)+", m=m_ref="+str(m_ref))
for n in range(3, N + 1):
    y = 2*random()-1  # Input y_2 ... y_(N-1)
    A += y/(n**2 + n) - 3*A/(n + 1)  # Recursive A_3 ... A_N
    B += y/(n**2 - n) - 2*B/n        # Recursive B_3 ... B_N
    m = 12*A - 6*B;                  # Recursive m_3 ... m_N
    sum1 += (n-1)*y  # Direct (not needed for recursive)
    sum2 += y        # Direct (not needed for recursive)
    m_ref = 12*sum1/(n**3 - n) - 6*sum2/(n**2 - n)  # Direct m_3 ... m_N
    print("n="+str(n)+", A="+str(A)+", B="+str(B)+", sum1="+str(sum1)+", sum2="+str(sum2)+", m="+str(m)+", m_ref="+str(m_ref))

The recursive algorithm seems stable, with the state variables A and B remaining small with a total of N=10000 independent random samples from a uniform distribution $-1 \le y < 1$. At the last sample the state variables and the output are:

n     = 10000
A     = 6.246471489186524e-07
B     = 5.20978725548091e-07
sum1  = 624647.1426721811
sum2  = 52.09266276755386
m     = 4.369893433735283e-06
m_ref = 4.369893433735271e-06

where m is $m_n$ calculated by Eq. 9 using the recursions of Eq. 10 for $n \ge 3$ and m_ref is the direct calculation of $m_n$ by Eq. 9. The two ways of calculating $m_n$ agree to 12 decimal digits. The direct sums sum1 and sum2 seem to be growing without bound. Of course, n is also growing without bound, and I don't think there's anything that can be done about that if this is what you want to calculate.

I did some testing against an arbitrary precision implementation using Python's mpmath (data not shown), and found that the root mean square error of the direct and the recursive methods did not differ significantly for $N = 10^5$ random samples. Personally, I would prefer the direct summation method using a large (in number of bits) fixed-point accumulator, if necessary. With direct summation, especially if it is not necessary to compute $m_n$ for all $n$, many slow division instructions are avoided, compared to the recursive method.

Recursion without auxiliary state variables

For what it's worth, I also found a recursion that doesn't require auxiliary state variables (such as $A$ and $B$), but requires remembering the previous two $m$ and the previous $y$:

$$m_n = \frac{2m_{n-1}(n - 2)}{n + 1} - \frac{m_{n-2}(n - 2)(n - 3)}{n^2 + n} + \frac{6y_{n-1}(n - 3)}{n^3 - n} - \frac{6y_{n-2}(n - 3)}{n^3-n}.\tag{13}$$

Slope from the most recent $N$ samples

Earlier I had thought that you want to calculate the slope based on $N$ most recent samples. I will keep the following stub of an analysis about that, in case it might be useful for someone.

For a least squares fit polynomial $f(x) = c_0 + c_1 x$ with $x=0$ at the center of the length-$N$ window, you can calculate the "intercept" $f(0) = c_0$ using a Savitzky–Golay filter, and the "slope" $f'(0) = c_1$ using a Savitzky–Golay differentiation filter, with both filters parameterized to have a polynomial degree of 1. These are linear time-invariant filters each described by an impulse response, easy to obtain for example using the function sgolay in Octave:

pkg load signal;
N=9;  # Number of samples in the window
sg0 = sgolay(1, N, 0)((N+1)/2,:);  # Savitzky-Golay filter
sg1 = sgolay(1, N, 1)((N+1)/2,:);  # Savitzky-Golay differentiation filter

The filters obtained have the following coefficients (which when reversed give the impulse response):

sg0 = 0.11111   0.11111   0.11111   0.11111   0.11111   0.11111   0.11111   0.11111   0.11111
sg1 = -6.6667e-02  -5.0000e-02  -3.3333e-02  -1.6667e-02   3.5832e-18   1.6667e-02   3.3333e-02   5.0000e-02   6.6667e-02

Figure 1. Savitzky–Golay (blue crosses) and Savitzky–Golay differentiation filter (red o's) coefficients, with $N=9$ and polynomial order 1. Swap the sign of the horizontal axis for impulse responses.

The first filter is a moving average. The second filter has a truncated polynomial impulse response. I think such filters can be implemented recursively (see my question about that), perhaps elaborated in this pay-walled paper: T.G. Campbell,  Y. Neuvo: "Predictive FIR filters with low computational complexity" IEEE Transactions on Circuits and Systems (Volume: 38, Issue: 9, Sep 1991).


Scilab Image and Video Processing toolbox provides imread() function.

The related How To entry gives example and it's similar to Matlab:

im=imread(filename);

See ATOMS article to learn how to install additional modules in Scilab.


The first thing I would enjoin you to do is to measure the actual accuracy of your ADC. The 12-bit performance announced by the manufacturer can often be reached only in ideal conditions, such as:

Using an external precision voltage reference rather than the built-in one.
Using a distinct power supply for the analog section of the MCU.
Using all the necessary EMI isolation/filtering parts as specified in the app notes... that your target dev. board might not have because of cost cutting (it can be as simple as a ferrite between VCC and AVCC).
Using slow sample rates.
Preventing any I/O pin to switch during acquisition (eg, disabling a SPI or UART spitting data in the background, disabling interrupts causing a GPIO to be toggled, disabling PWM peripherals..), or even putting the entire chip in sleep mode during acquisition.

In non-ideal conditions, I have seen ADCs lose up to 3 bit of precision. All this to say that your assumption that doing the processing with 12-bit data all along is correct as long as you truly have a resolution of 12-bit on your input; but if it's actually 8 bits because you are far away form the "sweet spot" of the ADC, there's no point keeping and processing 4 bits of noise all the way - you can just do the reduction at the source.

Converting from 12 bits to 8 bits through rounding or shifting is going to cause quantization noise correlated with your input signal. Dithering algorithms can be used to solve this problem and spread the quantization error. You also have the option of shifting the quantization noise in a frequency band where it is less harmful, through noise shaping methods.


How to design very narrow band, very sharp low pass filters - Only DC needed?

If all you want to do is get an estimate of the DC component of your data then you can simply obtain the average value of your signal over a unit of time. This can be easily shown by looking at the expression for calculating the Fourier coefficients for 'f=0' (DC component).

The easiest way to do this is via a moving average filter (http://en.wikipedia.org/wiki/Moving_average) which would work over a buffer.

In terms of an FIR filter, this looks like 'h = ones(1,Nbuff)./Nbuff' (so, 'h=[0.5,0.5]' for a 2 sample buffer). 'Nbuff' stands for the order of the filter. This can be applied directly on the buffer you obtain from your CODEC through a convolution function as provided by your processor's command set.

You could also use an 'online algorithm' (an IIR filter essentially) which maintains an updated estimate of the mean of the signal on a per-sample basis. For this purpose please see: http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm You could also design a "typical" IIR filter that can be as strict as you like, provided that your processor and sampling period can accommodate the computations required.

(Of course you can adapt between working in terms of frames and working in terms of individual samples)

Just a note, I have a feeling that you are struggling with the concept of Group Delay. The cut off frequency of your filter does not exactly mean that "new data can be read every 100 seconds". New data will be fed to your filter constantly, at the sampling rate of the system, but the phase delay caused on each harmonic component of your signal between the input and the output of the system might vary. How this varies, depends on the type of the system (whether it is FIR / IIR) and its form. Here is a starting point: http://en.wikipedia.org/wiki/Group_delay_and_phase_delay

Hope this helps.


A Delay Between Two Filtered Chaotic Signals

I assume that your experiment where you saw no effect was with signals that underwent a linear phase vs frequency process in their delay (constant group delay). Under many situations the phase of the channel in the frequency domain is non-linear, causing different frequency components to experience different delays (group delay variation). When you frequency select portions of such signals, you will then see the result of the delay that is dominant for those bands.

This is the same reason it is not best practice to simply use cross-correlation peak to evaluate the group delay. See further details specific to the challenges of doing this, and the better solution for those cases, at this posting here:

How determine the delay in my signal practically

And here:

How to reliably compute the group delay of a comb filter

At the bottom of the first post is a plot of group delay that illustrates my first point. This is an extreme case of delay variation where we see as duplicated below there are frequency selective regions where the delay is as much as 15 samples, while generally elsewhere it is only around 2-3 samples. If we further filtered this signal (with a linear phase filter to not cause further delay variation), we change the weighting between a 15 sample delay and a 2-3 sample delay.


Doing resampling, how can we choose lower cutoff frequency for LPF filter?

One thing to consider is the existing spectral content of the signal being downsampled.   If there is no significant energy above 15kHz, then the anti-alias filtering requirements are certainly relaxed since it won't have to do much for a 30kHz output sampling rate.

So where to put the filter cutoff frequency also depends on assumptions about the spectral content of the input signal as well as the transition band rolloff characteristics of the filter.  For a 30kHz sampling rate the 3dB cutoff point of the filter may need to be much lower than 15kHz to prevent aliasing from 15.5kHz or 16kHz (for example).

Much of this depends on your tolerance for the level of aliasing and at what frequencies.


When should we using moving average in algorithm design?

So when will you think you need moving average in an algorithm design?

If you mean moving average filters; moving averages as the name suggest are computed as averages of samples say,  $M-1$ previous samples (+ current sample) from input $x[n]$ to get an average output $y[n]$, repeating the process to get all $y$ samples. Computing their mean to get an averaged $y[n]$ signal. What you get is a smoothing effect.

As explained by Sanjit Mitra in the chapter on Discrete-Time Systems in his book. When your data is corrupted with noise and you happen to have multiple measurements of the same data samples, a good estimation of your signal is doing the ensemble average. But if you only have one set of measurement, then moving averages is a way of doing it. The result is a smoothing of your noisy signal. An M-point moving average is defined as a FIR system as:

$$
y[n] = \frac{1}{M}\sum_{k=0}^{M-1} x[n-k],\quad 0\leq n-k \leq M-1
$$

  And if we need one, how we should select the weight?

As seen in the equation, all the weights are equal and dependent on the length of the chosen filter:

$$
h[n] = \frac{1}{M}
$$

The choice of the length is critical. An efficient implementation however is done by introducing a recursive form and requires only 2 additions and 1 division.

Again, depending on your application, there are other variants  of the moving average that may be more suited to your signals.


It looks like there is a problem with your scaling.

The scaling for the DWT has the same interpretation as it does for the CWT.

For the DWT, the scale of the analyzing function (wavelet) is increased using a dyadic scale (increasing by factors of 2) in non-overlapping time intervals.  Increasing scale corresponds to narrower frequency distributions in the frequency domain.  Basically think of the analyzing function (wavelet) as a bandpass filter for which you are both narrowing the bandwidth and shifting down the center frequency by increasing the time scale.

As an example (see images below) I've transformed an audio signal using a CWT and then again using a DWT.  Notice that the results look similar, they scale in the same way, except in the DWT, the regions of analysis are discrete non-overlapping intervals (sharp boundaries).  Notice also that as the time scale increases (heading south on the CWT diagram), you have fewer DWT coefficients, so the number of blocks in the image decrease.  I don't see this kind of scaling in your data.


Objeckt tracking - filter training (MOSSE)

The solution to the problem was, to rearrange the quadrants of the filter. After rearranging the quadrants, my filter looks like this:

A lot more like the one in the paper. However, it still lacks in tracking precission.   To tweak this I played arround with the 2D gaussian. And obtained the filter below. Still not optimal but it's an improvement.

This is the formula for the 2D gaussian
$$
g_{i}(x,y) = e^{-\frac{(x-x_{i})^2+(y-y_{i})^2}{\sigma ^{2}}}
$$

In the paper I mentioned in the question, Bolme et al. used $\sigma = 2$. Instead of this I obtained $\sigma$ now dependent on the size of my template

$$
\sigma = \sqrt{{imgSize_{x} \cdot {imgSize_{y}}}} \cdot \frac{1}{16}
$$

like used by Galoogahi et al. (2015). Correlation Filters with Limited Boundaries


What Is the Purpose of Bundle Adjustment?

If all you want is to reconstruct a scene from a pair of images from a pair of calibrated stereo cameras, and your calibration is sufficiently accurate, then you do not need bundle adjustment.

You do need bundle adjustment if you want to reconstruct a scene from a sequence of images or for a sequence of stereo pairs, where the camera poses in each view are unknown. Search for "structure from motion" and "visual odometry".


What kind of result ? Are all the floating point data positive ? Anyway, I will list down some of the Quality Metrics and when to use them for which quality of data.

$1$. SNR - If the difference between the original signal and the reconstructed signal can be interpreted as a zero-mean noise process. SNR measurements relating to signal power or energy are less usual, because the result will be biased by the mean signal value. For audio and speech, zero-mean property applies anyway.

$$
SNR[dB] = 10 \log\frac{\mathcal{E}(s(\textbf{n})-m_{s})^{2}}{\mathcal{E}(s(\textbf{n})-\tilde s(\textbf{n}))^{2}}
$$

Sensitive against phase shifts ( whereas linear phase shift are uncritical for high visual quality ). For example, an observer would hardly find any disturbance in an image shifted by half a pixel using a high quality interpolator while the SNR criterion might indicate
a poor reconstruction quality.

$2$. PSNR - For image and video signals which are strictly positive and bounded within a maximum value $A_{max}$. PSNR is widely used for image and motion distortion comparison.

$$
PSNR[dB] = 10 \log\Sigma \frac{A_{max}^{2}}{(s(\textbf{n})-\tilde s(\textbf{n}))^{2}}
$$

$3$. Weighted SNR -  A more general form of SNR which additionally allows individual position-dependent weighting of errors is the weighted SNR.

$$
PSNR[dB] = 10 \log\Sigma \frac{r(\textbf{n})w(\textbf{n})}{w(\textbf{n})(s(\textbf{n})-\tilde s(\textbf{n}))^{2}}
$$

The weight $w(\textbf{n})$ could be used dependent on local structure (e.g. different weights for edge pixels or textured areas). It has value of one for relevant and
zero for irrelevant positions and allows comparison of two images over an arbitrary-shaped area.

Then there is,

$4$ Segmental SNR -

  SNR criteria are average measures, which do not take into account  fluctuations of quality. Often, a signal with only one highly distorted segment or region will be judged to be of bad quality, while the SNR computed over the entire signal may be extremely high. As possible workaround solutions, additional criteria such as SNR measurements over smaller segments, computation of minima, maxima or other criteria which  express fluctuations have
  been proposed. The  segmental SNR, which is the mean or expected value of SNR results from equal-size segments [Jayant, Noll 1984], penalizes fluctuations in SNR.

And also there is Perceptually Weighted Criteria.


Behaviour of microphone with high volume input

An overview of the analog signal path in a system such as yours:

Clipping can take place:

In the microphone as mechanical soft clipping, distorting the waveform rather than clipping exactly at some value.
In the built-in impedance converter field effect transistor (FET) of the electret microphone capsule. (I think this is very much like hard clipping)
In a separate preamplifier as analog hard clipping.
At the analog-to-digital converter of the audio interface, as hard digital clipping.
On the digital signal if any digital processing is done on integer data, as hard digital clipping.

Because the clipping is hard at numerical values -1 and +1, digital clipping must be taking place. Because it looks like there is hard clipping also at other numerical values than -1 and +1, there is probably hard clipping taking place at more than one stages separated by an alternating current (AC) coupling. Blowing directly into the microphone is quite an extreme disturbance, so there is probably also the 1st and the 2nd kind of clipping going on.

There are a few places in your system where there can be an AC coupling or equivalent:

In the microphone, as the pressure in the back of the omnidirectional closed-back capsule slowly equalizes with that of the outside world due to a (by design) leaky seal.
On the analog signal. AC coupling may be built in into the plug-in power supply to the microphone, into the preamp and into the audio interface.
Digitally in the audio interface or in the digital audio workstation (DAW). Check the settings!

If you blow into the microphone long enough then the AC coupling begins to see the average pressure level of that as direct current (DC), centering the signal to zero during the disturbance. When you stop blowing it takes some time for the coupling to adapt to the normal atmospheric pressure.

If you want the electronics to recover faster after a disturbance, try adding a capacitor in series with the signal output of the plug-in power supply as yet another stage of AC coupling. The capacitor should be small enough to make the time constant of the resulting resistor–capacitor (RC) filter short enough.


Book recommendations on Gabor filter and Fourier series

This is a book that I would recommend as a starting point: Richard Szeliski: Computer Vision: Algorithms and Applications


Audio EQ in by FFT?!? Please don't do that. CPU usage will be the least of your worries. http://blog.bjornroche.com/2012/08/why-eq-is-done-in-time-domain.html There was also a related discussion recently on the linkedin dsp group.

This post covers how to do a bell filter in the time domain, but it's extremely similar to create a low pass (actually a low pass is easier), and will be MUCH easier than trying to do an FFT (which won't do what you want anyway): http://blog.bjornroche.com/2012/08/basic-audio-eqs.html

You won't have to worry at all about memory (it only needs 2 memory slots) and CPU usage is very low. Plus latency is minimal. The final implementation is like 10 lines of code.

The original "cookbook" reference is here: http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

Now, that said, this filter may not do for your purpose, but you haven't really made it clear what that is, so please feel free to clarify.


Calculating Power of Sinusoidal Term and a Convolution expression

The power spectrum $P(e^{j\omega})$ is the Fourier transform of the autocorrelation function of $x(n)$:

$$R_x(k)=E\{x(n)x(n+k)\}\tag{1}$$

If you evaluate $(1)$ you should get

$$R_x(k)=\sigma_v^2\delta(k)+\frac{A^2}{2}\cos(k\omega_0)\tag{2}$$

Taking the Fourier transform of $(2)$ leads to the given expression for $P(e^{j\omega})$.

The convolution is obtained by noting that

$$u_0(\omega-\omega_0)\star W_B(e^{j\omega})=W_B(e^{j(\omega-\omega_0)})\tag{3}$$

and

$$\frac{1}{2\pi}\sigma_v^2\star W_B(e^{j\omega})=\sigma_v^2\frac{1}{2\pi}\int_{-\pi}^{\pi}W_B(e^{j\omega})d\omega=\sigma_v^2\tag{4}$$


3-tap FIR filter: simple expression for $H(e^{j\omega})$ using trigonometric identities

HINT:

Write $H(e^{j\omega})=e^{-j\omega}G(e^{j\omega})$ and use $(e^{j\omega}+e^{-j\omega})=2\cos(\omega)$. $G(e^{j\omega})$ should turn out to be real-valued.


How to compare the qualities of two interpolation (image resizing) algorithms?

Usually to test interpolation algorithm (for resizing, but also rotation and all sorts of "reversible" image transforms, you would apply the forward algorithm several times (ie resize to twice the size, take the result and resize it to 1.5x the size, etc .) and then apply all the backward transformations.

The more transforms you apply the more errors you will see and the more sensitive your comparison is. It can also represent your use case better (people pick an already resized image, resize it, and publish it)

Then you compute the max difference (or the standard deviation, or just visually compare) between the result and the original image.

In your case you can also apply all the forward transforms and then resize it once to the original size.