Is it fast time or slow time for the delay in the ambiguity function domain for radar signal?

Your question is a little bit unclear, but I'm going to try to answer according to my understanding of your question:

So it sounds like you're asking about pulse-Doppler radar processing, where we have two time dimensions: fast-time and slow-time.  Let us review those terms really quick.

In pulse-Doppler processing, multiple pulses are sent out by the radar at a steady rate/frequency (PRI/PRF, pulse repetition interval/frequency).  Each time a new pulse is sent out, the radar receives signals for a period of time, and then it switches back to a transmit operation and starts the cycle over again.  The data that is received is stored in a 2D matrix.  Each time the radar receives a new pulse, a new row is added to the matrix, and the data that is recorded at each receive interval is placed in this row.

This means that, like you said before, we have two time dimensions: fast-time and slow-time.  Fast-time corresponds to a small-scale time measurement of radar data at each pulse; slow-time is a (typically) much larger time scale, and corresponds to the dimension where each pulse was sent out.

Note that in this explanation so far, we have not taken the DFT of anything, so our 2D matrix represents time and time.  Like you alluded to before, if you take the DFT of the slow time dimension's IQ data, what you will get is a transformation to the Doppler domain. This happens because over the slow-time dimension, your target signal's phase will be changing if it is moving.  Assuming the target is moving at a constant rate, this means that the phase of the target signal will linearly change, and a liner change in phase is equivalent to oscillation at some frequency; this is  why we use the DFT for this: it allows for an easy way/convenient way to measure phase.

So now that we've reviewed that, lets talk about ambiguity functions:

Think of how range-Doppler data recording is structured.  There are multiple transmit/receive options, and in order to keep the target's phase stable, every transmit pulse is (typically) exactly the same.  Therefore, how can I tell whether a target signal appeared in range during my reception of the first pulse, or if the target appeared at any other range?  This is the concept of range ambiguity.  When you look at an ambiguity function, typically the x-axis is labeled as Tau, representing a time delay: this is equivalent to the range dimension.

Similarly, let's think about how Doppler is being measured: with pulse-Doppler processing, we're counting on being able to adequately measure the rate of change of the phase in question using a set sampling interval (our PRI/PRF from earlier).  It turns out that our PRF/2 is our effective Nyquist frequency for our slow-time sampling operation  Well, just like in sampling say, digital audio signals, we have to worry about aliasing any signals we record.  For audio signals, if my Nyquist frequency is, for example, 1kHz, any signals with frequency content above 1kHz will be aliased to look as if they are at another frequency.  This is exactly what can happen in the Doppler domain when using pulse-Doppler processing: if I have target velocity that is greater than my PRF, my DFT cannot accurately represent it, and that target signal "appears" somewhere else in the Doppler domain ambiguously.  This is called Doppler ambiguity.

So how does the range-Doppler map equate to the ambiguity function? Well, they've both effectively measuring the same thing!  The ambiguity function is a measurement of Doppler speed and range/time.  The range-Doppler map (RDM) from pulse-Doppler processing is a measurement of... Doppler speed and range/time!

There are two types of ambiguity functions: single pulse, and multi-pulse.  A multi-pulse ambiguity function simple shows you the "full" picture of both range/Doppler that you cannot see on the RDM due to measurement ambiguities.

Now a single pulse ambiguity function is also a 2D image, but we have no slow-time dimension.  Instead, our Doppler dimension is derived by using multiple hypotheses of the target speed and baking our guesses in to the pulse-compression operation that gets us our range dimension.

So hopefully all of that helps answer your question a bit.  If you'd like to learn more and have access to MATLAB, try walking through the examples on this MATLAB page


how to eliminate frame boundary discontinuity?

Windowing is the first line of defense.  you may need the windowing to be complementary, which is that all window functions (at their various frame delays) add to 1.  the Hann window is an example.  the Hamming window (which might have better sidelobe behavior than the Hann) is sorta complementary, but it adds to 1.08 (just a simple scaling issue).

if your window use is only for analysis (and not for synthesis) you might not need it to be complementary, so you could use a Kaiser window or a Gaussian window, both which have some useful mathematical properties in analysis.

the second line of defense is to use your 2048-sample as just buffers and space your frames (and the analysis windows that come with), by the detected period and adaptively advance the frame locations by a period each frame process.  that's a little weirder because you need to have, a priori, some knowledge of the period which what your pitch synchronicity is all about.

the most weirdest line of defense, is to advance each frame by a well-chosen "random" value between $L$ and $2L$ (where $2L$ is longer than the longest anticipated period) to scramble up whatever effect of the boundary discontinuity.

as we might say on comp.dsp, "you'll have to pay for more information".


Why does diagonal loading of a covariance matrix make an adaptive beamformer more robust in the case of a perturbed array?

I had worked on an array processing problem where I had used diagonal loading of the measurement covariance matrix. But I had used diagonal loading as a solution to what I thought was a numerical issue with my eigen values being too small. Since I had to invert the covariance matrix with small eigen values, I had numerical issues with the result.
Diagonal loading of the covariance matrix helped me get rid of the issue as the eigen values were now bounded by the diagonal loading value.

I am not sure, but one place where diagonal loading helps is in numerical stability of inverting the covariance matrix which is most likely encountered in adaptive beamforming.


Deconvolution Using Complex Division in The Frequency Domain

It's not about the index, it is about the Filter :-).
Think of your coefficients as filters and what would happen to a Sine Signal with frequency of $ \frac{N}{2} {F}_{s} $ that would be filtered (Convolved) with your samples of signal $ a $.
You will sum samples with the same absolute value yet one is negative and the other is positive, namely their sum is zero which means nothing will be left from this signal.
In practice we never use complex division for Deconvolution.
The minimum you can do is use some kind of regularization.
Zeros in Frequency Domain
In the noiseless case, having zero magnitude at some indices means those indices were zero either in the input signal or the filter.
The problem is we can never know hence when we want to recover any of them the data is lost. This is due to convolution being multiplication in frequency domain so at some index $ {C}_{k} = {A}_{k} {B}_{k} $ and if either $ {A}_{k} $ or $ {B}_{k} $ is zero then $ {C}_{k} $ is zero without being able to recover $ {A}_{k} $ or $ {B}_{k} $.
It means, that, potentially if there are no zeros for the data or the filter (At least in their overlapping domain) data can be restored.
In practice, this is very hard problem given noise.


Hough lines and Convex Hull methods give jagged lines

Let's go step by step.

First of all you can remove perspective distortion without camera calibration. e.g.

  Robust Radial Distortion from a Single Image, Faisal Bukhari and Matthew N. Dailey

  Robust Line Based Calibration of Radial Distortion from a Single View, Thorsten Thormahlen, Hellward Broszio, Ingolf Wassermann

For any given image, Hough transform might find multiple lines, many of those would be an indicator of the same edge. But luckily you could easily merge the lines which overlap. To detect the overlap you could well use the angle between the lines and the distance of lines. This way, you obtain a clustering of the lines. In any case, finding a short line on each edge would be sufficient for you to recover the projective quad.
If you are always sure that the orientation of this quad is more or less fronto parallel, then you could as well use directed gradients (X or Y) to emphasize the vertical or horizontal lines, during the edge detection step.

I agree with the other comments on Gaussian smoothing and anti aliasing. I would not recommend the use of convex hulls though, as they are very noise sensitive.


How may samples per symbol in ieee 802.11 g

If you oversample the signal by a factor of 2 as you've suggested, then you will have a bunch of empty subcarriers that are outside the specified signal bandwidth. When you modulate the OFDM frame to a time-domain signal, therefore, you'll have to use an IFFT that is twice as long. The result will be a discrete-time signal that is interpolated by a factor of 2 compared to using the nominal sample rate.


Reconstruct images from PCA reduced dimensions with NN

I'm not sure I understand your question. You want to create a neural network able to reconstruct an image thanks to its PCA ? This could be some sort of autoencoder that uses PCA for its latent representation. The encoder of the AE would simply be the PCA (AE isn't actually relevant here, nothing trainable in the "encoder", but I like the parallel since we're talking about generative machine learning/up-scaling 2D representations).

  If possible, what approach is the one you think would work?

Since an image will have 2D eigenvectors of already good shape, you could simply concatenate them, flatten, use a mean, use a concatenation of different 2D statistics/transformations, combine the results of several sub-networks starting from different 2D statistics of the eigenvectors (or the eigenvectors themselves), etc. as an intermediate representation. This is the tricky choice, following which you'll have to deal with typical AE decoder architecture design.

An illustration of the sub-networks idea (replace 1D inputs and 1D output with 2D data in your case):

  Is it possible to recontruct the images from the reduced dimensions?

If you're starting from reduced dimensions, it's not a problem: you can use either the so-called "deconvolution" (transposed convolution is the actual name) for direct 2D upscaling (see here for Keras implementation, and here for the equivalent PyTorch function), or simply if working with a flattened input/intermediate representation, re-shape a 1D large enough neural network layer, to go back to the original image shape.

Finally, note that if you use a neural network as it is commonly defined, your transformations will involve non-linearities.

That being said, not sure it's going to give you a great result, worth a try ! IMO you should try feeding each (or the n-largest) eigenvector into their own very simple sub-network, before combining their respective intermediate representations in a common final network that will output the right shape. Use simple architectures since you probably want to compare with linear combinations of eigenvectors.


2D Gaussian Seperation into 1D Gaussian components

Normalization is not "required". It only serves to have scale-consistent results, which a not so useful for visualization, but mostly for measurements: if the Gaussian kernel is "sum normalized", the result of the filtering of a constant image is the same constant image. Constant images are "invariant" for smoothing.
You need (or not) to do that for exactly the same reason as above. As it is separable, you can, if you want, normalize first on the 1D sum on the rows, second on the sum onn the columns.
They do not scale the 1D component: the 1D kernel g1d7 $[   0.006  ,  0.061    ,0.242 ,   0.383  ,  0.242  ,  0.061  ,  0.006]$ has almost $1$ average. They do a tensor product of $x$ ($x'. x$) to build a $7\times 7$ kernel, and only keep the center kernel with $273$ normalization.

So,  here is, in pseudo-code with Matlab syntax, a sequence of operations that leads to the $5\times 5$ array g2d5ci from Figure 3.

% 1D 7-tap Gaussian kernel
g1d7 = [0.006,0.061,.242,.383,.242,.061,0.006]';

% 2D 7-tap Gaussian kernel, normalized to 1
g2d7 = g1d7*g1d7';g2d7 =g2d7/sum(g2d7(:));

% 2D 5-tap Gaussian kernel cropped from the central 5x5 part of g2d7, normalized to 1
g2d5c = g2d7(2:6,2:6);g2d5c =g2d5c/sum(g2d5c(:));

% 2D 5-tap Gaussian kernel with integer values
g2d5ci = round(273*g2d5c);

The 1D kernel is not obtained by sampling the Gaussian, but by integrating it over  sub-intervals. For instance, instead of computing $g(x) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left({-\frac{x^2}{2\sigma^2}}\right)$ at $t=3$, you integrate it over the interval $[-3.5,-2.5]$.

The result is in the following figure (pretty close, isn't?):

The Matlab code goes like this:

increment = 0.001;sigma = 1;
time = [-3.5+increment/2:increment:3.5-increment/2]';
data = exp(-0.5*(time/sigma).^2)/(sigma*sqrt(2*pi));
plot(time,data)
data7 = reshape(data,length(data)/7,7);
data7p = sum(data7)';
data7p = data7p/sum(data7p);

figure(1);clf;hold on
plot(time,data)
plot([-3:1:3]',g1d7,'ob')
plot([-3:1:3]',data7p,'xr')
legend('Gaussian kernel','Given values','Recomputed values')
axis tight


Kalman Filter - Non Linear Measurement Model

So let's assume your model is:
$$
 \textbf{x}_{k} = \textbf{x}_{k-1} +  \textbf{n}_k
$$
where
$$
\textbf{n}_k =
\left(\begin{matrix}
n_x \\
n_y \\
n_z \\
\end{matrix}\right)
$$
and each of the $n_x,n_y,n_z$ are independent, zero-mean Gaussian noise terms with variances $\sigma^2_{nx},\sigma^2_{ny},\sigma^2_{nz}$.

The output equation is then:
$$
\textbf{z}_k = \textbf{f}(\textbf{x}_k) +  \textbf{m}_k = \left(\begin{matrix}
p_x^2 \\
p_y^2 \\
p_z^2 \\
\end{matrix}\right)
+  \textbf{m}_k
$$
where
$$
\textbf{m}_k =
\left(\begin{matrix}
m_x \\
m_y \\
m_z \\
\end{matrix}\right)
$$
and each of the $m_x,m_y,m_z$ are independent, zero-mean Gaussian noise terms with variances $\sigma^2_{mx},\sigma^2_{my},\sigma^2_{mz}$.

So following the EKF Wikipedia page, that means:
$$
\textbf{F}_{k-1} = \textbf{I}\\
\textbf{H}_{k} = {\tt diag}(2\textbf{x}_k)
$$
where ${\tt diag}$ creates a square matrix from the vector argument.

Implementing this in R (code below) shows the follow example of the true state and the state estimate.  As you can see, it has some problems because the output equation masks the sign of the state (that's why the top plot shows a mirror image around the zero axis).

R Code Below

#30103

# First, construct the signal model and generate the data.
T <- 1000

sigma_n <- 0.1
n_k <- array(rnorm(T*3,0, sigma_n), c(3,T))

sigma_m <- 0.1
m_k <- array(rnorm(T*3,0, sigma_m), c(3,T))

F <- array(c(1,0,0, 0, 1, 0, 0, 0,1),c(3,3))

x_0 <- array(1,c(3,1))
x <- array(0,c(3,T))
x[,1] <- x_0

z <- array(0,c(3,T))

for (k in 2:T)
{
  x[,k] <- F %*% x[,k-1] + n_k[,k]
  z[,k] <- `^`(x[,k],2) + m_k[,k]
}

# Next, form the EKF
H <- 2* array(c(1,0,0, 0, 1, 0, 0, 0,1),c(3,3))

Q <- sigma_m^2*array(c(1,0,0, 0, 1, 0, 0, 0,1),c(3,3))
R <- sigma_n^2*array(c(1,0,0, 0, 1, 0, 0, 0,1),c(3,3))

library("MASS") # For pseudo inverse ginv()

xkm1km1 <- matrix(rep(0,3*T+3),3,T+1)
xkm1km1[,1] <- x_0

xkkm1 <- matrix(rep(0,3*T),3,T)
K <- array(rep(0,3*3*T),c(3,3,T))
Pkm1km1 <- array(0,c(3,3,T+1))
Pkm1km1[,,1] <- array(c(1000,0,0 ,0,1000,0, 0,0,1000), c(3,3))

zhat <- matrix(rep(0,3*T),c(3,T))
err <-  matrix(rep(0,3*T),c(3,T))

for (k in 2:T)
{
  xkkm1[,k] <- F %*% xkm1km1[,k-1]
  Pkkm1 <- F %*% Pkm1km1[,,k-1] %*% t(F) + Q

  H <- 2*diag(xkkm1[,k])

  K[,,k] <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
  err[,k] <- z[,k] - H %*% xkkm1[,k]
  xkm1km1[,k] <- xkkm1[,k] + K[,,k] %*% err[,k]
  Pkm1km1[,,k] <- (matrix(c(1,0,0,0,1,0,0,0,1),3,3) - K[,,k] %*% H) %*% Pkkm1
  zhat[,k] <- as.numeric(H %*% xkkm1[,k])
}

p1 <- x
p2 <- xkm1km1
lims <- c(-5,5)
par(mfrow = c(3,1), pty="m")
plot(p1[1,], col="grey",ylim=lims)
lines(p2[1,], col="red")
title("True (grey) and Estimated (red) State p_x")
plot(p1[2,], col="grey",ylim=lims)
lines(p2[2,], col="red")
title("True (grey) and Estimated (red) State p_y")
plot(p1[3,], col="grey",ylim=lims)
lines(p2[3,], col="red")
title("True (grey) and Estimated (red) State p_z")


In short, yes: the key phrase that you want to search on is optimal smoothing.
You didn't ask how to do it, and the "how" is chapter-length, in a good book about Kalman filtering.  So I'm not going to try to describe that.
In general hand-wavy terms, though, optimal smoothing takes into account that if you know both past and future values of a system output, then you can do a better job of estimating its state than if you just know past values.
So, some close study of optimal smoothers should get you the "how".  As at least a partial benefit, you should find that the end effects that you see with filtfilt aren't as pronounced, because filtfilt uses a time-invariant filter run both directions, while an optimal smoother -- if done right -- adjusts its gains based on where in the vector it's doing its estimation.
Given that your system is nonlinear, you're going to have to do a bit of mathematical tap-dancing to make an extended (or unscented) optimal smoother, or you'll need to dig in the literature to find nonlinear optimal smoothing schemes.  I'm sure they're out there, they'll just be more interesting and obscure than optimal smoothers for linear systems.


Hydrophone kit to transmit (and receive) underwater

Yes, may be you are confused between Hydrophone and Modem.Hydrophone is just a passive device that records sound underwater. And what you just describe in your question is modem, which emit and receive sound waves underwater. there are definitely so many options you have. There are number of factors you need to consider before buying, mostly the operating frequency as most of the equipment in underwater works with KHz frequency. There are some amazing Software Defined Modems are available in the market but definitely those are expensive. If you want to just design and develop some communication protocols then you can use the UnetStack, basically convert your laptop sound card into a modem and have pretty good functionality. If you’re planning to buy some cheaper uw-modem, you can try buying popoto development kit, which includes transmitter and receiver. If you’re interested in some of the hydrophones then it worth give it a shot to Oceansonic icListen SC2, i must say this is an impressive military grade performance Hydrophone. This guys also had icTalker, but this can not generate custom signals, but definitely generate some chirp, continues sweeps, tones etc.
Also there are good market players in underwater domain such as,

Subnero
Benthos
EvoLogics GmbH
WHOI Micromodem
Popoto


FSK analysis: only ever find one frequency

You can make a very simple frequency discriminator by multiplying your received signal with a delayed version of itself, followed by a low pass filter. Adjust the delay to set the frequency range of your discriminator. This works because a multiplier followed by a low pass filter is a phase detector (multiply two sinusoids with same frequency and different phase and you will get the sum and the difference of the frequencies: the sum is the doubled frequency to filter out with the low pass filter and the difference is the phase). A change in phase versus a change in time is frequency, therefore the low pass filtered multiplication of the signal with a delayed version of itself will be proportional to the frequency over that particular interval. Choose your low pass filter and delay based on your symbol rate.

Optionally you can make a tuned discriminator with two pass band filters with your signal fed into the input of each one. The BW of the filter depends on your symbol rate. Follow each filter with a magnitude detector, and subtract one detector output from the other. The output of the subtractor will be your demodulated FSK.


Is the sampling period arbitrary for a bilinear transform, and why?

The sampling frequency is not arbitrary and important to specify for pre-warping. This is clear when reviewing how each analog frequency in the bilinear transform is mapped one-one to a digital frequency according to:
$$f = \frac{F_d}{F_s} =  \frac{1}{\pi}\tan^{-1} \bigg(\pi\frac{F_a}{F_s}\bigg) \tag{1} \label{1}$$
Where:
$f = F_d/F_s$ is the normalized digital frequency in units of cycles/sample,   $F_a$ is an analog frequency in Hz
$F_s$ is the sampling rate in Hz
Observe that $tan^{-1}(\infty) = \pi$ and thus the Bilinear Transform "warps" the entire analog frequency domain from $F=0$ to $F= \infty$ to the digital frequency domain of $f=0$ to $f=1$ (where $f=1$ is the sampling rate in units of normalized frequency), but the specific mapping of each analog frequency over that range is dependent on the value of $F_s$ that is used, as given by \ref{1}.
Normalizing the discrete-time frequency domain is commonly done such that the resulting spectrum associated with DC to the sampling rate extends from $f =0 \text{ to } 1$ cycles/sample, which may be the source of the OP's confusion.


Not an easy job. A loudspeaker produces a complicated 3 dimensional sound field, which heavily interacts with the environment. Even at a distance of a about 3 feet the reflected energy is significantly bigger than the direct energy at most frequencies. Typical loudspeaker measurements are

On axis response in an anechoic chamber
Full 3D response in an anechoic chamber
Radiated power in the far field in a highly reflective environment.
Room responses taken in "real world" setups, deriving a perceptually relevant spectral metric and averaging over a large number of environments and placements

Each of these measurements require specific measurement methods and techniques which are also not trivial.

Anything you can do with a single microphone in real room will most likely be either useless and misleading. The box is also not a good idea since any box has very strong internal acoustic modes and what you will be measuring will be dominated by these modes regardless of what the loudspeaker is doing.

Depending on the physical size of the loudspeaker and the size of your measurement space, you may be able to pull off a near field measurement with careful time gating of the impulse response (loosing a lot of low frequency resolution in the process). Getting this right is a lot of work and in my opinion not worth the bother.


Database into training and testing images

Some of the possible reasons would be:

Object recognition accuracy is not really well-defined. Instead accuracies are for particular datasets.
What you are working with is called cross-validation, where you split your own dataset into testing and training datasets and then use the trained model to make predictions on the test dataset.
Depending on how you do the cross-validation and the final model, the accuracy tests may vary. For example, if (even by a random process) your training set contains well-conforming images and the test dataset contains images which are not even closely conforming to your model, then the accuracy would be low.

As for a solution to your dilemma, you should probably consider looking at the dataset and sampling in such a way that this situation can be avoided. Alternatively you could use a more generalizable algorithm and check for overfitting of the data via cross-validation again.


Colorimetry: Gamma/gamut colors Color/color light

As you might know, there are different colour schemes/models available (RGB, HSL, HSV etc) and the important thing to notice is that these schemes are equivalent, which means you can achieve
exactly the same colour manipulations using any of the schemes. Why do we have different schemes then? It's because some things are (much) simpler to do using a scheme designed for a particular purpose.

The hue, saturation and lightness/luminance parameters are a part of the HSL colour scheme. The three parameters together define a colour. It's the luminance parameter that defines the lightness.
This particular scheme  makes it extremely easy to create effects such as colour gradients. The only parameter to change is,of course, the luminance.

As an example, let's say that we want to create a gradient from red to white. In the RGB colour scheme, red is encoded as (255,0,0) and white is (255,255,255). As you can see, both the green and blue components have to be changed
in order to create the gradient.

In contrast, in the HSL scheme red is encoded as (0,255,0) and white is (0,255, 255), which means that we change only the lightness to create our gradient.

To answer your question, all the three parameters constitute a colour whereas only the lighness/luminance defines a colour light.


If you have three or more antennas, simply don't put them in a linear array. Any other array setup doesn't suffer from the ambiguity. Linear array is the worst choice here.

Of course, that breaks the periodicity constraint on the linear array that makes the autocorrelation space method MUSIC So, you'd need to change that algorithm.

Luckily, you can simply put them in a circular array and change from cartesian to polar coordinates, and keep a very similar structure for your MUSIC implementation. Just remember that the resulting estimate you get might not be in angles, but in trigonometric function of angles, and needs to be converted back.

Another option is to use two different linear arrays and to combine the observation.

As long as you stick to the linear array, you simply never get any information that can tell one side from the other, and that's the end of the story :(


There are independent sub-carriers in both the positive and negative frequencies.  In other words, yes, they have completely different data.


In the end I used the Hilbert transform to compute the minimum-phase phase of my FIR.
I have used the following formula:

Where H is the spectrum of the desired filter.
Which has been used in:

Individualisation d’indices acoustiques pour la synthèse
binaurale
Phase Unwrapping for Spherical
Interpolation of Head-Related
Transfer Functions
Transformation characteristics of the external human ear

However I am not sure of the results I got.
A follow up question discussing the results and implementation can be found here: minimum-phase phase via Hilbert transform returned values


Bias instability in Gyroscopes : AVER / ADEV

According to this app note by Freescale: http://cache.freescale.com/files/sensors/doc/app_note/AN5087.pdf

which came from this recent related question: How to interpret Allan Deviation plot for gyroscope?

with a copied graphic below where they give the bias instability for each axis.  What a specific manufacturer does on their datasheet, you would be best asking them directly.

Yes you would want to see a low bias instability, and I would think that the same number at a lower $\tau$ would actually be worst unless you were only concerned with doing a measurement over that short duration in time (implying the turn-up occurs sooner so that unit would then be drifting sooner and over the longer duration would be a lot worst). You ultimately want the lowest value at the further averaging time indicating the best long term stability, but depending on the process and use, there may be reasons to be concerned with shorter term stability as well such as if you were making difference reading measurements at 1 second intervals (then the 1 sec ADEV would be of interest).


This is an idealized case of the echogram that simply depicts the times of arrival and energy of reflections. You can see this type of plot in literature and as an output in modelling software such as CATT-A and ODEON.

If you have a real-life impulse response then you can obtain information about reflections by representing the IR in logarithmic scale, i.e.:


Finding PHASE between two sinusoids with FFT/DFT (ArduinoFFT w/ C)

You can't do better than this for accuracy in my experience.

A Two Bin Solution

Select a frame size of 1 1/2 cycles or 2 1/2 cycles.
You only need to calculate two bins, either 1 and 2, or 2 and 3.  No need for a full DFT.  So it is very efficient as well.  By placing the frequency squarely near the center of the two bins you minimize noise effects and avoid needing to deal with near bin frequencies which call for a different set of equations for the fitting vector.
The alternative form is for when the frequency is really close to a bin.

DFT Bin Value Formulas for Pure Real Tones
An Alternative Form of the Pure Real Tone DFT Bin Value Formula

There are other tricks to do it point by point.  I'll let others provide those.


Particle filter, resampling step implementation

Your question is about the resampling step, and let's focus on that.

Resampling is used in particle filtering to counteract "sample impoverishment" that is the fact that some particles may have very low weights. This is a waste of resources as you want to describe your probability distribution ate best. You should note that it should not interfere (a priori) with the measurement of prediction steps.

The optimal situation is to have all particles having roughly uniform weights and that's the goal set by most resampling algorithms. I propose you use this function which uses something similar to histogram equalization for a ratio resampleof particles (a nd-array were the last line represents the weights):

def resampling(particles, resample):
    """
    Resample weighted particles.

    The particles are grouped in N_pop batches to represent multiple velocities.
    The probability is given for all particles but corresponds to the proba-
    bility of one batch.

    Input
    -----
    N particles as a particles.shape array

    Output
    ------
    N particles as a particles.shape array

    Parameters
    ----------
    resample : gives the ratio of particles that get resampled at every frame (between 0 and 1)

    """
    particles_out = particles.copy()
    N = particles.shape[1]  # number of particles

    # resample a percentage of the particles (the total number keeps to N)
    if resample > 0:
        # draw resample (in %) random addresses of particles that we will reassign
        N_resample = int(N*resample)
        address_resample = np.random.permutation(np.arange(N))[0:N_resample]

        # draw from this subset some addresses uniformly over their pdf using histogram equalization
        proba_resample = particles_out[-1, address_resample]
        proba_resample /= proba_resample.sum()
        address = np.interp(np.linspace(0, 1, N_resample, endpoint=False)+1/2./N_resample,#np.random.rand(N_resample),
                            np.concatenate(([0.], np.cumsum(proba_resample))),
                            np.arange(N_resample+1))
        address = [int(k) for k in address]

        # reassign these particles and set their weight to a uniform value
        particles_out[:, address_resample] = particles_out[:, address_resample[address]]
        subset_weight = proba_resample[address].sum() # should be \approx N_resample / N, that is, resample
        particles_out[4, address_resample] = subset_weight / N_resample #1/N_resample #particles_out[4, address_resample].sum()

        # ultimately, normalize weights
        particles_out[4, :] /= particles_out[4, :].sum()

    return particles_out


Computing the partial fraction expansion,

$$H (s) = \frac{2.302 s + 0.3548}{s^{3} + 0.739 s^{2} + 3.223 s + 0.3548} \approx \frac{2.2861 - 0.0309306 s}{s^2 + 0.626454 s + 3.1525} + \frac{0.0309306}{s + 0.112546}$$

For the time being, let us neglect the step response of the 1st order subsystem. Note that the 2nd order subsystem has the following form

$$\pm \gamma (s - z) \left(\frac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}\right)$$

where $\gamma \in \mathbb R$ is the gain (or attenuation) and $z \in \mathbb R$ is a (finite) zero. However, the formula you used to calculate the step response overshoot is only applicable for 2nd order systems of the form

$$G (s) := \frac{\omega_n^2}{s^2 + 2 \zeta \omega_n s + \omega_n^2}$$

whose zeros are at infinity. To summarize, the overshoot is smaller than you expected because:

you neglected the 1st order subsystem.
you neglected the attenuation $\gamma \neq 1$.
you did not subtract $\gamma$ times the derivative of the step response of $G (s)$ at the peak time.

In any case, not neglecting the above would only provide one estimate of the overshoot. To compute a satisfactory approximation of its exact value, you can take the inverse Laplace transform, differentiate and find where the derivative vanishes.


Difference between recursive and system with memory

you can have memory of past outputs.  thems would be $y[n-1], y[n-2]...$  and you can have memory of past inputs which would be $x[n-1], x[n-2]...$  and you know your current input $x[n]$.

if the output you're trying to compute, $y[n]$ will depend on any of the past outputs, it is recursive.  if not, it is not recursive.

if the output, $y[n]$, depends on any past outputs or past inputs, then it is not a "memoryless" system.


Effect of number of data segments in Periodogram-based PSD

If you increase the number of segments for a fixed total signal length, each segment becomes shorter. This means that you'll have a trade-off between the estimate's variance and the frequency resolution. Averaging many segments will reduce the variance but by using shorter windows, the ability to resolve closely spaced narrow band components is decreased. Shorter windows also decrease the peak amplitudes of any narrow band component. Shorter windows basically average out any details in the spectrum.


Removing white noise by taking the mean of many samples

The variance of the white noise will go down by $N$ through averaging (or standard deviation will go down by $1/\sqrt{N}$, where the standard deviation as a magnitude quantity will be visually consistent with the "spread" of noise on the graphic the OP is looking at. Thus to make it visually reduce by a factor of 10, the OP would need to average 100 times (if it is indeed white noise).
In practice with actual signals there is usually an upper limit on how much we can average, until non-stationary noise contributions (such as phase noise, 1/f noise etc) limit further improvement. If further averaging does not reduce this further, then such noise processes may be limiting further improvement.
Also be very careful in a simulation of white noise that the random number generator isn't actually creating the exact same sequence each time (in which there will be no change through averaging)- given the OP's plots I suspect that is what may be occurring here. This would be clear by inspecting each sample in successive runs. The other confirmation of a white noise process is to take compute the autocorrelation function  on the waveform, which for the white noise component will be an impulse at $\tau=0$ (or experimental approximation of one, limited by the actual time duration of the signal, so a large spike proportional to duration of the signal, and very small everywhere else).
One approach to remove white noise from a black box system with an unknown signal confined to a specific band of frequencies is to pass the output of the signal through a bandpass filter which performs a form of weighted averaging over previous samples of the one signal (rather than what the OP has done which is essentially averaging a system with N outputs, which by definition of ergodicity would be equivalent if the system is ergodic.)


Is there a counterpart for zero forcing in communications in DSP?

If you understand zero forcing from the perspective of vector spaces, you would be able to identify it different systems including DSP. In essence you are projecting your equivalent (received vector in communications) in the null space of a specific subspace which makes signals lying in that subspace to be zeroed. In communications, this would be the null space of interfering users (in a multi user scenario).

What's more, you can extend this idea from vector space to Hilbert spaces. For ex: consider the vector/Hilbert space of N degree polynomial, now you might want to project your polynomial on the null space  of the subspace spanned by some polynomials which are not of your interest.

For ex: you might want to approximate a N degree polynomial by a M degree polynomial, $M<N$, then you can find the null space of the polynomials from degree $M+1$ uptil N (using gram Schmidt procedure if possible) and then project the N degree polynomial along that to zero force those terms. So you see the idea is generic. This last example of polynomial approximation I quoted is also an example of zero forcing in DSP, less trivial one though.

Mbaz has quoted a reference for a more trivial explanation.


Does the cv2.waitKey argument determine the frame-rate while capturing a video?

It does not determine the framerate, but it has an influence on the framerate.

The first caveat is that although you're only waiting 1ms between frames, the program has other work to do in order to render each frame, so the delay between frames is not purely caused by that 1ms wait.

The next caveat is that (from the doc)

  Since the OS has a minimum time between switching threads, the function will not wait exactly delay ms, it will wait at least delay ms, depending on what else is running on your computer at that time.

So, you may be thinking, okay, I won't get 1/delay framerate, but I will get at most 1/delay framerate.  Alas!  Caveat 3 strikes:

What if you spam keyboard input, e.g. by holding down a key and letting it repeat?  Then the wait function will get cut short.  Again, you may be led astray by thinking, but if I use delay=1 (millisecond) this will be trivial because no keyboard at a rate anywhere near 1 KHz!  But remember that this delay is a minimum amount of time it will wait, in reality the wait is much longer.

In an application I wrote (more complex than this though) using cv2.waitKey(1) after showing each frame, I get a 50% increase in fps by simply holding down a key on my keyboard.

 Hands-off FPS: 64.09
Key-repeat FPS: 90.99

The moral is that if you want stable, controlled FPS, you should pick another route, because there are so many caveats and unpredictabilities here.


Different state-space representations for Auto-Regression and Kalman filter

I unfortunately don't know a whole lot about Kalman filters, but I think I can help you out with the state space stuff.

In Example 1, the AR model is exactly your good old DSP recursive definition of output:

$$ y_t = \alpha + \phi_1y_{t-1} + \phi_2y_{t-2} + \eta_t$$

In this case we write down the state-space model with direct correspondence with the above equation:

$$\begin{pmatrix}y_t \\ y_{t-1}\end{pmatrix} = \begin{pmatrix}\phi_1 & \phi_2 \\ 1 & 0\end{pmatrix}\begin{pmatrix}y_{t-1} \\ y_{t-2} \end{pmatrix} + \begin{pmatrix}\alpha \\ 0 \end{pmatrix} + \begin{pmatrix}1 \\ 0 \end{pmatrix}\eta_t$$

Note that in this case, the states of the system are current and previous values of the output.

In the second example, you're separating your states $c$ from your output values. This means that the states can now be anything, even though they still directly map onto output values. This way we get

$$y_t = \mu + c_t$$

$$ c_t = \phi_1c_{t-1} + \phi_2c_{t-2} + \eta_t$$

And therefore

$$\begin{pmatrix}c_t \\ c_{t-1}\end{pmatrix} = \begin{pmatrix}\phi_1 & \phi_2 \\ 1 & 0\end{pmatrix}\begin{pmatrix}c_{t-1} \\ c_{t-2} \end{pmatrix} +  \begin{pmatrix}1 \\ 0 \end{pmatrix}\eta_t$$

You should also recognize this as the standard state-space representation of a linear system, because you equations for state evolution and state-dependent output are two different equation. This separation is trivial in case of an AR model, but this latter notation is how we think of all linear state-space models in general.

The third example is a curious one. If you multiply out all the coefficients you will realize that it is actually equivalent to the first and the second examples. So why do it? I turns out that example 2 (being the proper state-space representation of the system) is called the Controllable Canonical Form of this system. If you do some reading or simply analyze the system carefully, you will realize that we can put this system into any state we like provided well-behaved values for $\phi_1$ and $\phi_2$ with the single input $\alpha$. Therefore we call such systems controllable, and it's very easy to see from this form of the state-space equations.

You should notice that two linear systems can be identical up to a change of basis. This means that we can pick a different basis to represent the same linear system. You can convince yourself that that's exactly what we've done to go from second to third example. Particularly, we like this linear transformation to transpose the state transition matrix, so that we would get for some unknown state $\boldsymbol{s}$

$$y_t = \begin{pmatrix}1 & 0\end{pmatrix} \boldsymbol{\alpha_t}$$

$$\boldsymbol{\alpha_t} = \begin{pmatrix}s_t \\ s_{t-1}\end{pmatrix} = \begin{pmatrix}\phi_1 & \phi_2 \\ 1 & 0\end{pmatrix}\begin{pmatrix}s_{t-1} \\ s_{t-2} \end{pmatrix} + \begin{pmatrix}\alpha \\ 0 \end{pmatrix} + \begin{pmatrix}1 \\ 0 \end{pmatrix}\eta_t$$

Now we can use the change of basis to find out what this state $\boldsymbol{s}$ has to be with respect to state $\boldsymbol{y}$. And we can calculate it to be

$$\begin{pmatrix}s_t \\ s_{t-1}\end{pmatrix} = \begin{pmatrix}y_t \\ \phi_2 y_{t-1}\end{pmatrix}$$

This form (transpose of Controllability Canonical Form) is called the Observability Canonical Form because if we can put a system in this form, we can easily deduce which states of the system can be observed by simply looking at the output. For some description of the canonical forms, you can read this document, and of course look around on the web. Note that in the document the states are flipped upside down, which does not change anything about the system representation, simply reordering the rows / columns of the matrices.


You don't really put any constraints on the family of codes (not even the alphabet). For real applications, you do need to put more constraints on what you want for a code. And distance may not be so important (e.g. the minimum distance of many LDPC codes is quite low, but you need to decode well beyond distance/2 to come close to capacity on many channels).

Brute force searching a code is a bad idea, since most codes are hard to decode (you need some structure to exploit to decode efficiently). Also, computing the minimum distance is NP hard in general (Vardy, 1997).

Here is an example of a code, which has worked well for many applications (CD's, RAID arrays, TV, Voyager mission, etc.). Reed-Solomon Codes are a family of maximum distance separable (MDS) codes, i.e. for a given length n and number of information symbols k, they have the largest distance possible (n-k+1).

They are also easily decodable, say by the Peterson-Gorenstein-Zierler algorithm or Berklenkamp-Massey algorithm (more efficient). You can apply these to generalized reed solomon codes as well (See R.M. Roth, Introduction to coding theory, for details). Many other codes are related to this family. For example, BCH codes are alternant codes from a RS code.


Most of the time, it's the ratio of a physical quantity to a reference. The reference can either be an industry standard (Example: dBSPL, dBm, dBu etc.) or some other suitable quantity (Example: voltage at the input of a filter).

"Physical quantities" in audio are things like pressure, volume velocity, particle velocity, displacement, voltage, current, force, etc. (20 log) or sometimes also power, energy, or intensity (10 log).

For example: a microphone measure "sound pressure" as a physical quantity and to make it logarithmic you use

$$P = 20 \cdot log_{10} (p/p_0)$$

where $p$ is the original sound pressure, $p_0$ a suitable reference pressure, (20 micro Pascals is industry standard) and P the resulting sound pressure level in dBSPL.

Taking a ratio is very important, since it also takes care of the units. In digital signal processing, you often just deal with numbers (without units), but a proper choice of reference is still important. Typically choices are either digital full scale, some suitably calibrated physical quantity or just simply "1" . The reference is important: a dB value without a well defined reference (implicit or explicit) is meaningless.


Verifying the computation of a convolution

I would have liked to give a longer answer, but yes, your approach is correct. I don't see any problems with what you laid out in your question.


Guitar Effects/Filters in real time using low cost DSP?

I'm sure you can implement a simple fuzz or overdrive effect on that board. Apply asymmetrical clipping to you your data, and filter the result with a biquad IIR filter. You'll need a lowpass filter to smooth out some of the nasty high frequencies after clipping (try a cut-off frequency $f_c\in [4..7]$kHz), and you might want a parametric EQ to boost some mid frequencies. You can compute the filter coefficients using the Audio EQ Cookbook. You can probably do much more than that but this should get you started.


Problems calculating 8-point FFT of an 8-point sine wave by hand

I will start from very beginning. So having signal: $ x[n] = [0; \ 0.7071; \ 1; \ 0.7071; \ 0; \ -0.7071; \ -1; \ -0.7071] $ and given exponent: $ w[n] = [1; \ 0.7071 - 0.7071i; \ -i; \ -0.7071 - 0.7071i;\ -1;\ -0.7071 + 0.7071i; \ i; \ 0.7071 + 0.7071i]$.

First we need to inverse the bits and rearrange our signal. For 8 samples we will obtain following change in indices: $ [0;\ 1;\ 2;\ 3;\ 4;\ 5;\ 6;\ 7] \Rightarrow [0;\ 4;\ 2;\ 6;\ 1;\ 5;\ 3;\ 7] $, thus our signal becomes: $ x[n] = [0; \ 0; \ 1; \ -1; \ 0.7071; \ -0.7071; \ 0.7071; \ -0.7071] $.

Now according to following schematic:

We calculate "butterflies" in following manner:

Where: $ W_{N}^{k} = e^{\dfrac{-2\pi i k}{N}} $ (rotating vector on our complex plane with step angle of $2 \pi/N $ )

Starting from a first butterfly in first stage ($x[0]$ and $x[4]$ as the inputs) we have following:

$ A = 0, B = 0, W_8^0=w[0]=1$, giving output: $ [0; 0] $, next butterfly:

$ A = 1, B = -1, W_8^0=w[0]=1$, giving output: $ [0; 2] $, next butterfly:

$ A = 0.7071, B = -0.7071, W_8^0=w[0]=1$, giving output: $ [0; 1.4142] $, next butterfly is the same with output:  $ [0; 1.4142] $.

So now, after first pass our vector is as follows (same to the one obtained by Peter Griffin):

$x'[n]=[0; \ 0; \ 0; \ 2; \ 0; \ 1.4142; \ 0; \ 1.4142] $

At second stage, first butterfly consists of samples: $x'[0]$ and $x'[2] $, where $ W_8^0=1 $. We can therefore calculate the output as: $[0+1\cdot 0; \ 0-1\cdot 0] = [0; \ 0] $. These are outputs $x''[0]$ and $x''[2]$ of the second stage.

Next butterfly is given by inputs $x'[1]$ and $x'[3] $, where $ W_8^2=w[2]=-i $. We calculate output $x''[1]$, $x''[3]$ as: $[0-i\cdot 2; \ 0-(-i)\cdot 2] = [-2i; \ 2i] $

Two following butterflies are calculated in the same manner:

$[0+1\cdot 0; \ 0-1\cdot 0] = [0; \ 0] \Rightarrow x''[4]$, $x''[6]$

$[1.4142-i\cdot 1.4142; \ 1.4142-(-i)\cdot 1.4142] = [1.4142-1.4142i; \ 1.4142+1.4142 i] \Rightarrow x''[4], \ x''[6]$

So our vector after second stage is following:

$x''[n] = [0; \ -2i; \ 0; \ 2i; \ 0; \ 1.4142-1.4142i; \ 0; \ 1.4142+1.4142i] $ - you've made a mistake here.

Finally, we calculate results for third stage. By looking at the main schematic, first butterfly consists of samples: $x''[0], x''[4] = [0; 0]$. The $W_N^k$ coefficient is given by $W_8^0=w[0]=1$.

Calculating butterfly we obtain:

$[0+1\cdot 0; \ 0+(-1)\cdot 0] = [0; \ 0] $ - these are final FFT values $X[0], X[4]$.

Moving to next butterfly based on samples $x''[1], x''[5] = [-2i; 1.4142-1.4142i]$ with coefficient $W_8^1=w[1]=0.7071 - 0.7071i$. Result: $[-2i+(0.7071 - 0.7071i)\cdot (1.4142-1.4142i); \ -2i-(0.7071 - 0.7071i)\cdot (1.4142-1.4142i)] = [-4i; \ 0] $

And next butterfly based on samples $x''[2], x''[6] = [0; \ 0]$, with  $W_8^2=w[2]=-i$ will give $[0+(-i)0; \ 0-(-i)0] = [0; \ 0] $

And the last butter fly based on samples $x''[3], x''[7] = [2i; \ 1.4142+1.4142i]$, with $W_8^3=w[2]=-0.7071 - 0.7071i$ will result in: $[2i+(-0.7071 - 0.7071i)\cdot (1.4142+1.4142i); \ 2i-(-0.7071 - 0.7071i)\cdot (1.4142+1.4142i)] = [0; \ 4i] $

So our final vector is following:

$X[k] = [0; \ -4i; \ 0; \ 0; \ 0; \ 0; \ 0; \ 4i] $

Which means MATLAB is calculating FFT correctly ;)

One final remark - always take care about $W_N^k$ factor - look on the schematic and it's pretty obvious.

Good luck!


This publication may be of interest for you:

Jonathan Driedger, Stefan Balke, Sebastian Ewert, and Meinard Müller: Template-Based Vibrato Analysis of Music Signals, in: Proceedings of the International Conference on Music Information Retrieval, New York, USA, 2016, pp. 239–245. [PDF]

Furthermore, there is an accompanying website for this paper.


Why images need to be padded before filtering in frequency domain

Multiplication in the frequency domain corresponds with circular convolution in the spatial domain.  This means that without padding the image properly, results from one side of the image will wrap around to the other side of the image.

You can think of 2D filtering as a sliding window that is centered over each pixel in the image and the center output pixel is a weighted sum of the pixels in the window.  With circular convolution, when the window hangs over the right edge of the image, it is actually wrapping back around to the left side of the image.  This means that output pixels on the right edge of the image will be affected by pixels on the left edge, which is almost never what is actually desired.

Zero-padding allows space for this wrap-around to occur without contaminating actual output pixels.


Good methods for hardware based 1-2 word speaker independent speech recognition?

There is really no need to do this in hardware. You can do the same on any cpu and it will be more flexible, energy efficient and updatable. You can pack simple computing core into your chip and run on it if you want. I don't think you'll gain anything from hardware implementation.

For a good implementation of keyword spotting you can check CMUSphinx http://cmusphinx.sourceforge.net, it has recently introduced good efficient keyword spotting. Fixed point arithmetics, noise suppression, tunable threshold. All that you need in single package.


How are these two discrete time signals equal?


Is the STFT a good solution to guitar monophonic pitch detection?

No.  An STFT is not only not overkill, it is not enough to find guitar "notes" reliably.  An STFT measures energy at selected spectral basis frequencies.  But a note usually refers to psychoacoustic perceived pitch, or pitch frequency, which is very often not the same as the spectral frequency peak, especially for guitar notes.  Guitars can produce more sound power at various (potentially slightly inharmonic) harmonic frequencies than at the perceived pitch frequency.  Also, the STFT basis vector frequencies may not accurately correspond to any of the musical frequencies of interest (being between them instead).

Try looking up pitch detection/estimation algorithms instead of using an STFT or bare FFT magnitude.


Calculate Min/Max value of multiple sine waves added together

The maximum value will occur when the frequencies all constructively interfere, which will certainly happen periodically if they are different frequencies.  Normalize by the sum of all the individual magnitudes.


How to construct a preemptive digital filter to neutralize an unwanted filter

If $S_x(\omega)$ is the desired input power spectral density (PSD), and $S_y(\omega)$ is the PSD at the output of the black box, which has the (unknown) frequency response $H(\omega)$, then the following relation holds:

$$S_y(\omega)=S_x(\omega)|H(\omega)|^2\tag{1}$$

Now you want to cascade the black box with a filter to recover the original PSD $S_x(\omega)$ from the distorted PSD $S_y(\omega)$:

$$S_y(\omega)|G(\omega)|^2=S_x(\omega)\tag{2}$$

From $(2)$ it immediately follows that

$$|G(\omega)|=\sqrt{\frac{S_x(\omega)}{S_y(\omega)}},\quad S_y(\omega)\neq 0\tag{3}$$

You need to design a filter with magnitude response $|G(\omega)|$ given by $(3)$. Note that for this application, the filter's phase response is irrelevant.


Note that a QPSK signal using a band-limited pulse has an envelope that passes through zero every time there is a phase transition of $\pi$. So its envelope is not constant if all symbol transitions are allowed.

Offset QPSK (OQPSK) doesn't have phase transitions of $\pi$. By staggering the $I$ and $Q$ signals by half a symbol interval, the maximum phase transition is $\pi/2$, leading to an approximately constant envelope. Slight droops in the envelope (occurring at phase transitions of $\pi/2$) can be eliminated by hard-limiting.

Note that with an appropriate choice of the transmit pulse, phase transitions in OQPSK can be completely avoided. This results in a modulation scheme with an exactly constant envelope. One example of such a continuous-phase modulation is minimum-shift keying (MSK).

For more information on OQPSK take a look at this question and its answers. Also browse this site for other questions on OQPSK.


Low analog lowpass + ADC vs High analog lowpass + ADC + low digital lowpass

I'm curious whether there is any difference between this strategy and, say, setting the lowpass to something really high like 250Hz and sampling at 500Hz, then once digitized and doing analysis, applying a digital lowpass at 30Hz

Good idea.  This is how oversampling ADCs work, sampling at a high rate using a poor quality converter (even a 1-bit converter), and then using digital filtering to downsample to the desired value with better accuracy.  Digital filters can have a much sharper cutoff than analog filters for essentially no cost, reducing aliasing that would otherwise get into your recording.  It can also improve your signal-to-noise if the analog signal isn't already restricting this.

  Would the double filtering induce artifacts?

Not unless you design the filters poorly.  :)


Separating excitation from filtering for sound sythesis

Does the above make any sense or where did I go wrong?

Yes.

  Is this idea of separating excitation from filtering widely used in sound synthesis? (I believe there is some audio compression algorithm, which uses LPC)

Yes it does.

For decades it was the strongest approach for speech - and by extension singing voice - synthesis.
As for musical applications, in the framework of subtractive synthesis, filters or resonator banks (Moog 914) were used to model the coloration brought by the body of an instrument, or to create vocal sounds.
The vocoder uses this principle to apply the formants (filters) extracted from a vocal signal onto a synthetic excitation signal coming from a synth.
It is not uncommon to find a filter simulating the body of an instrument as the final step in a physical model (for example a simple bowed instrument model).
Modal synthesis is another sound synthesis technique based on the representation of sound as an excitation shaped through a filter.
One could also consider that the virtual amp/speaker modelling techniques used to build today's electric guitar sound are inspired by this approach - the signal coming straight from the guitar being considered as a raw excitation shaped by the footprint of an amp & speaker.
Actually, any convolution processor (such as a reverberator) might be seen as an application of this technique - when you sing in a room, your voice is the excitation, the room the resonator.
This idea is also relevant to pitch-shifting. When a speaker speaks at a higher or lower pitch, the fundamental is modified, but the formants do not move (ideally...) - so a good pitch-shifting technique should preserve formants. To pitch-shift vocals or musical instruments sounds in a meaningful way, it might be a good idea to pitch-shift only the excitation part of the signal, while keeping the filter representation unchanged. In practice, though, we use techniques that do not explicitly make the actual separation (PSOLA in the time-domain, or spectral envelope manipulations in the frequency-domain when using a phase vocoder).

  Are there other methods which separate excitation from filtering?

Cepstral analysis. In some situations, adaptive filtering.


How exactly does the compression take place behind the simple ffmpeg commands?

MP3 is a perceptual coder. It creates a time/frequency representation and quantizes in the frequency domain primarily by calculating the spectral and temporal  masking thresholds of human auditory perception. It maximizes the masking of the quantization noise to minimize the audibility of the quantization noise.

Are spectrograms the best way to analyse the change in bitrates of an audio file?

No. Depends a bit on what you want to analyze, but a spectrogram doesn't correlate well with the audibility of compression artifacts. Most of problematic artifacts are temporal anyway, the loss of a few high frequencies is a comparatively smaller perceptual problem.

The STFT used, gives us the phase and the amplitude from which we take the magnitude and plot the spectrogram. So how different can I expect the phases of the higher and lower bitrates to be?

Phase will be fairly close for all spectral components that are above the masking threshold. Not so much for low energy components.

Is it possible to somehow recover the lost frequencies because ffmpeg based bitrate upscaling doesn't produce accurate results on the Spectrogram

No. If that would be useful in any way, the original inventors of MP3 (or AAC, DOlby, DTS, WMA, etc.) , would have included it in the codec. Many people have tried but it turns out the harm outweighs the benefit.


As I wrote in the comments, this can be solved using cross-correlation.

Take your mask of length $L$, invert it and start performing the convolution between this inverted mask and the original signal (or, directly, use scipy.signal.correlate() in Python with the two signals, the original one and the mask). If you get a peak in your correlation at sample $m$, this means that if you take your mask and overlay it over the original signal such that the first sample of the mask is at the $m$-th sample of the original signal (using the function with the parameter mode='same'), then you will have a good match (i.e., for that instant of time, the signal is pretty similar to the mask). Then you can conclude that from $m$ to $m+L−1$ the signal has a shape that looks like the one of the mask.

I checked the signal you've uploaded and, as you pointed out, following the procedure above one gets some non-desired peaks due to the fact that the original signal has some high values that don't correspond to a similarity in shape with the mask, but nevertheless return a high value of cross-correlation because the sum of the elementwise product is high. This could be solved "filtering out" these non-desired values.

The most similar signal to the mask will be... the mask. Thus, calculate the correlation between the mask and itself (the autocorrelation at $0$). This will give you a number that would represent the maximum possible similarity to the mask. Then, you can ignore those values in the original cross-correlation that are greater than the autocorrelation at $0$, as they surely correspond to those high peaks the original signal has. This should do the trick.


Unit delay operator as parameter in function

With $A(q^{-1},t)$ as given in the formula you get

$$A(1,t)=\prod_{k=1}^n(1-2\cos(k\hat\omega_0(t))+1)\tag{1}$$

and

$$A(\rho,t)=\prod_{k=1}^n(1-2\cos(k\hat\omega_0(t))\rho+\rho^2)\tag{2}$$

So $K(t)$ is simply the quotient of (1) and (2).


The Gaussian distribution is separable. Apply your transformation to each coordinate separately and you will get a 2D Gaussian.

If $G(x)$ is a 1D Gaussian, then $G(x) G(y)$ is a 2D Gaussian.

Thus, transform the distribution first according to one coordinate, creating a distribution that has normal distribution in one direction and uniform in the other. Then transform it again according to the other coordinate.


Ok, after some head scratching and re-re-re-reading I realize what I was missing. When creating the chirp signal used in convolution, the values of (k - n) which result in a negative need to be wrapped around and placed at the end of the signal rather than in the beginning.

So, replacing:

// Create y[n] by performing the following:

value = 1 - N

for n = 0 to n = 2N - 1

y[n] = e^i(PI / N * value * value)     // W^((k - n)^2 / 2)
value = value + 1

next n

With the following:

// Create y[n] by performing the following:

// Take care of values of k - n which are greater than zero.
for n = 0 to N - 1
y[n] = e^i(PI / N * n * n)

// Add the zero padding.
for n = N to M - N
y[n] = 0

// Now place the negative values at the end of the signal
for n = 1 - N to -1
y[M + n] = e^i(PI / N * n * n)

Yields the expected results.


Exponential moving average cut-off frequency

The discrete recurrence relation given on the linked page is
$y[n] = (1-\alpha)y[n-1] + \alpha x[n]$
with $x[n]$ being the input samples and $y[n]$ being the output samples.
So taking the Z-transform
$Y(z) = (1-\alpha)z^{-1}Y(z) + \alpha X(z)$
$H(z) = \dfrac{Y(z)}{X(z)}= \dfrac{\alpha}{1-(1-\alpha)z^{-1}}$
To find the -3 dB corner of the filter, in terms of the normalized digital radian frequency, $\Omega$ in the range $[0, \pi]$, one solves
$\left|H(z=e^{j\Omega_{3dB}})\right|^2 = \dfrac{1}{2} = \left|\dfrac{\alpha}{1-(1-\alpha)e^{-j\Omega_{3dB}}}\right|^2$
$ = \dfrac{\alpha^2}{\left|1-(1-\alpha)\cos(-\Omega_{3dB}) -j (1-\alpha)\sin(-\Omega_{3dB})\right|^2}$
which results in this equation
$\left[1-(1-\alpha)\cos(\Omega_{3dB})\right]^2+\left[(1-\alpha)\sin(\Omega_{3dB})\right]^2 = 2\alpha^2$
After some more algebra and trigonometry, I come up with
$\Omega_{3dB} = \cos^{-1}\left[1-\dfrac{\alpha^2}{2(1-\alpha)}\right]$
and since
$f_{3dB} = \dfrac{\Omega_{3dB}}{\pi} \cdot \dfrac{F_s}{2}$
I come up with
$f_{3dB} = \dfrac{F_s}{2\pi}\cos^{-1}\left[1-\dfrac{\alpha^2}{2(1-\alpha)}\right]$
Here's some Octave/MatLab code so you can verify the result on a plot; just change your Fs and alpha as required:
Fs = 10000.0
alpha = 0.01

b = [alpha];
a = [1 -(1-alpha)];

[H, W] = freqz(b, a, 1024, Fs);

f3db = Fs/(2*pi)*acos(1-(alpha^2)/(2*(1-alpha)))

plot(W, 20*log10(abs(H)), [f3db, f3db], [-40, 0]);
grid on

Update in response to comment:
The graph is a plot of the filter's magnitude response in dB (on the y-axis) vs. input frequency in Hz (on the x-axis).  The vertical line designates the -3 dB corner frequency.
Since you want to specify your -3 dB corner frequency to find the value for $\alpha$, let's start with this equation from above:
$\left[1-(1-\alpha)\cos(\Omega_{3dB})\right]^2+\left[(1-\alpha)\sin(\Omega_{3dB})\right]^2 = 2\alpha^2$
and after some algebra and trigonometry, one can get a quadratic equation in $\alpha$
$\alpha^2 +2(1- \cos(\Omega_{3dB}))\alpha - 2(1- \cos(\Omega_{3dB})) = 0$
which has solutions
$\alpha = \cos(\Omega_{3db}) - 1 \pm \sqrt{\cos^2(\Omega_{3dB}) -4\cos(\Omega_{3dB}) +3}$
of which only the $+$ answer from the $\pm$ can yield a positive answer for $\alpha$.
Using the above solution for $\alpha$ and the relation
$\Omega_{3dB} = \dfrac{\pi}{F_s/2}f_{3dB}$
one can find $\alpha$ given $f_{3dB}$ and $F_s$.
Here is an updated Octave script:
Fs = 40000
f3db = 1

format long;
omega3db = f3db * pi/(Fs/2)

alpha = cos(omega3db) - 1 + sqrt(cos(omega3db).^2 - 4*cos(omega3db) + 3)

b = [alpha];
a = [1 -(1-alpha)];

[H, W] = freqz(b, a, 32768, Fs);

figure(1);
plot(W, 20*log10(abs(H)), [f3db, f3db], [-40, 0]);
xlabel('Frequency (Hz)');
ylabel('Magnitude Response (dB)');
title('EWMA Filter Frequency Response');
grid on;

W2 = [0:75] * pi/(Fs/2); % 0 to 75 Hz
H2 = freqz(b, a, W2);
W2 = W2 / (pi/(Fs/2));

figure(2);
plot(W2, 20*log10(abs(H2)), [f3db, f3db], [-20, 0]);
xlabel('Frequency (Hz)');
ylabel('Magnitude Response (dB)');
title('EWMA Filter Frequency Response near DC');
grid on;

Additional information in response to Sep 11 '20 comment
For a stable, low pass filter, the pole, $(1-\alpha)$, must be inside the unit circle of the z-plane and be greater than $0$ on the real axis.  So we must satisfy
$$0 < 1 - \alpha < 1$$
thus
$$0 < \alpha < 1$$
These are bounds on $\alpha$ (and not a "definition") that come from requirements for stability and a filter passband centered at DC.
However, not all valid values of $\alpha$ will yield a -3 dB frequency that exists in the interval from DC to the Nyquist frequency, $\Omega \in [0, \pi]$ or equivalently $f \in \left[0, \dfrac{F_s}{2}\right]$.
When the pole, $(1-\alpha)$, is very close to the origin of the z-plane, the filter response will be so gradual or flat, that it will never drop to -3 dB or below.  Specifically, when
$$\Omega_{3dB} = \cos^{-1}\left[1-\dfrac{\alpha^2}{2(1-\alpha)}\right]$$
is undefined, due to
$$\left|1-\dfrac{\alpha^2}{2(1-\alpha)}\right| > 1$$
the filter response is so gradual, it will not have a drop down to -3 dB or below before the Nyquist frequency.  The $\alpha$ at which this happens is given by the positive solution to
$$ \alpha^2+4\alpha-4 =0$$
which is
$$\alpha = -2 + 2\sqrt{2} \approx 0.82843$$
EWMA filters with an $\alpha$ that is this high, have a pole so close to the origin of the z-plane, that the filter response is so gradual, there is no frequency attenuated by 3 dB or more.
The Wikipedia derivation of the 3 dB cutoff point is flawed and comes to an incorrect/inexact result.  Specifically it directly applies a continuous time, Laplace s-domain relationship:
$$f_c = \dfrac{1}{2\pi RC}  \quad \text{(Don't use this!)}$$
into a derivation of a result for a discrete time, z-domain, filter.  That step invalidates the derivation of an exact result in the Wikipedia article.
As a consequence, the result of the derivation in the Wikipedia article is only an approximation and (subjectively) only good for cutoff frequencies up to about $0.1 F_s$.


Try something like this:

$\begin{eqnarray}
\mu(n) &=& (1 - \alpha_1) \mu(n) + \alpha_1 x(n) \\
\bar{x}(n) &=& x(n) - \mu(n) \\
s(n) &=& (1 - \alpha_2) s(n) + \alpha_2 \bar{x}(n)^2 \\
\sigma(n) &=& \sqrt{s(n)} \\
\end{eqnarray}$

This is equivalent to sending your input signal to a 1-pole DC blocking high-pass filter with a pole at $\alpha_1$, then sending the result to a RMS detector whose response time is set by $\alpha_2$ - to measure how much the signal is wiggling once its central trend has been removed.

Observe that $\sigma(n)$ is consistent with the definition of the standard deviation.


Well you could try to estimate a transformation matrix that maps between the coordinates in the two images. The idea is that

\begin{equation}
\begin{bmatrix}
x_a\\
y_a
\end{bmatrix}
=
\begin{bmatrix}
A_{xx} & A_{yx}\\
A_{xy} & A_{yy}
\end{bmatrix}
\begin{bmatrix}
x_m\\
y_m
\end{bmatrix}
\end{equation}

where $x_a, y_a$ are the coordinates (in pixels for example) of the aerial image and $x_m, y_m$ are the coordinates of the photograph. Writing this out leads to the system of equations:

\begin{equation}
x_a = A_{xx}x_m + A_{yx}y_m
\end{equation}
\begin{equation}
y_a = A_{xy}x_m + A_{yy}y_m
\end{equation}

That way you can represent scaling, rotation and shearing operations between the images. You need at least four points per image to solve that system but I recommend more and then calculate the least squares solution using linear regression. This will also give you an indication on how well everything works. If there are large residual errors after the fit, then the linear transformation might not be sufficient.

If you then want to overlay the images, just apply the transformation matrix to your map and the result should be an image that has the reference coordinates on top of each other and the areas in between are linearly transformed accordingly.

There are many things that can go wrong here. For once, it only works for roughly flat terrain  (even though that effect becomes weaker if the plane is further away). Also, if there is lens distortions in the photograph, that will affect the accuracy of the fit (usually close to the image borders). And finally, map inaccuracies might throw you off.

In the end, it's just a matter of trying. If the 2D-transformation model is not capable of delivering accurate results, you need to model more effects. But estimating those will probably also require more referene points.

One other idea would be to put the image on a plane (like a texture) and place the plane in a 3D-space and match it with the plane of the photograph. That would end you up with a 3-dimensional transformation matrix for this plane. Might work as well.


Proof of complex conjugate symmetry property of DFT

Hint:

According to Euler's formula we have
$$e^{-j2\pi k}=\cos(2\pi k)-j\sin(2\pi k)=\ldots$$


Showing that the sum of zero-mean noise is zero. Then computing the convolution of zero-mean noise with a given function

Strictly: No such $B$ exists. You could simply have a sufficient streak of "bad luck" and draw positive $\eta>\epsilon>0$ continuously, for example. Obviously, $B<\lim_{N\to\infty}\sum_{n=-N}^{N-1}\epsilon_n<\lim_{N\to\infty}\left\lvert\sum_{n=-N}^{N-1}\eta_n\right\rvert\,\forall B\in \mathbb R$.
Granted, the event that every $\eta_n>0$ has probability 0 ($=\lim_{N\to\infty}\prod_{n=-N}^{N-1} P(\eta_n>0)=\lim_{N\to\infty}\left(\frac12\right)^N$). But you asked for the existence of such $B$, not for a proof of you being able to be sure that you'll be able to pick a finite $B$ which your sum-absolute never exceeds.
Sadly, there exist uncountably infinitely many cases for which the $B$ can't exist, and all have probability 0. Whether or not that makes a non-zero sum probability isn't trivial to say, far as I can tell.
What you could show is that there's a 100% probability that your sum is below some $B$; i.e. show that:
\begin{align}
1&=P\left(\lim_{N\to\infty}\left\lvert\sum_{n=-N}^{N-1}\eta_n\right\rvert<B\right)\\
&=P\left(\lim_{N\to\infty}\left\lvert\underbrace{\sum_{n=-N}^{N-1}\frac{\eta_n}{2N}}_{=:X_N\sim\mathcal N(0,\sigma_\eta^2)}\right\rvert<\frac{B}{2N} \right)\\
&=P\left(\lim_{N\to\infty}\left\lvert X_N \right\rvert<\frac{B}{2N} \right)\\
&=P\left(\lim_{N\to\infty}X_N <\frac{B}{2N} \right)-P\left(\lim_{N\to\infty}X_N <-\frac{B}{2N} \right)\\
&=F_{\mathcal N(0,\sigma_\eta^2)}\left(\lim_{N\to\infty}\frac{B}{2N}\right) - F_{\mathcal N(0,\sigma_\eta^2)}\left(\lim_{N\to\infty}-\frac{B}{2N}\right)\\
&=\frac12-\frac12\\
&=0\text,
\end{align}
which really isn't a probability of 1. Quite the opposite is true: You can be certain that if you let a random walk go on for infinity, that you're infinitely far from where you started.
By the way, the above proof doesn't use the fact that $\eta$ is uniform – it just "needs" that all $\eta_n$ are i.i.d. and that they have bounded variance.
You could construct the very same proof even without knowledge of the central limit theorem, just using Chebyshev's inequality.


Self study question about Padé Approximation and transfer functions

The following might help. Given the data of $5$ samples,

$$ x[n] = [1, 2, 3, 4, 5]$$

You will have possible choices for the orders $(p,q)$ of $a[k]$ and $b[k]$, such as $\{(4,0),(3,1),(2,2),(1,3),(0,4)\}$.

You can solve or each of the possibilities and check whether they yield stable systems or not, by looking at the roots of $a[k]$... afaik they should all be causal as this is Pade modeling assumption (unless equations modified to handle otherwise).

Note that for certain data sets (including this one) some tail coefficients might turn to be $0$ and the actual order can be less than that of indicated by $p$ or $q$.

Also note that a model with one pole is not an all-pole model. All-pole model requires that $q=0$, and does not depend on $p$.


Please refer to Parseval's theorem.
If your impulse response is $h(n)$ and its DFT is $H(k)$, then
$$\sum_{n=0}^{N-1} |h(n)|^2 = \frac{1}{N} \sum_{k=0}^{N-1} |H(k)|^2 $$
So the square root of the sum of the amplitudes in the frequency domain, divided by the number of points in the DFT should give you the RMS. It is common practice to calculate RMS on a frequency band basis to get more meaningful results.


Why does the separable filter reduce the cost of computing the operator?

Assume you have a $N\times M$ sized image.

If you know take what is classically used, a square filter kernel, of let's say size $L\times L$, you'd need to convolve that with the picture – which gives you $N\times M$ pixels, each needing $L^2$ multiply-accumulates. So you end up with $A_{2D}=L^2MN$ operations.

Now, if you can decompose that filter into an $L$-sized horizontal and an L-sized vertical 1D-filter, you could first do all rows – that's $M$ values per row, each needing $L$ operations, so $LMN$ for all rows – and then you'd do the same with the vertical filter, so $LNM$ for all columns – and you end up with $A_{1D}=2LMN$, and you'd only need to show that

$$\begin{align}
A_{1D} &< A_{2D}\\
\iff\\
2LMN &< L^2MN &&||:LMN, \text{ legal since $L,M,N >0$}\\
\iff\\
2 &< L
\end{align}$$

most filters are larger than 2.


How is the lowpass to bandpass transformation derived?

In the $s$-domain, the LPF-to-BPF transformation doubles the order of the filter.  that is because the LPF has one transition from passband to stopband, but the BPF has two such transitions.

Remember that when we do frequency response, we substitute $s=j\omega$.  That is, we evaluate the $s$-plane transfer function $H(s)$ as $H(j\omega)$.  The magnitude of $H(j\omega)$ or $|H(j\omega)|$ is the gain of the filter when the input frequencey is $\omega$.  Also keep in mind that the reciprocal of $j$ is its negative: $ \frac{1}{j} = - j $.

In the LPF, the middle of the passband is at DC, $\omega=0$.  The transformation to BPF is such that maps, using a nice continuous function, the middle of the BPF passband at $\omega=\omega_0$ to $\omega=0$.  A simple, mathematically continuous mapping as such is:

$$ \omega \leftarrow A ( \omega^2 - \omega_0^2 )  \quad A \ne 0 $$

So whenever the $\omega$ on the right (which is the BPF $\omega$) is equal to either $+\omega_0$ or $-\omega_0$, the $\omega$ on the left (which is the LPF $\omega$) is zero.

now suppose we set $A \triangleq \frac{Q}{\omega_0 \omega} $ which is still non-zero (and non-infinite for $\omega \approx \omega_0$).  now you have

$$ \omega \leftarrow Q \left( \frac{\omega}{\omega_0} - \frac{\omega_0}{\omega} \right) $$

$$ \begin{align}
j\omega & \leftarrow Q \left( j\frac{\omega}{\omega_0} - j\frac{\omega_0}{\omega} \right) \\
& = Q \left( \frac{j\omega}{\omega_0} + \frac{\omega_0}{j\omega} \right)  \\
\end{align} $$

or

$$ s \leftarrow Q \left( \frac{s}{\omega_0} + \frac{\omega_0}{s} \right)  $$

so this means, regarding frequency response:

$$ \begin{align}
H_\text{BPF}(j\omega) & = H_\text{LPF}\left(j Q \left( \frac{\omega}{\omega_0} - \frac{\omega_0}{\omega} \right) \right) \\
& = H_\text{LPF}\left(Q \left( \frac{j\omega}{\omega_0} + \frac{\omega_0}{j\omega} \right) \right) \\
\end{align} $$

or, regarding transfer function:

$$ H_\text{BPF}(s) = H_\text{LPF}\left(Q \left( \frac{s}{\omega_0} + \frac{\omega_0}{s} \right) \right) $$


Why a Large Diaphragm Microphone is limited to lesser frequency range

There are typically two reasons: First, a larger diaphragm is also heavier and hence it's inertia is bigger, which slows it down at higher frequencies (but also provides for significantly higher max SPL handling).

The main reason, however, is the relationship between wavelength and diaphragm size. The microphone basically "integrates" the sound pressure over the spatial area of the diaphragm. For a plane wave, the sound pressure is the same everywhere, but in a diffuse field that's not the case anymore. The wave length at 10 kHz is about 34 mm. So there are points on the diaphragm that are half a wavelength apart and it's entirely possible that one point sees positive pressure and the other negative pressure. These can cancel each other out (to some extent).


Prewitt operator and central difference?

Yes, the $1/2$ factor correction could be present, so that the magnitudes between a) the continuous derivative and b) the approximated gradient remain consistent in some way: a (continuous) line with a unit slope will have its discretized version get a unit gradient with the $1/2$ factor.

However, as long as the factor is nowhere applied, all gradient magnitudes remain comparable in a relative way. I suspect that, Prewitt being an ancient operator, computations were of concern. As is only involves simple additions and subtractions, one probably did not care about "costly" divisions, and most of all returning non-integer values.

The same happens with other kernels, see using a sobel operator - edge detection. The normalization is somewhat arbitrary, using a sum of absolute values, or in energy, as for Fourier transform normalizations.

However, if one uses different kernels, or different orders of derivations, it is safer to perform a consistent normalization, in order to be able to compare magnitudes.

The idea behind the constant columns is that a vertical edge should spread over several pixels, so you perform an additional averaging in the direction parallel to the edge, to reduce noise effects for instance, providing you with a more robust estimate of the actual gradient.


Assume that signal that comes to you through your antenna is:

$$ s_r[n] = I[n] \cos(\omega_0 n) + Q[n] \sin( \omega_0 n ) $$

which is a bandpass, modulated, narrow band signal where the components $I[n]$ and $Q[n]$ are the lowpass baseband in-phase and quadrature-phase components respectively.

Then in order to obtain the I and Q components from $s_r[n]$ you can apply the following just like the analog case:

$$I[n] = h_{LP}[n] \star ( s_r[n] \cos(\omega_0 n) ) $$
and
$$Q[n] = h_{LP}[n] \star ( s_r[n] \sin(\omega_0 n) ) $$

Where the lowpass filter $h_{LP}[n]$ is will have a cutoff frequency larger than the bandwidth of the signals I and Q and less than any adjacent band (if there is any). You may also apply a scaling such as $2x$ to adjust the magnitudes.

If you would like to form the analytic lowpass signal from the I and Q then you can easily form it like:

$$x_a[n] = I[n] + j Q[n]$$

where it can be shown that the analytic signal is also given by

$$ x_a[n] = A[n] e^{ j\phi(n) } $$

where $A[n] = \sqrt{ I[n]^2 + Q[n]^2 } $ and $\phi(n) = \tan^{-1}(\frac{Q[n]}{I[n]})$

Depending on the application, information in $A[n]$ (the envelope) or in $\phi(n)$ (the phase) is used.


First of all, the amplitude of a signal is expressed in arbitrary units, so there is no "right" answer. It's unlikely you'll be able to compute an absolute acoustic energy value - this depends on the gains of the recording chain, and the characteristics of the transducers.

Some audio analysis/feature extraction software consider that the full scale of the input audio signal is [-32768, 32767] (due to the predominance of the 16-bit format), some other implementation expect their input to be in the [-1.0, 1.0] range. Matlab's wavread function can return both formats when reading data, with the 'double' or 'native' argument.

I suspect the multiplication by 32767 in your code is simply there to keep the output consistent with a reference implementation of the algorithm which expected data to be in the [-32768, 32767] range, or maybe to keep the results consistent with a previous version of the software that used an input signal with the 'native' instead of 'double' scale.


the main idea behind BLIT is that these analog synth waveforms that we are trying to generate digitally can be thought of as the integral (over $t$) of impulse trains.

a sawtooth can be thought of as the integral of the sum of a little bit of DC and an impulse train.  a square wave is the integral of impulses of alternating signs.  the  triangle wave is the integral of the square wave.

so, to create bandlimited waveforms of the above, the impulse trains are bandlimited which means that each impulse $\delta(t-t_n)$ is replaced by a $\operatorname{sinc}(t-t_n)$ function, which is that impulse bandlimited through a Nyquist brick-wall LPF.  that sequence of bandlimited impulses is a BLIT.

then, since integration is a filter with s-plane transfer function of $H(s)=\frac{1}{s}$ and is LTI (Linear, Time-Invariant), integrating the BLITs will introduce no new frequency components.  if your BLITs are bandlimited, so are the other waveforms that are derived from filtering the BLITs.


IIRC, Karplus-Strong is usually a single linear harmonic waveguide.  Plucking a stringed instrument starts several things in vibration other that the string (sound board, air in cavity, etc.), many of them producing inharmonic overtones due to material stiffness, thickness and non-uniformity, some of which exchange energy and/of change in a non-linear fashion over time. The sharp peak high in the Cepstrum could indicate some kind of tiny but sudden frequency modulation(s) of the fundamental pitch.


Does speaker recognition system need the same sample rate audio?

Does it need the same sample rate for training and testing speaker recognition system?

I suppose that you are asking whether or not the difference in sampling frequency ($Fs$) will affect your classification result.

The short answer is "It depends on the data". The slightly longer answer is:

The $Fs$ has a direct impact on the captured bandwidth. At 8kHz, you are capturing (nominally) 0-4kHz. At 16kHz, you are capturing (nominally) 0-8kHz.

Speaker recognition is essentially trying to discriminate between the spectral differences in the sound of different speakers.

Therefore, a trivial example is one where if your $Fs$ is 8kHz and the key variation of spectral differences occurs between 6-7kHz, then your classifier will completely fail to discriminate between its classes.

Another detail is being sure that the $Fs$ is suitable, i.e. it includes the frequency range of interest, but, the recordings are awful and some other sound source is masking the range of frequencies the classifier depends on.

Therefore, "It depends on the data".

But, if we assume that the recordings are not problematic and a lower $Fs$ does indeed contain the frequency band of interest, then what might be worth considering more is the size of the transforms that are involved in the derivation of the MFCC that your classifier depends on. In other words, lowering the $Fs$ might require that you increase the length of the transforms to maintain the same (or similar) frequency resolution. Otherwise, adjacent frequency bins might get lumped together and reduce the clarity of the observed spectrum.

Finally, a more practical note would be to make sure that any existing code you might be using does not contain hard coded $Fs$. If the original author has assumed a fixed $Fs$ throughout the implementation of the model, it is likely that the classifier will have unpredictable behaviour if you feed it with sound files of different characteristics.

Hope this helps.


Direct path and reflections of ultrasound transmission

That sounds like a usual FMCW sonar.
It's technically identical to FMCW radar. So, multiply your receive with your transmit signal, you get a signal at the beat frequency, which is proportional to the range.
Standard FMCW approaches; there's a wealth of info out there.


What filter to use in audio analysis filterbank instead of FFT?

This has already been addressed in depth here:
https://stackoverflow.com/questions/5901483/simple-audio-filter-bank

and i found some c# code on the subject here:
https://waveletstudio.codeplex.com/

Copied answer from SO here

Using FFT to split an Audio signal into few bands is overkill.

What you need is one or two Linkwitz-Riley filters. These filters split a signal into a high and low frequency part.

A nice property of this filter is, that if you add the low and high frequency parts you get almost the original signal back. There will be a little bit of phase-shift but the ear will not be able to hear this.

If you need more than two bands you can chain the filters. For example if you want to separate the signal at 100 and 2000Hz it would in pseudo-code somewhat like this:

low  = linkwitz-riley-low (100, input-samples)
temp = linkwitz-riley-high (100, input-samples)

mids = linkwitz-riley-low (2000, temp)
highs = linkwitz-riley-high (2000, temp);

and so on..

After splitting the signal you can for example amplifiy the three output bands: low, mids and highs and later add them together to get your processed signal.

The filter sections itself can be implemented using IIR filters. A google search for "Linkwitz-Riley digital IIR" should give lots of good hits.

http://en.wikipedia.org/wiki/Linkwitz-Riley_filter


My answer is not from speech processing angle but from a generic signal processing involving Linear and Circular Convolution.

  Is zero-padding necessary for speech processing to satisfy linear
  convolution property?

If the FFT size is $N$, and the length of result of linear convolution size is $\gt N$, then you need to do zero-padding. Otherwise it will result in time-aliasing in time domain. The effect of this is, at the output of FFT you will get an under-sampled representation of the theoretical continuous time DFT signal.
$$
(Y.*H)(e^{j2\pi k/n}) = (Y.*H)(e^{j\omega})|_{\omega = 2\pi k/N}
$$
but
$$
(Y.*H)(e^{j\omega})|_{\omega = 2\pi k/N} \ne (\sum_0^{N-1}p[n]e^{-j\omega n})|_{\omega = 2\pi k/N}
$$
because of time aliasing. Here $p$ is the sequence obtained by linear convolution.

  What happens if I just stay with the circular convolution (i.e. NFFT
  length = signal length)?

Since your FFT length $N$ is equal to signal length $N$, the above time aliasing will happen.

  What if I just want to maintain the same signal length before and
  after processing?

Your FFT size has to be atleast equal to the resultant linear convolved signal length ($2N-1)$. Because that will be the length of signal you will be working with after the point-wise multiplication. So your FFT size has to be at least $2N-1$ before hand. See this question and its answers too Linear and Circular Convolution in Fourier Domain (DFT)


Confusion regarding frequency of Discrete Time Periodic Signal

Your assumption that the fundamental period has to be an integer since we are in discrete time is wrong. It can be any real number, nobody says, there has to be an integer number of periods in a signal of finite length. Also, if there is no $k$ representing this exact period, the peak in the spectrum will smear across neighboring bins.


Does the output of this system have a closed form expression?

You figured out the solution by yourself, so this is just about how to show it formally. First, split the sum into two sums:

$$y[n]=\sum_{k=-\infty}^{\infty}(-1)^ku[k]u[n-1-k]+\sum_{k=-\infty}^{\infty}(-1)^ku[k]u[n-2-k]\tag{1}$$

The unit steps give you upper and lower summation indices, but note that the resulting sums are only valid if the upper summation index is not smaller than the lower one:

$$y[n]=u[n-1]\sum_{k=0}^{n-1}(-1)^k+u[n-2]\sum_{k=0}^{n-2}(-1)^k\tag{2}$$

Now note that you can increase the upper summation index by any even number without changing the value of the sum. The second term on the right-hand side of $(2)$ can thus be rewritten as

$$u[n-2]\sum_{k=0}^{n-2}(-1)^k=u[n-1]\sum_{k=0}^{n}(-1)^k\tag{3}$$

In $(3)$ the argument of the unit step can be changed to $n-1$ because the value of the sum equals zero for $n=1$. Using $(3)$, Eq. $(2)$ can be rewritten as

$$\begin{align}y[n]&=u[n-1]\left(\sum_{k=0}^{n-1}(-1)^k+\sum_{k=0}^{n}(-1)^k\right)\\&=u[n-1]\left(2\sum_{k=0}^{n-1}(-1)^k+(-1)^n\right)\\&=u[n-1]\left((1+(-1)^{n+1})+(-1)^n\right)\\&=u[n-1]\tag{4}\end{align}$$


Cut-off frequencies for fractional sample rate adjustment

Don't be confused; you're doing everything correctly here:

$$\begin{align*}
f_{sample,in} &&\overset{\text{interpolate}}{\rightarrow}&& f_{sample,intermediate} &&\overset{\text{decimate}}{\rightarrow}&& f_{sample,out}\\
22 \mathrm{kHz}&&\overset{\uparrow 20}{\rightarrow}&& 440 \mathrm{kHz} &&\overset{\downarrow 11}{\rightarrow}&& f_{sample,out}
\end{align*}$$

To avoid images after interpolating, you will have to suppress all signal above 11kHz, which demands low pass with a cut off of $\frac1{2\cdot20}f_{sample,intermediate}$.
To avoid aliases during decimation, you will have to suppress all signal above 20kHz, which demands a low pass filter with a cut off frequency of $\frac1{2\cdot11}f_{sample,intermediate}$.

Since the filter from 1. implies 2., using one filter is sufficient.

This is a very commonly applied optimization.


Have you normalized the waveforms differently? The lowest point of the phase-shifted saw should be at amplitude -2 ln(2) / pi = -0.44 if the amplitude of the saw is normalized to -1 .. 1. The positive peak of the phase-shifted saw approaches infinity as you increase the bandwidth. So you could also have a clipping problem. Nonlinearities of the physical world affect the phase-shifted saw more strongly. The two waveforms have the same power spectrum and "complex envelope", but human hearing is more complex and may be able to detect some difference.


Compute dynamic range in a linear quantization system

Just to summarize the comments, for the next person who comes looking:

decibels are $20 \log_{10}{x}$ (or $10 \log_{10}{\left( x^2 \right) }$ when power levels used instead of the amplitude levels)

When looking at bit noise limited N-bit system, there are $2^{N}$ measurement levels, and the smallest distinguishable difference is $1$. So

DR$_{dB}$ = $20 \log_{10}{\frac{2^N}{1} } = N*20\log_{10}{2}\approx 6 N$


As Matt stated, you should use the correlation coefficient!.

Points 1 to 4 calculate the crosscorrelation. From that you have to find the highest peak (or lowest, if it has a higher absolute value). This value is the value of the nominator.

The denominator consists of the two autocorrelation values. Those are obtained by using the same algorithm where both signals are the equal. Here the peak should be in the middle (t=0) as already stated by welcomedungeon. Taking the square root of both autocorrelation values and multiplying them, gives the denominator.

Edit:
Maybe this description is more clear:

$\frac{max(abs(ifft(fft(x_1)*fft(x_2)')))}{sqrt(max(abs(ifft(fft(x_1)*fft(x_1)'))))*sqrt(max(abs(ifft(fft(x_2)*fft(x_2)'))))}$

The apostroph means conjugate complex.

Edit: Two examples with Matlab code:

Using the same signal:

x = rand(1000,1)-0.5;
max(abs(ifft(fft(x).*fft(x)')))/(sqrt(max(abs(ifft(fft(x).*fft(x)')))).*sqrt(max(abs(ifft(fft(x).*fft(x)')))));

gives 1;

Using a sine and a cosine, should also give 1 because they are delayed versions of each other:

x = sin([0:pi/100:10*pi]);
y = cos([0:pi/100:10*pi]);
max(abs(ifft(fft(x).*fft(y)')))/(sqrt(max(abs(ifft(fft(x).*fft(x)')))).*sqrt(max(abs(ifft(fft(y).*fft(y)')))))

gives approximately 1

Using windowing before transformation to frequency domain should improve results.


The answer is simple: if your system is linear, then a (regular) Kalman filter will do just fine. A very brief summary of the differences between the two:

The extended Kalman filter (EKF) is an extension that can be applied to nonlinear systems. The requirement of linear equations for the measurement and state-transition models is relaxed; instead, the models can be nonlinear and need only be differentiable.

The EKF works by transforming the nonlinear models at each time step into linearized systems of equations. In a single-variable model, you would do this using the current model value and its derivative; the generalization for multiple variables and equations is the Jacobian matrix. The linearized equations are then used in a similar manner to the standard Kalman filter.

As in many cases where you approximate a nonlinear system with a linear model, there are cases where the EKF will not perform well. If you have a bad initial guess of the underlying system's state, then you could get garbage out. In contrast to the standard Kalman filter for linear systems, the EKF is not proven to be optimal in any sense; it's merely an extension of the linear-system technique to a wider class of problems.


Why is there a dip in this Chebyshev filter response?

You are confused about which parameters are what.  The $S_{11}$ parameter is the power that's reflected back from the network; the $S_{21}$ parameter is the desirable bit that's passed forward.
Just by conservation of energy, if the filter is lossless then $\left | S_{11} \right | = 1 -  \left | S_{21} \right |$.  So that horrible deep null in $S_{11}$ actually just means that $S_{21}$ is at unity, which is desired.


Using Matlab:

You might try converting your image to greyscale first with

greycells = rgb2gray(rgbcells);

where rgbcells contains the color image of your cells (see here).

Then use Otsu's threshold on the greyscale image like so

threshcells = graythresh(greycells);

like so (see here).

EDIT
After OP explained what image should be the resulting image.

Using: FIJI (formerly ImageJ)

I am still not 100% sure as to which images you finally want to obtain. It seems like you only want to see the cells with the many purple dots (sorry, I am not a biologist myself) inside them and completely remove the others. For this attempt at a solution I use FIJI, but the principle can be easily applied in Python / Matlab, you name it...

Your images are somewhat different, because in the first image the other cells pretty much have the same color as the one in the middle (that you want to segment probably). In the second image the other cells are much lighter.

Let us look at the second image first:

Converting to HSI

Use the Color Transformer 2 Plugin to convert the images to HSI using the options "From Colorspace: sRGB" and "To Colorspace: HSI". Looking at the channels we can see that S and I carry pretty much the same information. So let's stick with the I channel. It looks like this: .

Selecting a Threshold

Now we use the useful "Auto Threshold" (Image->Adjust->Auto Threshold) feature of Fiji on this channel. Clip the image from the HSI stack and convert it to 16 bit first. Use the option "Try All", also check "Ignore White" and uncheck "white objects on black background". Now we can see a whole panel of different thresholding algorithms applied to this image. It looks like this:

You can see that many thresholding methods do what you want. Among them Otsus method (as you suggested) and MaxEnt. Now that we know MaxEnt or Otsus thresholding works for this image lets do the same thing again (meaning AutoThreshold), but instead of "Try All" we select "MaxEnt" as our threshold of choice (you'll see below why MaxEnt). Then we invert the threshold and overlay it on the image and we end up with:

Now for the first image

Here we can do exactly the same. Repeating the steps, we end up with an I channel like this:

The other cells are very well visible here. This is also true for the other channels. I have also plotted the histogram. The left hump corresponds to the wihte background, the right hump to pretty much all cells. This will be important in the next step.

Now try all thresholds again like above and we will see this:

You see that Otsus threshold separates the cells from the background, but does not give you the cell in the middle that you want. Look at the Histogram: due to the way Otsus method is designed (see here) it will seperate these two clusters. Which is not what you want (as I understand) in this case. But MaxEnt(Maximum Entropy) does. So use the MaxEnt threshold as above.

You will end up with this:

SUMMARY
This is my processing chain:

Do HSI transform
Apply MaxEnt threshold to the I channel
Use the binary mask on the RGB image again to clip out the relevant parts

OUTLOOK
If you are not satisfied with the results, then play around with different thresholds. Ultimately changing the colorspance and thresholding are powerful but very low tech ways of image segmentation. You might be interested in more sophisticated ways. For a software suite catered especially to cell biologists check out Cell Profiler. Former colleagoues have used it and were very happy with it.

LITERATURE
Maximum Entropy Threshold: Wong et al: "A gray-level threshold selection method based on maximum entropy principle",  https://ieeexplore.ieee.org/document/35351


Designing observer with non-observable system

If the unobservable dynamics is stable, you should be able to separate the unobservable states and the observable states via a similarity transform. You then design an observer for the observable part.


it seems like the bigger problem is that subsampling the signal would result in aliasing frequencies

No, subsampling in frequency domain corresponds to aliasing in time domain. So the idea here is to purposefully alias in time domain so that you get a sub-sampled FFT. That is exactly why 'Claim 3.7' of paper mentions $y_i=\sum_{j=0}^{n/B-1}x_{i+Bj}$. These are $n/B$ copies of $x[n]$ shifted by $B$ and overlapped. Each copy is positioned at $0,B,2B,\ldots,(\frac{n}{B}-1)B$. This aliased signal can be used to compute FFT ($B\log_2 B$ operations) to arrive at subsampled FFT. Total number of operations will the operations need to create the aliasing (summing of non-zero values of $x$ which equals $O(\text{supp}(x))$ in addition to the operation needed to compute FFT of aliased $x$ ($O(B\log_2 B)$).

For example, if a time domain signal is $1024$ samples long, ideally you need $1024$ point FFT. But if you sub sample by taking every 8th sample, $B=128$, $n/B=1024/128=8$, the time domain signal at $\tilde{x}[0]$ will be having addition from $x[0]$,$x[128]$,$x[256]$,$\ldots$,$x[896]$. Like this for $\tilde{x}[n]=\sum_0^{7}x[128k+n]$. You can now use this $\tilde{x}$ to compute $B$-point FFT to arrive at the Subsampled FFT.