How does window length influence cut-off frequency of filter and group delay?

Half average length. What is it?

Primarily something the author made up. They design a lowpass more or less as a weighted moving average filter. They simply use the window itself as the impulse response of the filter. That's NOT the same as designing an ideal low pass filter at the cutoff frequency you want and then windowing it down to the filter order. This seems to work reasonably well, you get frequency response of a nuttall window (which isn't bad) and the cut off frequency is roughly the sample rate divided by the window length.
Ir your sample rate is 44.1kHz and your window length is, say, N=32, the cut off will be around 1320 Hz or thereabouts.

Compensation of the delay. .... What is wrong with my interpretation of this formula?

Your interpretation is correct. If the window length is even, you can't correct for the group delay by simply shifting samples. The code only "approximately" compensates for group delay and is off by half a sample.
There are a few more things wrong about the code:

If someone calls it with half_average > fft_size/4 , bad things will happen
Multiplication in the frequency domain is circular convolution (not linear convolution) in the time domain. The input must be properly zero-padded before FFT to avoid time domain aliasing.


Correlation between frequency and initial phase estimates using FFT

It's a spurious artifact due to not referencing your phase generation to the center of your data window, and not doing an fftshift before the fft to place your phase reference point at fft input index 0.

e.g. For varying frequencies, either you have to know, generate, or want the phase at location T/2 of your continuous sinusoid.  The phase result of an FFT after an fftshift will produce an estimate of the phase at T/2 for comparison.  For continuous signals, when using this fftShift, interpolation to improve phase estimation referenced to the center point is possible (Sinc kernel better than linear or polynomial).

This has to be done anytime the input signal isn't absolutely locked to be perfectly integer periodic within the fft length in order to get reasonable phase measurement results.  e.g. for any sweeping, variable, or unknown frequencies.

Otherwise the phase will toggle during the frequency sweep.  This is because non-integer periodic sinusoids have a circular discontinuity at the window boundaries, and this discontinuity changes sign during frequency sweeps, causing the FFT phase result to have a discontinuity.  An fftshift moves this discontinuity away from the beginning or phase measurement point of an FFT, which is the reference point where both the cosine and sine, or real and imaginary components of the FFT basis vectors start with a phase of zero.


Would the ratio of Autocorrelation value and autocorrelation of noise for different data types be different? in Matlab

I am not really sure what you want to do so please comment and I will modify my answer.

The case that $x$ and $u$ seems statistically identical is that $x$ is standard Gaussian random variable and $u$ has $\pm 1$ values.

Use your notations, $x$ is Gaussian random variable with zero mean and unit variance (deducted from the fact that you used the standard rand function, even though the randn() function must be used for Gaussian pdf).

Let's stay with real $x$ for the sake of simplicity :

$$p(x) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$

with $\mathbb{E}[x] = \mu = 0 \textrm{ and } \mathrm{var}(x) = \sigma^2 = \mathbb{E}[x^2] = 1$.

Thus the probabilites $\mathrm{Pr}(x >= 0) = \mathrm{Pr}(x < 0) = 0.5$.

The processed data $u = 1 \textrm{ if } x \ge 0 \textrm{ and } u=-1 \textrm{ otherwise}$. Thus $\mathbb{E}[u] = (1) \times \mathrm{Pr}(x >= 0) + (-1) \times \mathrm{Pr}(x < 0) = 0.$

Then $\mathrm{var}(u) = \mathbb{E}[(u - \mathbb{E}[u])^2] = \mathbb{E}[u^2] = (1)^2 \times \mathrm{Pr}(x >= 0) + (-1)^2 \times \mathrm{Pr}(x < 0) = 1$.

The means and variances are just coincidentally identical.

MATLAB/Octave code

octave:6> N = 65536;
octave:7> x=randn(1,N);
octave:8> u=2*(randn(1,N)>0) - 1;
octave:9> var(x)
ans =  1.0020
octave:10> var(u)
ans =  1.0000

The case of OP question $x$ is (0,1)-uniform (by rand()) and $u$ is binary:

$\mathrm{Pr}(x >= 0.5) = \mathrm{Pr}(x < 0.5) = 0.5$.

$\mathbb{E}[x] = 0.5$

$\mathrm{var}(x) = \int_0^1 (x-0.5)^2 \mathrm{d}x = \frac{2}{3}\times 0.5^3 \approx 0.8$

The processed data $u = 1 \textrm{ if } x \ge 0.5 \textrm{ and } u=0 \textrm{ otherwise}$. Thus $\mathbb{E}[u] = (1) \times \mathrm{Pr}(x >= 0.5) + 0 \times \mathrm{Pr}(x < 0.5) = 0.5$.

$\mathrm{var}(u) = \mathbb{E}[(u - 0.5)^2] = (1 - 0.5)^2 \times \mathrm{Pr}(x >= 0.5) + (0 - 0.5)^2 \times \mathrm{Pr}(x < 0.5) = 0.25$

MATLAB/Octave code

octave:13> N = 65536;
octave:14> x=rand(1,N);
octave:15> u=rand(1,N)>0.5;
octave:16> var(x)
ans =  0.083424
octave:17> var(u)
ans =  0.25000

The variances are different !!!


Why can FFT only operate on images with specific properties?

The FFT can only operate on grayscale images, yes. Each color images contains multiple channels and only if they are displayed together (depending on which color model you use or rather what channels you use), the image will become a full color image. A good example is an RGB image. It contains three grayscale channels with varying intensity. Each channel is colored and R+G+B added together. So if you want to use the FFT on a color image, you have to take it on each color separately or convert the image to grayscale.

The Radix-2 FFT can only be used with signals with $2^n$ samples, yes. The FFT takes advantages of the number of samples in the signal by reordering the signals samples and then taking the fourier transform: Even numbered samples and odd numbered samples. By cutting the size of the transform in half, you only need a quarter of multiplications as opposed to using the DFT on the full signal. But what stops us from reordering the samples again? - Nothing. So we reorder again and again and again until we end up with a "one bit DFT", do the transform and reorder all of the samples recursively. This allows us to keep the transform simple, but it only works with powers of 2.

Using a different Radix, for example 4, the same principle applies (Radix-4 further reduces the number of multiplications compared to Radix-2), but now the signal has to be $4^n$, a power of four: 4, 16, 64, 256, 1024, ...


I think you are trying to do this:

y = np.zeros(N*len(s))
for i in range(0, len(s)):
    y[i*N] = s[i]

return y

There are other ways in python that avoid explicit loops, for example:

y = [0] * (N * len(s))
y[::N] = s


Question about normalization before estimate the spectrum

Most standard time-frequency-scale transforms are consistent with mean/standard-deviation corrections. In a more pedantic ways, linear transforms can cope well with affine transformations, in the shape of:

$$x' =  \frac{x-a}{b}$$

of which the $z$-score is a special case. The potential bias can easily be inverted.

So as long as the representation remains linear, and somehow preserves energy, such a normalization does not sound so harmful to me. With CWT tools, you can compare shapes. However, you loose amplitude, stuff becomes relative


Why adaptive filter does not work in my application

I tried your code change adaptfilt.lms to adaptfilt.nlms
also decrease the step size to 0.0001
These conditions gave me better results.
nlms is better than lms as there is stability in learning filter coefficeints.The lms algorithm could change the filter coefficients drastically.


How do I calculate the frequency of the highest amplitude wave?

This shows the subset of the data you included and the absolute value of the DFT (discrete Fourier transform) of a zero-padded version of it.

R Code Below
    #Q31040

data <- read.csv('Q31040/data.csv')

par(mfrow=c(2,1))
plot(data, type='l', col="blue")

zero_padded_data <- c(data$Volts, rep(0,1000))

omega <- seq(0,2*pi - 2*pi/length(zero_padded_data),2*pi/length(zero_padded_data))
fft_data <- abs(fft(zero_padded_data))
plot(omega, fft_data, type='l', col="blue")

peak <- which.max(fft_data)
points(omega[peak], fft_data[peak], col="red", lwd = 5)
title(paste("Peak location is", omega[peak]))


Wrong estimation of derivatives with an extended Kalman filter

Finally, I found my mistake. It was a combination of unfavorable parameters and the wrong computation of the prior error covariance. Yu et al. and Welch and Bishop actually give the same equation. I just got confused by the indices used by Yu et al. So the correct listing for the EKF should be:

% Apply extended Kalman filter equations
for k = 1:T

% Apply time update
if k == 1

    % Equation (15) in [Yu2004]
    X_prior(k, :) = f(x0);
    A = dfdx(x0);
    W = dfdw(x0);

    % Equation (12) in [Yu2004]
    P_prior(:, :, k) = A * P0 * A' + W * Q * W';

else

    % Equation (15) in [Yu2004]
    X_prior(k, :) = f(X_posterior((k - 1), :)');
    A = dfdx(X_posterior((k - 1), :)');
    W = dfdw(X_posterior((k - 1), :)');

    % Equation (12) in [Yu2004]
    P_prior(:, :, k) = A * P_posterior(:, :, (k - 1)) * A' + W * Q * W';

end

% Apply measurement update
H = dhdx(X_prior(k, :)');
V = dhdv(X_prior(k, :)');

% Equation (11) in [Yu2004] using the pseudo inverse to handle
% badly scaled matrices
K = P_prior(:, :, k) * H' / (V * R * V' + H * P_prior(:, :, k) * H');

% Equation (14) in [Yu2004]
X_posterior(k, :) = X_prior(k, :) + (K * (Z(k, :)' - h(X_prior(k, :)')))';

% Equation (10) in [Yu2004]
P_posterior(:, :, k) = P_prior(:, :, k) - K * H * P_prior(:, :, k);

end

Thank you for your comments and votes!


If you want to figure out if an image is actually a stretched version of a smaller image, you can probably assess if there is any information lost at the transformation. What you cannot do is assess spatial resolution because you don't really know how much actual space the image covers in the first place.

The idea is this: When you scale-down an image, you throw away some of the information it carries. If you were now to scale-up the scaled-down version of the image, you don't exactly get the original image because that would mean having a way to fill in the information you threw away during the first step. Therefore, the spatial frequency spectrum of a scaled-down->scaled-up image is identical to that of the scaled-down image but interpolated too. That spectrum actually stops at a specific spatial frequency component and from that information you can probably also work out how small this image used to be before it got scaled up.

For example:

This image:

Has this spectrum (notice pixelation):

And this image:

Has this spectrum:

The scaling here is extreme to illustrate the point. We go from 64 pixels in the horizontal direction to 64 pixels blown up to 1024 using cubic interpolation.

Notice how the large image can accommodate many more spatial components until its periphery but that space remains "dark" (i.e. low power at high spatial frequencies).

This is the signature of a stretched image (scaled-down->scaled-up). If the image happens to be simply blurred, then you would still detect less power in the higher frequencies but you would have to come up with a way of discriminating between pictures that are blurred and pictures that genuinely only use less spatial components than what would be expected.

To figure that out, you can integrate the spectrum radially in an attempt to recover the average envelope of the frequency spectrum. In the case of a blurred image, this envelope would tend towards known low pass frequency responses, compared to the spectrum of an image that is genuinely not blurred.

This latter part is a bit more challenging to get right and would have certain limits beyond which you could say confidently that it is a blurred image rather than a small but clear image.

The image that you bring as an example seems to be a simply blurred one. Its spectrum looks like this:

Hope this helps.


Why do digital images use the concept of ‘Stride’?

"Stride" is something completely different then width and height of an image; it's a concept related to the logical storage layout (memory addressing) of the image data array on the physical memory accessed by a program, and not to the geometric properties of the image.

Stride means "the logical memory address distance between two successive pixels of the image on a given axis" (it's the value that should be added to the pointer (which points into the buffer of the image data on memory) that would then access the data of the next sample on the chosen image axis)

For example, consider a row-major image with width $w=3$ and height $h=2$:

$$\mathbb R^{2\times 3}\ni A = \begin{bmatrix}
a_{11} & a_{12} & a_{13}  \\
a_{21} & a_{22} & a_{23}  \end{bmatrix}$$

with the internal storage looking like

$$
\begin{matrix}
\text{storage loc.} & \text{pixel} \\
0 & a_{11}\\
1 & a_{12}\\
2 & a_{13}\\
3 & a_{21}\\
4 & a_{22}\\
5 & a_{23}\\
\end{matrix}
$$

Here, the stride between elements of the same row (i.e. along the horizontal axis) is simply $s_x=1$.

The stride along the vertical axis is exactly the width ($s_y=w=3$) of the image here, since the storage is row-major and compact.

Now, it's pretty common to see formats specify a step size for row lengths, because that typically has computational advantages. For example, if we took the same $A$ and said "the row length is padded to multiples of $4$", then the storage format for $A$ would look like
$$
\begin{matrix}
\text{storage loc.} & \text{pixel} \\
0 & a_{11}\\
1 & a_{12}\\
2 & a_{13}\\
3 & \text{padding (don't care)}\\
4 & a_{21}\\
5 & a_{22}\\
6 & a_{23}\\
3 & \text{padding (don't care)}
\end{matrix}
$$

Now, the stride (i.e. how far are two elements successive on the same column apart?) is $s_y=4$, but the image width is still $w=3$, hence $s_y \mathbf \ne w$. I think this illustrates nicely why stride and dimensions are very different things.

This gets even more important for higher dimensional data.

For example, imagine a video. We can think of it as sequence of images like $A$. Let's call it $B$:

$$B = \left(A_i\right)_{i=1,\ldots N}\quad\text,$$

with $A_i \in \mathbb R^{2\times 3}$ being our frames, and thus, $B$ having the first axis signifying "frame number", the second "column", the third "row index",

$$ B \in \mathbb R^{N\times 2\times 3}\,\text.$$

Now, let's again assume the padded storage format.

The row- and column-strides stay the same.

But, what is the storage distance between two pixels with the same position, but in successive frames?

$$\begin{align}
s_t &= h s_y\\
&= 2 \cdot 4\\
&=8
\end{align}
$$

Notice that we had to use height times y-stride here because we've defined our pictures to be row-major. In column-major format, we'd use width times x-stride. (and in time-major format: I'll leave this up for you as exercise in basic understanding of storage formats).


Separating a signal into its constituent components

You do a QR decomposition of $X^T$. Then the system reads as $$Y=AR^TQ^T.$$ Now multiply with $Q$ from the right, $$YQ=AR^T,$$ which already should give you hints on the reliability of the eventual solution, since ideally, the block of $YQ$ right of the leading $m\times n$ block should be zero. If that is true, then $$A=[YQ]_{[1:m,1:n]}[R^T]_{[1:n,1:n]}^{-1}$$.

Note that in the QR decomposition using Householder reflections, you can apply the reflections computed for $X$ simultaneously to $X$ and $Y$, if $$S=I-2 vv^T,$$ then $$X_{next}=XS=X-2(Xv)v^T\text{ and }Y_{next}=YS=Y-2(Yv)v^T$$ and $$X_{final}=R^T,\; Y_{final}=YQ$$.


Frequency modulation of audio signals. Does increasing the pulse width of frequency tone result in larger operating distance?

It is very hard to give a single, correct answer because there are so many variables involved. Some things to keep in mind:

The main variable that determines the number of errors in the receiver is the ratio of energy per digit and the noise (this is the signal to noise ratio, or SNR). When you decrease the pulse duration, you reduce its energy and in consequence reduce the SNR. Increasing the pulse duration will increase the SNR and will allow you to increase the distance, all else equal.
Received power decreases with the square of the distance; for relatively large distances, this will dominate the SNR.
The technique you are using is called frequency shift keying (FSK). It is easy to implement but requires large bandwidth. Your audio setup has limited bandwidth (20 kHz in the most optimistic case). Reducing the pulse duration increases the bandwidth; at some point, your speaker/microphone combination will not have the bandwidth to properly reproduce and detect your signal.
Bandwidth considerations also suggest that you should shift your frequencies to around 10 kHz: at 17 kHz, you will soon hit the system's limit of 20 kHz. Operating at 10 kHz, you have 10 kHz of bandwidth to either side of your carrier.


Why do I hear a "blip" at the end of an audio clip of a sinusoid?

There's still the first derivative discontinuity. That gives a blip with a -12 dB/octave decaying spectrum. It's less audible (because it has less high frequencies) than the -6 dB/octave decaying spectrum of a proper discontinuity, but still audible, especially with a subwoofer...

Do a (long enough a) fade out like Fat32 says and there will be no audible click.


Similarity between two Signals: cycle by cycle

so i'll throw the cross-correlation and/or matched-filter thing at you.  a little bit similar to Karthik but maybe a little more rigor.

so you have a reference signal, $y(t)$, and a measured signal $x(t)$.  suppose they have something in common and some differences and may one of the differences is a misalignment in time.  but if you investigate what these two signals have in common (or what they don't) for a variety of different offsets in time, that dissimilarity due to time misalignment can be made non-salient (or maybe you want to know the time offset).

so first consider measuring the dissimilarity (at a variety of different offsets or "lags") as:

$$ Q_{xy}(\tau) \triangleq \int\limits_{-\infty}^{+\infty} \left( x(t) - y(t+\tau) \right)^2 dt $$

so you're looking at how much $x(t)$ and $y(t+\tau)$ are different, for different lags, $\tau$, squaring that difference or "error" so that it's always positive (this is analytically easier than if you use the absolute value), adding up the total squared-difference over all of time and you have an overall measure in how dissimilar $x(t)$ and $y(t+\tau)$ are.  then picking the lag $\tau$ so that $Q_{xy}(\tau) \ge 0$ is minimum might be considered to be a measure on how dissimilar $x(t)$ and $y(t+\tau)$ are.

the smaller this measure (which is always non-negative), the more similar $x(t)$ and $y(t+\tau)$ are.  this is the Average Squared Difference Function (ASDF) and has a similar motivation in its definition as does the more well-known and historical AMDF function.

now the cool thing you get that you don't if it's AMDF is

$$ \begin{align}
Q_{xy}(\tau) & = \int\limits_{-\infty}^{+\infty} \left( x(t) - y(t+\tau) \right)^2 \ dt  \\
& = \int\limits_{-\infty}^{+\infty} x^2(t) + y^2(t+\tau) - 2x(t)y(t+\tau) \ dt \\
& = \ \int\limits_{-\infty}^{+\infty} x^2(t) dt \ + \ \int\limits_{-\infty}^{+\infty}y^2(t+\tau) dt \ - \ 2\int\limits_{-\infty}^{+\infty}x(t)y(t+\tau) dt \\
& = \ \int\limits_{-\infty}^{+\infty} x^2(t) dt \ + \ \int\limits_{-\infty}^{+\infty}y^2(t) dt \ - \ 2\int\limits_{-\infty}^{+\infty}x(t)y(t+\tau) dt \\
& = \ E_x \ + \ E_y \ - \ 2\int\limits_{-\infty}^{+\infty}x(t)y(t+\tau) dt \\
\end{align} $$

$E_x$ and $E_y$ are the energies of $x(t)$ and $y(t)$ and not functions of $\tau$.  so $Q_{xy}(\tau)$ will be minimum with a $\tau$ that makes the integral on the right maximum.  we'll give that integral a name:

$$ R_{xy}(\tau) \triangleq \int\limits_{-\infty}^{+\infty}x(t)y(t+\tau) dt $$

and we see that

$$  Q_{xy}(\tau) = E_x + E_y - 2  R_{xy}(\tau) $$

so so $Q_{xy}(\tau)$ is minimum when $R_{xy}(\tau)$ is maximum.  in fact, with a scaling constant of $2$ tossed in there $Q_{xy}(\tau)$ and $R_{xy}(\tau)$ are upside-down mirror images of each other.

we call $R_{xy}(\tau)$, the cross-correlation function of $x(t)$ and $y(t)$ and is the inner product (if we're pretending we're living in a Hilbert space) of $x(t)$ and $y(t+\tau)$.  if we're dealing with complex $x(t)$ and $y(t)$ (we're not, i hope) then $y(t)$ should be complex-conjugated.  where $R_{xy}(\tau)$ is maximum is the time offset that $x(t)$ and $y(t)$ are most similar to each other.

the above can be rearranged:

$$ R_{xy}(\tau) = \frac{1}{2} \left( E_x + E_y - Q_{xy}(\tau) \right) $$

or

$$ \frac{R_{xy}(\tau)}{\frac{1}{2}(E_x + E_y)} =  1 - \frac{Q_{xy}(\tau)}{E_x + E_y} $$

so the cross-correlation is getting normalized by the average energy of $x(t)$ and $y(t)$.  we're always subtracting a non-negative number from 1 so that normalized cross-correlation is never larger than 1

and pick the $\tau_0 \ne 0$ where $R_{xy}(\tau_0)>R_{xy}(\tau)$ for all other $\tau$ and calculate that normalized cross-correlation and you have a measure of similarity with

$$ -1 \le \frac{R_{xy}(\tau_0)}{\frac{1}{2}(E_x + E_y)} \le 1 $$

if it's -1, they are identical but with opposite polarity.  if it's 1, they're totally similar (at that offset of $\tau_0$).


What is a correct way to downsample an impulse train

The usual way to do this is with a low-pass decimating FIR filter.  If you are going to decimate by a factor of 50 then you will probably want to break it up into multiple LPF decimation filters.  I would probably do it as /10 followed by /5.


How to perform zero-padding to 3 dimensional array

A lovely cat ("concatenate arrays along specific dimension") can do the job:

a = randn(10000,1,41);b= cat(1,a,zeros(size(a)));whos b


Difference Between Two Forms of Equations of Auto Regressive (AR) Model

Auto Regressive Model means the current output is a linear combination of previous outputs and driving noise:

$$ y \left[ n \right] = \sum_{k = 1}^{p} {a}_{k} y \left[ n - k \right] + v \left[ k \right] $$

As you can see, the current value $ y \left[ n \right] $ depends on $ m $ values before it.
The parameter $ p $ is the order of the model while $ v \left[ k \right] $ is the driving IID White Noise.

Usually if we want to estimate the parameters of the model we have 2 main cases:

The Model Order $ p $ Is Known
In this case we need to estimate $ p $ parameters $ \left\{ {a}_{k} \right\}_{k = 1}^{p} $ from a set of $ m \geq p $ samples $ \left\{ y \left[ n \right] \right\}_{n = 1}^{m} $. the classic way to do so is using the Yule Walker Normal Equations.
The Model Order $ p $ Is Unknown
If $ p $ is unknown the intuitive a approach is estimate it and then go back to (1). One could estimate is by running on a grid of options (Try model order 3, 4, 5, 6, ... and chose the one which works best in MSE) yet this might do overfit.  A classic approach to evaluate the grid is using Akaike iInformation Criterion (AIC).

In your equations:

The connection is between the output to a different signal which means it is not an AR model.

If you tell us where did you get those form we might be able to clear things up.

Welcome to DSP!


so the sampling theorem says basically:

let   $ \ x[n] \triangleq x(nT) \ $  where $ \ f_s = 1/T \ $   is the sampling frequency.

then $ x(t) $ can be reconstructed from the samples $x[n]$  using

$$x(t) \ = \ \sum_{n=-\infty}^{\infty} x[n] \cdot \mathrm{sinc} \left(\frac{t-nT}{T} \right) \ $$

which is what you get when you pass this ideally sampled signal

$$  T \sum_{n=-\infty}^{\infty} x[n] \cdot \delta(t-nT) \ $$

through an ideal brickwall low-pass filter with cutoff at $f_s/2$ and a passband gain of 1.

but what comes out of a DAC is not that, but is this:

$$  \sum_{n=-\infty}^{\infty} x[n] \cdot \mathrm{rect} \left(\frac{t-T/2 -nT}{T} \right) \ $$

so it's like the DAC has a hypothetical filter with impulse response of

$$ h(t) \ = \ \frac{1}{T} \mathrm{rect} \left(\frac{t}{T}-\frac{1}{2} \right)
 = \begin{cases}
\frac{1}{T} & \mbox{if } 0 \le t < T  \\
0           & \mbox{otherwise}
\end{cases} \ $$

so the "sample-and-hold effect" really should be called the "zero-order hold effect" and what the effect is, is to pass your signal through a filter with the above impulse response and this transfer function:

$$ H(s) \ = \ \frac{1 - e^{-sT}}{sT} \ $$

and this frequency response:

$$ H(j 2 \pi f) \ = \ e^{-j \pi f T} \mathrm{sinc}(f T) \ $$

now, these equations came out of the wikipedia references, which is why i suggested to look there first.


How to compute the impulse response from a transfer function

Your filter is the all-poles IIR, this simplifies things a bit. Normally you can write transfer function in following form:

$H(z)=\dfrac{\sum_{i=0}^{P}b_{i}z^{-i}}{\sum_{j=0}^{Q}a_{j}z^{-j}} $

Going back to the discrete time domain you will get:

$y[n] =\dfrac{1}{a_0}\left( \sum_{i=0}^{P}b_ix[n-i]-\sum_{j=1}^{Q}a_jy[n-j]\right)$

Therefore in your case it is:

$y[n]=x[n]-0.1y[n-30]$

Calculation of particular values can be done by feeding Kronecker delta to your system ([1 0 0 0 0 0 0 ...] signal). Basically you will observe that something is happening every 30 samples ;) You can easily generate the plot in MATLAB/Octave, just remember that your filter has following coefficients:

b=[1]
a=[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1]

By calling impz(b,a) function you will get similar plot (first 100 samples):


For this specific waveform as described, the following would demodulate the signal from the Frequency Shift Modulated Input into a square wave output:

This works given the duration of the "0" symbol is 12 samples by using a "delay and multiply" frequency discriminator. The multiply will have a strong double frequency component that needs to be filtered out (and can be used for timing recovery if the timing was not synchronous to the symbol rate as further detailed below), and this is done with a simple delay of 3 samples and add, which would provide a null at twice the frequency of the zero symbol.  The "1" symbol as shown by the OP is a low frequency symbol in the FSK modulation, and observe how in this case it is approximately formed using the same half cycle of the higher frequency "0" symbol followed by a zero magnitude for a full one and a half cycle duration before continuing with the negative half cycle followed by a zero magnitude for a duration for one and a half cycles (thus it is 4 times the duration of the "1" symbol.) This is convenient as the output for the case of demodulating the "0" symbol would be naturally nulled given the block diagram above without the need for subsequent filtering.

Notice with this approach that at the start of each symbol, if the symbol is a 1 the initial first half cycle result will be a sin^2 always regardless of the next symbol, but if the symbol is a zero the result will be zero throughout the duration of the zero. So in this case with this particular waveform, the way the lower frequency is constructed actually helps for a cleaner demodulation, since it provides the nulled response if the delay is the duration of the zero symbol (Maybe the waveform designer had this in mind which is why it has the peculiar form for the lower frequency).

The "Square Wave" output would be further filtered through an "integrate and dump" process prior to making a final symbol decision: sum over 12 samples in the symbol duration and decide if the result if greater or less than a decision threshold. If greater a "0" was demodulated, and if less a "1" was demodulated. Given the atypical longer duration between the symbols, four "1"'s in a row is the presence of an actual "1" symbol.

This structure can then be readily used with standard timing recovery for recovering the symbol clock boundaries. If the symbol clock is synchronized (coherent) to the waveform sample clock, the timing recovery is quite trivial. Two options show below:

The strong double frequency component can be used for symbol timing recovery for the case of the symbols no being coherent to the sample clock. If they are coherent then symbol timing recovery is simply aligning the sample edge to a start of a symbol which can easily done with a simple threshold detection at the Demodulated FSK output as shown in the block diagram (for lower SNR conditions it would be more robust to average the estimated timing position even such a threshold approach was used, such that the timing position does not change abruptly based on the result of any one threshold detection). Once detected, a modulo 12 counter starts, and every time the counter rolls over the symbol is selected at the output. The more robust timing approach would compare the roll-over timing clock and the FSK output threshold detector impulses to establish a timing error (another multiplier would do this), which would then be accumulated in a timing loop to adjust the actual start time of the cyclical timing recovery counter (again only need to do that for a very robust solution in the presence of low SNR, for the waveform as shown all that would not be necessary).

See this link where I detail non-coherent FSK demodulation approaches that further how a delay and multiply forms a frequency discriminator:

https://electronics.stackexchange.com/posts/293723/edit

With this graphic in mind; a delay multiply is a frequency discriminator:

Optional coherent techniques would have a 3 dB SNR advantage but require more processing. This could be done by matched filter correlation to the two symbols (have replicas of the two symbols and multiply and accumulate to each one to decide which symbol is observed). The delay and multiply followed by delay and add is notably simpler.


Filter recognition from transfer function

From one of your comments it appears that you've actually already answered the question. Here are a few remarks and questions that should help you gain some more understanding:

If you have a transfer function $H(z)$ (and if you assume that the corresponding system is stable), then the frequency response is obtained by choosing $z=e^{j\omega}$.
For $H(z)=z^{-k}$ figure out the magnitude and the phase of $H(e^{j\omega})$.
If the magnitude of $H(e^{j\omega})$ does not depend on frequency, what kind of filter have you got?
What is the input-output relation of the given system?


which features fundamental frequency is correlated?

Here is the answer as responded here in github.


Some subtleties in the PSD of digitally generated noise with $S_X(f) = A(f)$

You are only looking at the positive frequencies. The negative frequencies are also there. So yes, the spectrum is indeed (with some abuse of notation)

\begin{equation}
A(f) \approx \frac{1}{1+\left(\frac{f-\mathrm{sgn}(f)f_0}{b}\right)^2}\end{equation}

such that the PSD is mirrored with respect to $y$-axis.

I don't speak Mathematica, but I can tell you how it works. Your FrequencySamplingFilterKernel is in fact a simple frequency sampling approach to build FIR filters with arbitrary-magnitude frequency response (similar to Matlab's fdesign.arbmag). Given you have $N$ points of the frequency response, it essentially calculates an impulse response (by taking the inverse  DFT) whose Fourier transform is equal to your desired spectrum when evaluated at the corresponding frequencies. That is:

$$h[n]=\frac{1}{N}\sum_{k=0}^{N-1}H(k)e^{j\frac{2\pi}{N}k}$$

where $H(k)$ is the desired DFT (That is, $N$ samples from $\sqrt{A(f)}$).

Notice that the default impulse response $h[n]$ in the documentation of FrequencySamplingFilterKernel (under details and options) is an odd-length FIR filter with even symmetry. From this point everything is just regular FIR filtering with $h[n]$, and when a real process is filtered it will generate a PSD with even symmetry:
$$S_Y(f)=|H(f)|^2S_X(f)$$
where $|H(f)|^2$ is equal to $|A(f)|$ at the given $N$ desired points.

If you want to do it in Matlab:

N = 100;
f = 0:0.01:1;
b = 1;
f0 = 0.75;
A =sqrt(1./(1+((f-f0)./b).^2));
d = fdesign.arbmag('N,F,A',N,f,A);
hd = design(d,'freqsamp');
w = -pi:.001:pi;
freqz(hd,w);
x=randn(1000000,1);
y=filter(hd,x);
psd(y);

with low-order FIR approximation (N=8):


It's possible to calculate average power (proportional to RMS amplitude) with the rfft, but it's more complicated because the DC and Nyquist bins are not doubled in the full fft (and the Nyquist bin doesn't exist if n is odd).

NumPy's rfft is preferred over SciPy's, because it outputs complex numbers, which are compatible with N-dimensional arrays, rather than SciPy's "packed" format.

def rms_rfft(spectrum, n=None):
    """
    Use Parseval's theorem to find the RMS value of an even-length signal
    from its rfft, without wasting time doing an inverse real FFT.

    spectrum is produced as spectrum = numpy.fft.rfft(signal)

    For a signal x with an even number of samples, these should produce the
    same result, to within numerical accuracy:

    rms_flat(x) ~= rms_rfft(rfft(x))

    If len(x) is odd, n must be included, or the result will only be
    approximate, due to the ambiguity of rfft for odd lengths.
    """
    if n is None:
        n = (len(spectrum) - 1) * 2
    sq = real(spectrum * conj(spectrum))
    if n % 2:  # odd
        mean = (sq[0] + 2*sum(sq[1:])           )/n
    else:
        mean = (sq[0] + 2*sum(sq[1:-1]) + sq[-1])/n
    root = sqrt(mean)
    return root/sqrt(n)

More details here: https://gist.github.com/endolith/1257010


how to get the phasor(frequency domain) representation of a digitalized signal?

This is simply about representing a sinusoidal signal $x(t)$ in the following way:

$$x(t) \ = \ \cos(\omega_0t+\phi) \ = \ a \cos(\omega_0 t) \ - \ b \sin(\omega_0 t)$$

The 'method' described in your question is about identifying the constants $a$ and $b$. The phasor of $x(t)$ is then $X=a+jb$, because

$$x(t)=\text{Re}\{Xe^{j\omega_0t}\}$$

Using basic trigonometric identities you get

$$a=\cos\phi\quad\text{and}\quad b=\sin\phi$$

But you can also compute these values by correlation (as shown in your question), which is equivalent to computing the Fourier coefficients of $x(t)$ (apart from a sign difference for the coefficient $b$). Of course, the only non-zero coefficients are the coefficients corresponding to the fundamental:

$$a=\frac{2}{T}\int_0^Tx(t)\cos(\omega_0 t) \ dt,\quad  T=\frac{2\pi}{\omega_0}\\
b=-\frac{2}{T}\int_0^Tx(t)\sin(\omega_0 t) \ dt$$

These last two equations correspond to the first two steps described in your question.


Why is Frequency Shift Keying more immune to noise than ASK?

Looking through some literature shows that a side by side comparison of FSK with ASK does not give FSK a performance advantage over ASK in the presence of noise (1); however you can do some nonlinear processing with FSK in the demodulator that does produce a signal that is more tolerant to noise than ASK.

In some instances of FSK modulation/demodulation schemes, you can force the receive signal to clip (i.e. put the signal through a Schmidt trigger) so the amplitude noise is essentially eliminated.  You can't do this with the ASK signal because this eliminates the information in the modulated signal.  For the FSK signal, the clipped signal still carries frequency information with a slight phase error due to the noise signal.  There is a threshold at which noise can cause enough phase error to affect the demodulation even after clipping the signal, but below that threshold, the signal is largely immune to noise.


How to calculate SNR with sample noise and noisy signal only

The key issue here is whether the signal and noise are uncorrelated.  Assuming two real random variables x and y are both zero mean, the power of the combined signal is

$E\{(x+y)^2\} = E\{x^2+2xy+y^2\} = E\{x^2\}+E\{y^2\}+2E\{xy\}$

If the two variables are uncorrelated, then $E\{(x+y)^2\}=E\{x^2\}+E\{y^2\}$  In other words, the power of the combined signal and noise is the sum of the signal power and the the noise power if the signal and noise are uncorrelated and zero-mean.  If this assumption is valid - as if often the case in practice - then the approach you outlined would work.  It would also mean that $P_{signal+noise}=A_{signal}^2+A_{noise}^2$.  If the signal and noise are correlated, then the approach would not work.

One thing you might consider is whether you are trying to measure the signal-to-noise ratio or the signal-to-noise-and-distortion ratio.  When the audio is present, the recorded signal will suffer distortion that will degrade the overall signal quality beyond what would be predicted simply via the noise power.  This would cause the SNR estimate obtained via the above approach to be optimistic in predicting the overall signal quality.


(Radiation) Pulse detection and height characterization when pulses are stacked

Can you use deconvolution to convert these decaying impulses back into impulses?
Proof of concept:
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

impulses = np.zeros(150)

times =   [3,  20, 30, 40, 45, 50, 55, 80, 90]
heights = [-8, -7, -1, -9, -1, -2, -1, -8, -1]
impulses[times] = heights

b, a = signal.butter(1, 0.04)

filtered = signal.lfilter(b, a, impulses)
fig, ax = plt.subplots(3, 1, sharex=True)
ax[0].plot(filtered[:151])
ax[0].set_xlim(0, 100)
ax[0].set_title('Original')

impulse = -10*signal.unit_impulse(50)
template = signal.lfilter(b, a, impulse)
ax[1].plot(template)
ax[1].set_title('Template')

deconvolved, remainder = signal.deconvolve(filtered, template)
ax[2].plot(deconvolved)
ax[2].set_title('Deconvolved')

Then they are separated from each other and you can count their heights and events more easily.
Depends what the properties of the original bumps are, but they look like lowpass filtered impulses, which would mean they are all from convolution with the same impulse response, so this should work?


How can I increase image size by zero padding?

I suppose the array names are not descriptive of the contents and you have simply repurposed them. To zero-pad from 8x8 to 32x32:

for (int i = 0; i < 32; i++)
    for (int j = 0; j < 32; j++)
        pspectrum[i][j] = 0;
for (int i = 0; i < 8; i++)
    for (int j = 0; j < 8; j++)
        pspectrum[i][j] = powerspectrum[i][j];

The purpose of this is that the discrete Fourier transform (DFT) of a larger image has a higher frequency resolution, so you can resolve for example the dominant frequency more accurately. It is the same as really well done interpolation in the frequency domain.


Viterbi decoder on a bit stream where the preamble is also coded

Is the 500-bit encoded message part of a continuous stream being decoded by the Viterbi decoder, or is it a 500-bit packet that is decoded individually?   If part of a stream, then decoding continuously and searching for the preamble indicates the beginning of the message.   If packetized, then the decoder memory needs to be managed so that the first decoded bit location is known.   This will be dependent on the decoder implementation.


For your non-real-time application, you do not need OS or OA.  Direct form convolution (such as Octave function conv()) is easy and gives the same result.
But since the title of your post is "Filter size vs. FFT size and Overlap add", this graphic directly addresses that issue.


How to Classify Image Filter as Low Pass Filter (LPF) or High Pass Filter (HPF)?

Welcome to DSP Community.

General

We assume 2 types of filters: LPF or HPF.

Classifying Filter Type

Usually, if it is a well planned LPF and well Planned HPF a simple test will do.
Calculate the sum of all coefficients.
The sum of the coefficients is the first element of the DFT of the signal.
It means it is the DC gain and well behaved LPF has gain of 1 and well behaved HPF have DC gain of 0.

So sum the coefficients and:

If the sum is $ 1 $ then this is an LPF.
If the sum is $ 0 $ then this is an HPF.

Pay attention that those method can't differentiate between HPF and Band Pass or between LPF and Band Stop.

The more general method is to look on its DFT Spectrum to see its behavior.

Converting Filter Type

There is very intuitive formula - ALL_PASS_FILTER = LPF_FILTER + HPF_FILTER.
It's not accurate but this is the idea.

Be more accurate in the Frequency Domain All Pass Filter is a filter which has value of 1 for any frequency.
Assuming we have LPF Filter then HPF_FILTER = ALL_PASS_FILTER - LPF_FILTER.
In the spatial domain it means $ h \left[ n \right] = \delta \left[ n \right] - l \left[ n \right] $.

Where $ \delta \left[ n \right] $ is the filter with 1 in its origin and $ 0 $ anywhere else.

In your case $ \delta \left[ n \right] = [0, 0, 0, 0, 1, 0, 0, 0, 0] $.
So the HPF Filter becomes $ h \left[ n \right] = \delta \left[ n \right] - l \left[ n \right] = [-1/9, -1/9, -1/9, -1/9, 8/9, -1/9, -1/9, -1/9, -1/9] $.


How can I detect pipeline cracks using OpenCV and Python?

Try this:

Don't convert the image to grayscale, instead convert to HSV
Then, by changing the range of H,S and V, try to mask the crack (using a trackbar can be very helpful)

For example:

lower_value_blue = np.array([80,106,5])    # creating a range of HSV colours
higher_value_blue = np.array([150,255,150])
mask1 = cv2.inRange(blur1,lower_value_blue,higher_value_blue)

This one mask blue colour, without using trackbar

Hope you will get the result


It's way easier to use the complex representation of cosine not involving the real part operator. In your case this translates into expressing

$$ \cos(h(t))=\frac{e^{jh(t)} + e^{-jh(t)}}{2} $$

Note that this yields an autocorrealtion of
$$ R_s(\tau)=\int_{0}^{T_s} \frac{e^{jh(t)} + e^{-jh(t)}}{2}\frac{e^{jh(t-\tau)} + e^{-jh(t-\tau)}}{2} dt $$


Why is the bandwidth for this QPSK system larger than expected?

In gnuradio, the throttle block only determines the speed at which your computer is processing samples - it does not determine what the actual sample or baud rates are. As an example - consider that you generate random PSK symbols at 1 sample per symbol, and then upsample by 2 with a root raised cosine w/ 35% exccess BW so we're now at 2 samples per symbol (this is what the constellation modulator block does in your flowgraph). What is the absolute baudrate that you're running at? Answer: Whatever you want it to be. If I hook that I/Q stream up to an RF front end, I could theoretically transmit at any baud rate using those exact same samples. E.g. if the DAC is running at 1MSPS (assuming this is an I/Q DAC for simplicity like you might consider the ad9361) then my baud rate would be 1MSPS/(2 samples/symbol) = 500k symbols/sec. Assuming my PC could keep up (and I increased my throttle block so that my PC isn't slowing me down from producing samples), I could instead run the DAC at 10MSPS to transmit at 5M symbols/sec.
So - your question will probably be - how do I simulate a certain baud rate without a DAC/front-end? Well - it's whatever you choose it to be. In your flowgraph, you set the FFT sink to think that the sample rate is 500 kHz. The data you have hooked up to that sink is coming in at 2 samples/symbol, but factoring in the excess bw, the relative bandwidth will be about (1.35/2.0)*500e3 = 337.5 kHz which is exactly what the spectrum is showing you. If you did further upsampling and then looked at the FFT, the signal would appear narrower and vice versa.


White noise is not "WSS by nature" whatever you mean by that phrase but it can be treated as a (zero-mean) WSS process insofar as its effects in linear systems are concerned.
For example, standard linear system theory ways when the input to an LTI system is an ordinary WSS process $\{X(t)\}$ with autocorrelation function $R_X(\tau)$, then the output of the LTI system is a WSS process
$\{Y(t)\}$ with autocorrelation function $R_Y(\tau)$ given by
$$R_Y = h\star \tilde{h} \star R_X \tag{1}$$
where $h(t)$ is the impulse response of the LTI system and $\tilde{h}(t) = h(-t)$ is the time-reversed impulse response of the LTI system. The power spectral densities are related as
$$S_Y(f) = |H(f)|^2S_X(f)\tag{2}$$
where $H(f)$ is the transfer function of the LTI system. If $\{X(t)\}$ is a white noise process with autocorrelation function $K\delta(\tau)$ and we pretend that $(1)$ and $(2)$ are still applicable, we get that $\{Y(t)\}$ is a zero-mean WSS process with autocorrelation function $R_Y = K\cdot h\star \tilde{h}$ and power spectral density $S_Y(f) = K\cdot |H(f)|^2$. Of course, mathematicians would laugh at this calculation but physical experiments using the naturally occurring thermal noise in electrical circuits as a stand-in for a white noise process show that these results are pretty close to reality. As engineers, we seek equations that match the universe as we observe it (physicists seek universes that match their equations while mathematicians don't care) and so we go blithely on our way treating white noise as a WSS process in linear systems and everything works out OK. The troubles start when we start treating white noise as a WSS process in nonlinear systems and the world comes crashing down about our ears and we need to start paying attention to what the math people are saying.
White noise is referred as white Gaussian noise if we pretend or claim or insist that $\{Y(t)\}$ is a Gaussian process which means that not only are all the random variables $Y(t)$ Gaussian random variables, but every finite set $\{Y(t_1), Y(t_2), \cdots, Y(t_n)\}, n \geq 2,$ of random variables is a set of jointly Gaussian random variables. Standard random process theory says that when a Gaussian process is passed through an LTI system, the output is a Gaussian process but this fact does not allow us to reverse-engineer the result and claim that all the $X(t)$'s are also Gaussian random variables.


Weighted Sum of Auto Correlation - Lower Bound

Defining $ a \left[ k \right] = {2}^{- \left| k \right|} $.
Moreover, the Auto Correlation function of $ v $ defined as $ {r}_{vv} \left[ k \right] = \left \langle {v}^{\left( 0 \right)}, {v}^{\left( k \right)} \right \rangle = \sum_{n = -\infty}^{\infty} {v}_{n} {v}_{n - k} $.
Pay attention that Auto Correlation is Hermitian Function.

Using the definition of Convolution one could write:

$$ \left( {r}_{vv} \ast a \right) \left[ 0 \right] = \sum_{k = -\infty}^{\infty} {2}^{- \left| k \right| } \left \langle {v}^{\left( 0 \right)}, {v}^{\left( k \right)} \right \rangle $$

Using the Convolution Theorem one could write that:

$$ \left( {r}_{vv} \ast a \right) \left[ 0 \right] = \int_{- \pi}^{\pi} {R}_{vv} \left( \omega \right) A \left( \omega \right) d \omega $$

Where $ R \left( \omega \right) $ and $ {R}_{vv} \left( \omega \right) $ are the DTFT of $ {r}_{vv} \left[ k \right] $ and $ a \left[ k \right] $ respectively.

One should notice the DTFT of $ a \left[ k \right] $ is defined only one sided. Yet since its symmetrical it can well calculated:

$$
\begin{align*}
A \left( \omega \right) & = DTFT \left\{ a \left[ k \right]  \right\} = \sum_{k = -\infty}^{\infty} a \left[ k \right] {e}^{-j \omega k} = \sum_{k = 0}^{\infty} {2}^{-k} {e}^{-j \omega k} + \sum_{k = 0}^{\infty} {2}^{-k} {e}^{j \omega k} - 1 \\
& = \frac{1}{1 - 0.5 {e}^{-j \omega}} + \frac{1}{1 - 0.5 {e}^{j \omega}} - 1 = \frac{1 - {c}^{2}}{1 - 2 c \cos \left( \omega \right) + {c}^{2}} = \alpha > 0 \quad \forall c < 1
\end{align*}
$$

In the above $ c = {2}^{-1} = 0.5 $ yet actually this will hold for any $ c < 1 $.

So the integral is given by:

$$
\begin{align*}
\int_{- \pi}^{\pi} {R}_{vv} \left( \omega \right) A \left( \omega \right) d \omega & = \int_{- \pi}^{\pi} {R}_{vv} \left( \omega \right) \frac{1 - {c}^{2}}{1 - 2 \alpha \cos \left( \omega \right) + {c}^{2}} d \omega \\
& \geq \alpha \int_{- \pi}^{\pi} {R}_{vv} \left( \omega \right) d \omega = \alpha {\left\| v \right\|}^{2}
\end{align*}
$$

As requested.

By the way the result must be real since $ a \left[ k \right] $ is symmetric and $ {r}_{vv} $ is hermitian function and hence its transform is real.


Hint: I think you want to determine the [Rx,Ry,Rz] required to make the left side of your equation look like [1 0 0], [0 1 0], and [0 0 1]. This will get you A,B,C expressed in terms of X,Y,Z.


how sobel edge detector as a first order derivative is turned to a 2 dim filter?

Background

The Sobel edge detector was introduced back in 1968 by Irwin Sobel and Gary Feldman as the Sobel-Feldman operator. In broad strokes, 'edges' in images are related to gradients, which motivated their development of a discrete differentiation operator. The Sobel-Feldman operator only computes an approximation of the gradient and not the actual gradient. However, this integer valued operator is simple and fast, which makes it useful.

  Directionality

The Sobel-Feldman operator is only capable of estimating gradients along mutually perpendicular directions separately and not jointly. Since an image is typically a matrix itself, the operator provides gradients along rows and gradients along columns. This limitation is important since we will use this property to extend the operator to multiple dimensions.

  Matrix Decomposition

The Sobel-Feldman operator in 2 dimensions is given as a kernel matrix. For image matrix $\mathcal{I}$, the Gradient operator (Sobel edge detector) for x- and y- directions is given as:

$\mathbf {G} _{x}={\begin{bmatrix}+1&0&-1\\+2&0&-2\\+1&0&-1\end{bmatrix}}*\mathbf {\mathcal{I}} \quad {\mbox{and}}\quad \mathbf {G} _{y}={\begin{bmatrix}+1&+2&+1\\0&0&0\\-1&-2&-1\end{bmatrix}}*\mathbf {\mathcal{I}} $

Sobel kernels can be decomposed as the products of an averaging and a differentiation kernel. For example,

$\mathbf {G} _{x}= {\begin{bmatrix}+1&0&-1\\+2&0&-2\\+1&0&-1\end{bmatrix}}={\begin{bmatrix}1\\2\\1\end{bmatrix}}{\begin{bmatrix}+1&0&-1\end{bmatrix}}$

This is the key idea in extending the Sobel operator from 1-D vector to 2-D matrix. The operator uses a triangle-kernel for averaging (smoothing) in the direction perpendicular to the desired gradient direction and the Sobel differentiator in the desired direction.

  Higher Dimensions

Sobel–Feldman filters for image derivatives in different dimensions with $x,y,z,t\in \left\{0,-1,1\right\}$:

vector: $h_{x}'(x)=h'(x); h_{x}'(x)=h'(x)$

matrix: $h_{x}'(x,y)=h'(x)h(y) h_{x}'(x,y)=h'(x)h(y)$

Tensor: $h_{x}'(x,y,z)=h'(x)h(y)h(z) h_{x}'(x,y,z)=h'(x)h(y)h(z)$

Tensor:  $h_{x}'(x,y,z,t)=h'(x)h(y)h(z)h(t) h_{x}'(x,y,z,t)=h'(x)h(y)h(z)h(t)$

  Tensor Sobel operator

For 3-D (x,y,z), the Sobel-Feldman operator in z-direction is:

$h_{z}'(:,:,-1)={\begin{bmatrix}+1&+2&+1\\+2&+4&+2\\+1&+2&+1\end{bmatrix}}\quad h_{z}'(:,:,0)={\begin{bmatrix}0&0&0\\0&0&0\\0&0&0\end{bmatrix}}\quad h_{z}'(:,:,1)={\begin{bmatrix}-1&-2&-1\\-2&-4&-2\\-1&-2&-1\end{bmatrix}}$

  Summary

The Sobel operator is only a poor approximation of gradient using a center-difference operator.
The operator is only informative about 1 direction at a time.
The operator is extended to higher dimensions by using a triangle shaped averaging filter for other perpendicular directions (dimensions).
There isn't a deep theoretical reasoning behind this operator. It is a discretized version of a differential (gradient) and Gaussian filter (smoothing). But it is simple and fast.


I think this overly conservative. If you apply it to an IIR filter, you would conclude infinite bit growth which isn't the case.
For a filter with impulse response, $h[n]$, the "worst case" gain the filter can provide to the input is given as the absolute some the coefficients, i.e.
$$G_{max} = \sum |h[n]| $$
This means that if $x[n] < x_{max}$ then $y[n] < x_{max} \cdot G_{max} $ for all $n$.
You are guaranteed not to overflow as long as you accommodate this maximum gain.
Even that's too conservative for most applications. The input signal that will create this maximum gain is basically $sign(h[-n])$ which is extremely unlikely to occur. So in practice a statistical approach works better: Run a sizable number of representative input signals through the filter and look at the probability distribution function of the output.
At the end of the day this is a trade off between the noise from occasionally clipping an output sample vs the quantization/rounding noise you create on every output sample.


Implementing headphone transfer function (HPTF)

Properly equalization improves the spectral accuracy of a binaural playback. It does not improve the spatial properties of the HRTFs other than that vertical localization is heavily dependent on proper spectral playback.

The best headphones for binaural playback tend to be diffuse field equalized and "acoustically open". This means that the headphone itself doesn't change the acoustic radiation impedance as seen from the ear canal into the outside world. It can actually be shown that headsets that are "acoustically open" do not require specific headphone equalization.

I'm not very optimistic about using Apple EarPods for this application. There is a huge amount of placement variability just through the way how deep and how hard you shove them into the ear canal. Small leaks between the pod and the walls of the ear canal make big changes in the frequency response and these vary greatly from person to person and even from placement to placement. So even if there was a set of functions published, it would not be a good fit for a specific person and placement.


Can signal processing be applied successfully to price data from financial markets?

Well of course you can apply it and many people try. If anyone has been really successful in doing this they probably are sitting now on a nice warm beach sipping Martinis, so this information may not be easy to come by.


I don't think filtering is the way to go. The spectra of the different signals (the desired voice, the screamy children and the cars) are probably overlapping, so there is no real way to get rid of all that noise without spoiling your signal.

Actually, the approach would depend in how loud the noise is in comparison to the signal. If the speaker is near the microphone, and therefore the relative noise level is ''low'', this sounds like a VAD kind of problem.

If noise level is quite low, an energy-based VAD may work: after windowing the input signal (let's say to 20ms windows), each window's energy can be computed. Then, a threshold can be applied in order to decide if that window is noise or voice. You would simply discard the windows that aren't load enough to meet the threshold. This is unarguably the simplest approach.

A simple improvement to this kind of VAD, assuming you only have one speaker, is to reduce the bandwidth in which you apply the stated algorithm. This way, if the kids voice spectrum begins at about e.g. 500Hz, and the cars induce some low frequencies up to e.g. 80Hz, you could use a band pass (e.g. 100Hz - 300Hz, assuming the speaker fundamental frequency is about 200Hz) filtered signal to decide wether each window is or isn't noise, and select those who aren't in the original signal. I know this is not ideal, as too many coincidences are needed to have a reliable system.

Unfortunately, this won't work when noise gets louder. Furthermore, some phonems may be discarded due to their low energy, such as /s/ or /z/. This just makes speech sound weird. There are some more advanced VAD techniques such as Sohn's A Statistical Model-Based Voice Activity Detection. You can find some Matlab implementations in the net. I have used this in addition to spectral substraction noise reduction (which doesn't seem to be needed in your case) in a noisy environment (babble + WGN) and results were susprisingly good.

Matlab code for using Voicebox's Sohn VAD.

in_file='yourfile.wav';
[s,fs,wmode,fidx]=readwav(in_file,'p',-1,3500); %3500 is the number of samples to discard at the beginning of your file.
[y1,zo]=vadsohn(s,fs); %y1 is 1 when there is voice, 0 when there isn't
y=s(1:length(y1)).*y1; %apply mask
y(1)=0;

for i=2:length(y)
    if y(i-1)<y(i)
        StartTime=(i-2)/fs
    elseif y(i-1)>y(i)
        EndTime=(i-2)/fs
    end
end


Is it possible to design linear-phase filters that sum to a flat frequency response?

Yes, of course.

If so is it practical to use them in real-time audio processing for as many as 10 bands?

That depends on your definition of practical. The problem at audio is that FIR filters at low frequencies get really long. It depends on the lowest frequency, desired steepness and amount cross band rejection. but a few thousand taps is typical. These can be fairly efficiently implemented using an FFT based scheme like Overlap Add, but that introduces a lot of latency: typically it's twice the filter length. That's prohibitive for many real time applications.
If you want to build a standard 10-band audio equalizer: this can be done very efficiently with cascaded biquads where each biquad is a "peakingEQ" filter and simply recalculate the filter coefficients whenever a band changes gain.
Parallel IIR bandpass filters are difficult to phase manage but it's possible. Odd Butterworth and Linkwitz Riley filter sum to a flat magnitude response but NOT a flat phase response, so you need to put compensation allpass filters in the parallel paths


Yes that is correct-- if you sampled the audio recording at 48 KHz and then played those same samples back at 44.1 KHz you would be "time-stretching" the recording and it would be slower. However you can also resample the samples by the ratio of 147/160 to then have the 44.1KHz samples play back at the same speed.


SNR estimation in segmented speech signals

The given equation should be used with caution. Not always we can separate powers of two mixed signals just like that. It will require them to be orthogonal (or independent and at least one of them zero-mean).

In your formula the power of signal $y[n]$ is defined as $$P=\frac{1}{N}\sum_{i=1}^{N}|y[i]|^2$$
where $N$ is the length of $y$. You shouldn't forget the normalizing factor $\frac{1}{N}$.

Just to make it clear, $x$ should be a recording of signal plus noise, and $n$ should be just noise (silence/ no speech). The oscillogram of $n$ should look like the background of oscillogram of $x$ for the formula to make sense.

Regarding the negative values, logarithm of a number less than one is negative. It implies that the numerator is smaller than the denumerator, or the signal power is less than that of noise. A negative SNR can be seen in practice, where the signal is weak in nature. GPS signal can be a good example. But I really think that your noise is polluted with some sort of harmonics that don't appear in $x$, that is why the overall result is less than one and the SNR becomes negative.


Do you really have to always satisfy Nyquist with all signals if your purpose is to reconstruct a signal?

No, because this is a sufficient condition (for regularly sampled signals), and not a necessary one. This condition restricts the space of all possible continuous signals to a subspace of discrete sequences that contain the same information.

Suppose that you can constrain the signal space, eg limited band-width, positivity, parametric models, sparsity, etc. Then many theoretical works ensure perfect reconstruction (in theory) under Nyquist-like conditions. Literature abounds, under the names of compressive or compressed sensing, sparse sampling, finite-rate of innovation, etc. Some are "exactly exact", some are exact within a high probability, under quantization bounds, etc.

The topic is very wide and active, I cannot provide more details unless the question is more precise. A couple of papers however:

A Systematic Review of Compressive Sensing: Concepts, Implementations and Applications, 2018
Sparse Sampling of Signal Innovations, 2008
A Review of Sparse Recovery Algorithms, 2018


Image cleaning - removing all regular shapes ( circles, squares, text, icons etc. )

I think the best way to get rid of those artifact shapes, is to implement an Shape Detection. Then, when you had coordinates of the unwanted shape, you can cut it out of your image. Then must fill its place with interpolating (through using background intensity values at the edges of the removed patch).

Now, there are various methods for shape detection and even OpenCV OpenCV Implementation  and Aforge  Aforge Implementation  support some of them. First you must get image edges, then using Hough transforms for lines and circles, detect shapes.

http://www.emgu.com/wiki/index.php/Shape_(Triangle,_Rectangle,_Circle,_Line)_Detection_in_CSharp

Finally as I said, cut out the image and fill out the removed patch through interpolation.


Using voice recognition algorithms to reduce noise in financial time-series?

First, could you please explain what you mean by "noise" in financial time series? Are you talking about extracting a trend from time series?

In itself, speech recognition consists in converting a signal into a higher-level semantic, symbolic representation of this signal. It does not achieve any denoising, so it would be surprising to use speech recognition techniques to perform any modification of a signal.

Maybe you have something else in mind?

Anyway, a speech recognition system is made of the following layers:

A feature extraction module, converting chunks of audio signal into a vector of features capturing a representation of sound which is somewhat noise-robust and very roughly pitch invariant.
An acoustic model (usually, GMMs) computing the match (a likelihood score) between the audio feature vector and a very large set of phone models.
Optionally, an adaptation module "closing the loop" between the first two steps - learning a better parametrization that makes the extracted features conform more to the known model.
Then a hierarchy of HMM-like models (phones models, word pronunciation models, language model) - all combined into a single giant finite state transducer.

The first step doesn't make sense at all for financial data, since the parametrization used for speech recognition is deeply rooted into the way speech sounds are produced, and perceived by the human ear.

There could be financial applications of the kind of tools (Hidden Markov Models and their variants) used for the following steps... But then, it's a very well known artillery that is not specific to speech recognition and they are indeed widely used in time series analysis - for segmentation, change detection, detection of outliers, or automatic labeling of time series with the kind of "patterns" used by technical analysts.


Check out the numpy docs ( https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfft.html )

To sum up, real DFT normally takes a real input of N samples and returns a complex output of N/2+1 values.

What you end up with are the positive frequencies (up to the Nyquist frequency) and the
negative frequencies, which are just complex conjugates of their positive counterparts, are discarded or are not even computed.


Help me understand this simple code snippet to get the magnitude frequency response

It looks like this mul-floor-div is an incorrect attempt to limit some numeric print-out to 3 or 6 decimal places.  Better to do that during the print formatting.


Leakage in power spectral density estimation

Could somebody explain what does this mean? Why does a window leak
  power into adjacent frequency bins?

You can think of frequency bins as imposing an artificial frequency grid which allows only certain frequencies to fit nicely within it (producing a sharp peak) and all other frequencies are smeared across adjacent bins thus producing spectral leakage.

This means that only those frequencies that are centered exactly on bin frequency will not produce leakage.

Since the bandwidth of a bin (i.e frequency resolution) is simply the
sampling rate divided by the dft length, we can conclude that only frequencies that are exact integer multiples of the frequency resolution will not produce spectral leakage.

For example, if the sampling frequency is 8000 Hz and dft length is 1024, the frequency resolution is computed as

 frequency_resolution = 8000 / 1024 = 7,8125 Hz

So, if your signal contains frequencies which are integer multiples of 7,8125, you'll get
no leakage (e.g 1*7,8125, 2*7,8125, 3*7,8125 and so on)

If, however, your signal contains a frequency that is somewhere between 2 adjacent bins, you'll get spectral leakage, which I think is intuitive because no single bin can
fit that frequency component exactly.

UPDATE (trying to answer the comments below):

It's true that windowing does introduce some additional leakage on its own. But even if you don't use windowing at all (which corresponds to a rectangular window), you'll still get spectral leakage (provided that the frequency under examination is not an integer multiple of the frequency resolution).
In fact, since using a rectangular window (no window at all) produces the worst leakage, using overlapping windows other than the rectangular helps reduce the spectral leakage.

So, why do we have spectral leakage in the first place?

Remember that DFT assumes that a signal is periodic. Since the frequency resolution is  fixed (e.g 7,8125 Hz in the example above) for every bin, which in this case corresponds to exactly 1024 samples, we can see that a single bin can hold exactly one cycle of a 7,8125 Hz frequency , 2 cycles of 15,625 Hz frequency, 3 cycles of 24,4375Hz frequency and so on. This is true for every bin, which means that there'll be no discontinuities either within a bin or between adjacent bins.

However, if we now change our frequency to something that fits between two bins, things change. This new frequency will take either more (or fewer) than 1024 samples (one or more cycles), which means that we'll get discontinuities either within a bin or between neighbouring bins and this in turn creates spectral leakage.

Hope this helps.


Cross correlation and cross power spectrum for signal + noise

$S_{xy}=S_{xx}$ is true but this does not mean $F_x*F_y = F_x*F_x$ when $y$ is random.

The actual definition is $$S_{xy}(f)=\lim_{T\rightarrow \infty}\textbf{E}\{F^*_xF_y\}$$ for random signals. Since $x$ and $v$ are uncorrelated, this would be equal to $$\lim_{T\rightarrow \infty}\textbf{E}\{F^*_xF_x\} = S_{xx}.$$ But your next assumption is not correct.

$F_x$, $F_y$ are only based on single realization of $x$ and $y$ respectively. For random signal, $S_{xy} \ne F_x^*F_y$


Using overlap add for drawing spectrogram

Apply window function to each overlapped window, and then simply draw
  the FFT magnitude response of each window as if the spectrum
  corresponds to the center time of the window.

Your approach is partially correct. But in the second half of sentence is not how you draw a spectrogram. You shouldn't align the FFT output in the same direction as time. It should be perpendicular direction as time. After you obtain FFT Magnitude, look at it from "above". Means the colour grades that you give depends on the magnitude value. This will be line corresponding to first window. Similarly, for second window, you place it next to it (since you took it after certain interval of time). An example below shows simple spectrogram having 3-bin FFT with middle bin having higher value compared to other bins.


How does one describe the “frequency response” of a mixer?

A mixer is time-varying, but an ideal mixer is linear.  In radio usage, a mixer is periodically varying, which actually does give you an avenue to evaluate it in terms of a frequency response.
In the time domain, in a superheterodyne or direct conversion receiver, an ideal mixer takes a received signal $r_{rf}(t)$ and multiplies it by a sine wave to produce an intermediate-frequency signal $r_{if}(t) = A \cos \left (\omega_{lo} t \right)r_{rf}(t)$.  Note that the fact that this is a simple multiplication of the input by some independent quantity makes the system linear; the fact that the independent quantity isn't constant makes it time-varying.
You can take the Fourier transform of the above equation:$$R_{if}(\omega) = \mathcal{F}\left \lbrace r_{if}(t) \right \rbrace = \mathcal{F}\left \lbrace A \cos \left (\omega_{lo} t \right)r_{rf}(t) \right \rbrace$$
The Fourier transform of a multiplication is the convolution of the Fourier transforms of the multiplicands, so $$R_{if}(\omega) = \mathcal{F}\left \lbrace A \cos \left (\omega_{lo} t \right) \right \rbrace * R_{rf}(\omega)$$
The Fourier transform of a cosine wave at $\omega_{lo}$ is a pair of Dirac impulses at $-\omega_{lo}$ and $\omega_{lo}$, so the output of the mixer is, in turn, $$R_{if}(\omega) = \frac{A}{2} R_{rf}(\omega + \omega_{lo}) + \frac{A}{2} R_{rf}(\omega - \omega_{lo})$$
You can extend this for a mixer that multiplies by something other than a sine wave, such as a chirp or (commonly for actual RF circuits) a square wave.
So all of this was to get you ready for -- tada! -- the frequency response of a mixer.  The ideal mixer responds as above.  A real mixer has a frequency response that depends on what it does to the RF signal before it hits the mixing circuit, and what it does to the LO signal after it hits the mixing circuit.  So -- at this level of analysis -- you can expect the output of a more real mixer to actually look something like $$R_{if}(\omega) = H_h(\omega) R_{rf}(\omega + \omega_{lo}) + H_l(\omega) R_{rf}(\omega - \omega_{lo})$$
In the nomenclature I'm using within this answer, the $H_h(\omega)$ or  $H_l(\omega)$ would be your "frequency response of the mixer", depending on whether you were interested in the high-side or low-side mixing product.
The one thing that you can't do on a simple spectrum analyzer is just get the $S_{21}$ figure, at least assuming that the analyzer is looking at the inputs and outputs around the same frequency.  You need to input one frequency, and then look for the output at the correct shifted frequency.


How to extract the maxima and minima of a noisy period input signal

If you have the entire signal, you can then run a Lowess smoothing. This algorithm is statistically robust and provides locally acceptable smoothing. Outliers will be tolerated as well. After smoothing your signal in that way, you might apply a derivative based (or even derivative free) peak detection in order to obtain the signal peaks. These will be your local minima/maxima.

You could use MATLAB or R and port smoothing algorithms to C, they are pretty straightforward and open source:

http://www.mathworks.com/help/curvefit/smooth.html

And here are your options in C:

1) https://code.google.com/p/variationtoolkit/ for Lowess/Loess.
2) For peak detection, there are many tools available in C, one of them being:

https://github.com/xuphys/peakdetect


Why can adaptive IIR filters result in unstable solutions?

The IIR filter doesn't have to be unstable, but it has the potential of being so; unlike the FIR case which doesn't have even the potential.

One reason for the (potential) unstability of an IIR (adaptive) filter is the numerical issues due to coefficient quantization. When the poles are closer to unit circle this will be critical.
This is especially important if you are using coarse quantization (like in an legacy 8 bit system?) or you are forcing the limits of your numerical precision.

Furthermore, during the adaptation process, the chaotic behaviour of the input (which is possibly nonstationary) may lead into unbounded updates of the coefficients. When the errors get quite large during a transient stage, the resulting update on the coefficients can also be quite large leading (depending on the adaptation algorithm) the filter to fall into an unstable region...


Analytical Solution for the Convolution of Signal with a Box Filter

DO NOT PASS GO; DO NOT Collect $200; DO NOT Take Fourier transforms, or worse yet, FFTs

You do the convolution exactly the way you would do any other convolution: start with the basic convolution integral (not what you wrote) and apply the properties of the signals that you are using to come up with an easier calculation.  Begin with
$$\int_{-\infty}^\infty y(x^\prime)h(x^\prime-x)\mathrm dx^\prime\tag{1}$$
and use the result that $h(\cdot)$ is nonzero only when its argument ($x^\prime-x$ in this instance) lies in the interval $\left[-\frac{\Delta}2, +\frac{\Delta}2\right]$. So, the integrand of the convolution integral is $0$ whenever $x^\prime$ is such that
$$x^\prime-x > \frac{\Delta}2 \implies  x^\prime > x + \frac{\Delta}2$$
or that
$$x^\prime-x < -\frac{\Delta}2 \implies  x^\prime < x - \frac{\Delta}2.$$
This allows us to simplify the convolution integral $(1)$ into
$$\int_{x - \frac{\Delta}2}^{x + \frac{\Delta}2} y(x^\prime)\frac{1}{\Delta}\mathrm dx^\prime\tag{2}$$
where I have substituted the nonzero value of $h$ in the integral to save a step.  Can you take it from here?


Quantifying the energy of a voiced frame

Are you being asked to write a simple vocoder?

Your question is somewhat vague and confusing, but it seems like you are being asked to determine parameter thresholds for deciding when a speech frame is voiced (has fundamental periodic component and harmonics), unvoiced or frictional (noise generated), and silent (do not encode).

It looks like you would be using the pitch frequency estimates to determine if the frame is voiced/unvoiced, and signal energy to determine if it is a silence frame.

If he said it's experimental, I think he is asking you to look at the energy and pitch content of an audio signal and see if you can see differences that will help you choose thresholds.  For example, get some audio of someone speaking.  When the person is vocalizing (i.e. vowel sounds like "aaaahh" and nasally sounds like "mmmmmmmm") the spectrum should show periodic components.  For unvoiced sounds (i.e. short "th" "t" "s" "f" type sounds) the spectrum will look more like noise.  And when there is silence, the signal should be relatively low energy.  Take what you find and use that information to determine your thresholds.


How to estimate and compensate for doppler shift in wireless signals?

This is a very common communications problem. Look in a textbook for "frequency synchronization"; whole books have been written on these and related topics. The technique that you would choose is a function of the specifics of your system. There are two common sources of frequency offset:

Differences in frequency between the reference oscillator at the transmitter and receiver. This error is usually small, depending upon the precision of the timebases available, and can be mitigated at some cost. Cheap crystal oscillators typically achieve 50 parts per million of error or better (although this does drift as the crystal ages). If you have a larger budget, you can use something like a Rubidium standard, which provides ~1 part per trillion of frequency error. A cheaper and increasingly common approach is to use a GPS receiver that has precision frequency output (typically 10 MHz). The highly accurate timebase available from the GPS constellation can be used to train the reference precisely on-frequency.
Physical dynamics effects between the transmitter and receiver. One notable example in which this comes into play are for satellite communications applications (especially in lower orbits), where the satellite is moving (and accelerating) very rapidly relative to any observer on the Earth. The high radial velocity of the satellite toward the receiver will cause a Doppler shift, and any change in that velocity caused by its orbit will cause that shift to change over time. In applications where you have these sort of dynamics, you typically can't mitigate it much, so you're left with building a receiver that tolerates the effects.

So how does a receiver synchronize with the transmitter in these cases?

One common approach useful for phase- or frequency-modulated signals is to use a phase-locked loop. Design of PLLs is a complex topic in itself, but in essence, they are feedback systems can be used to acquire and track phase and frequency offset while your receiver is in operation. If you only need frequency synchronization, then you can use a frequency-locked loop instead; while they will not provide phase synchronization for you, they often have better acquisition properties.
As an alternative to a feedback loop, there are also feedforward approaches to estimating frequency or phase offset. One feedforward approach would take advantage of your training sequence to estimate the frequency error based on how the phase offset changes over the course of the sequence. If the frequency offset changes as time goes on, however, you would need to repeat the estimation procedure to allow your receiver to catch up.
Another technique is to design your system to be robust to (reasonably small) frequency offsets. Differentially-encoded phase modulation is an example of this (although frequency error will show up as a phase offset after differential decoding, which must be handled). Frequency-modulated waveforms such as FSK also have some level of resistance to frequency offset, as long as the offset is small relative to the amount of frequency deviation used by the transmitter.

This very short summary really just scratches the surface of some of the more well-known approaches. Synchronization can be a difficult problem to solve in practical ways, and there has been much research over the years on different ways to do it. It's going to depend upon exactly how your system is structured, and one very important variable: the target SNR. There isn't a single "correct" answer. I will make one textbook recommendation; while it is very expensive, "Synchronization Techniques for Digital Receivers" by Mengali is a comprehensive text on timing, phase, and frequency synchronization.


Your value for $f_s$ is missing. Assuming you work above Nyquist sampling rate (i.e. $f_s > 2*f_{sig}$), I obtain the following results:

fs = 20000;
len = 16384;
fsig = 1000;
Sig = 1.5 * sin( (1:len) * 2 * pi * fsig / fs);
fsig = 700;
Sig = Sig + 0.9 * sin( (1:len) * 2 * pi * fsig / fs);

[matWelchD1, matWelchF1] = pwelch(Sig, length(Sig), 0, [], fs);
[matWelchD2, matWelchF2] = pwelch(Sig, ceil(length(Sig)/2), 0, [], fs);

subplot(2,1,1);
plot(matWelchF1, (matWelchD1), '-o'); grid;

subplot(2,1,2);
plot(matWelchF2, (matWelchD2), '-o'); grid;

df1 = matWelchF1(2) - matWelchF1(1);
df2 = matWelchF2(2) - matWelchF2(1);

df1
df2

Output:

df1 =

    1.2207

df2 =

    2.4414

As you can see, both frequencies are nicely represented in the figures.

One more comments on the window length (i.e. second parameter to pwelch): It determines how big the sections are, into which the signal is divided. Each section is multiplied with a Hamming window and then the FFT is taken. Afterwards, the values of all FFTs are summed together, yielding the PSD estimate.

I.e. putting the signal length as the window length, would result in a single section. For your stationary signal this is fine, but in reality you would want to have a smaller value, or leave the parameter out (or use the empty vector []), such that the PSD can be more accurately estimated.

As you can see, different window length create a different amount of frequency samples: The longer the window, the more frequency samples. NOte that pwelch calculates the Power Spectral Density, i.e. the power per Hertz. In order to get the energy of a frequency you need to multiply with the bandwidth of each bin:

>> sum(matWelchD1) * (matWelchF1(2)-matWelchF1(1))

ans =

    1.5300

>> sum(matWelchD2) * (matWelchF2(2)-matWelchF2(1))

ans =

    1.5300

Both PSDs contain the same overall energy. Since the frequency samples in the first PSD are closer to each other, the PSD has higher peaks for each frequency (to deliver the same overall power).


Is there a means to identify/separate multiple out of phase copies of the same frequency?

It is not possible. Let's consider the case of two cosines, with angles $\theta = 2\pi f_0 t + \alpha$ and $\phi = 2\pi f_0 t + \beta$. Then their sum is equal to $$2\cos\left(\frac{\alpha-\beta}{2}\right)\cos\left(2\pi f_0 t + \frac{\alpha+\beta}{2} \right).$$

Note that, if you measure the amplitude and phase of this signal, you can find plausible candidates for $\alpha$ and $\beta$, but in general there will be multiple possible solutions. Note that the signal could also be zero! As you increase the number of original cosines at $f=f_0$, the ambiguity will also increase.

As a final note, the DFT does provide phase information. The GnuRadio tools, however, only plot the magnitude of the DFT, and there the phase information is lost.


Correlation of two signals at different sample rate

If you are looking for the exact signal, then yes you need to resample one of them to the same rate as the other before performing a cross correlation. Otherwise you will be cross correlating against a time-stretched (or compressed) version of the desired signal (which probably will not yield accurate results). If your $Y$ signal bandwidth is less than $4000\textrm{ Hz}$, then you may want to use an anti-aliasing filter and resample it to $8000\textrm{ Hz}$ for reduced computational load in the cross correlation. If not, then you must resample the signal $X$ to $44100\textrm{ Hz}$ in order for a cross correlation to be meaningful.


LTI system without constant coefficient differential equation

As long as a system can be characterized by its transfer function, or, equivalently, by its impulse response, the system is LTI. For the LTI property it's irrelevant if the transfer function is rational or not.

If the transfer function is rational (and positive-real), the system can be realized by a circuit containing lumped elements, i.e., ideal R, L, and Cs. On the other hand, a simple example of a non-rational transfer function is an ideal delay line (e.g., approximated by a long cable), the transfer function of which is $H(s)=e^{-sT}$, where $T$ is the delay.

In sum, the exact expression for $H(s)$ is irrelevant for the LTI property of the corresponding system. However, $H(s)$ can only be realized by resistors, inductors, and capacitors if it is a rational function (and if it is positive-real).


Laminar air flow doesn't produce much in the way of sound.  Turbulence, which is usually excited by interaction that is, the disrupted flow by an interacting object.  Some of what is heard, is radiation from monopoles (dipoles, and quadrupoles) within the flow and some are resonances related to the object doing the disruption.  You can often look at a tree in high wind and visually correlate the intensity of sound with the intensity of motion of leaves and limbs.  If you stand near a steel tower in high wind, your will hear it "sing"

Part of what you record is the broadband excited vibration of the microphone's full structure.  This is why a windscreen on microphone works.  The screen dampens the wind incident on the microphones structure and is a poor radiator.  One model is a random stream of impulses exciting a linear filter.  The mic structure gives color to the wind sound.

This is a situation where a little care in data collection saves a lot of trouble down the line.

I read someplace once "Never invent gloves to fix a leaky valve"

You might also want to try to whiten the wind noise and see what happens.  It would be similar to inverse filtering on the structure resonance.

Another characteristic of the vehicle sounds are the doppler shifts at the time when the car is at closest distance to the microphone.


How to approach detection of anomalies in a video feed?

Here is some brainstorming:

You have two interesting variables here, I would say:

spatial changes in y-direction (left to right)
temporal changes (t-direction).

I would not directly include the x-direction (up-down direction) in the first try, since this is correlated to the t-direction.

Try the following:

Calculate the y-gradient
Calculate the t-gradient
Create a binary mask from both using two threshold values (you probably have to figure them out in a test run), so that the whole image is zero where there is no change and one where there is a change.
To suppress spurious noise in the binary mask, you could median-filter it with a 5-by-5 or 7-by-7 median filter (try out the best size)
Use a morphologic dilate filter to close "holes" in the binary mask and also to dilate the masks so that they do have a broader overlap

This way you are directly using the visual features of the image - changes in the gray values. Now a few things can happen with those two masks

Areas of the t and y masks overlap: You have a change in the image in y direction and there is also a change between to consecutive images - you found a localized change that could be a hole or a speedbump
There is a change in y direction, but no or just a small area in t direction: That could be a road marking (the algorithm might fail if your drone tumbles quite strongly and your camera has a low number of frames per second - then the t-gradient might be quite high in this case. To suppress those signal changes from lane markings, you could suppress lines in the binary mask by a Hough-transformation of the mask and kicking out pairs of dots)
There is change in t but not in y - probably a lane mark that started/stopped

If this really works with your data, I do not know. However, starting with simple things can lead you into a quite robust and reasonably fast direction for your image processing task.


The idea of applying interpolation works well when you have a very high sampling rate. In some sense, what you are trying to do here is fit the missing data between the second and third green bars.

Approach

This approach is motivated from the results discussed in [1].

Construct gapped IMFs (Intrinsic Mode Functions)
Fill in gaps in each mode, based on geometrical constraints
Add up all gap-filled IMFs

Though, I am not sure if this can be implemented in real-time. You can make assumption that you are processing 100 ms data in every batch.

Reference

[1] A. Moghtaderi, P. Borgnat and P. Flandrin, "Gap-filling by the empirical mode decomposition," 2012 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Kyoto, 2012, pp. 3821-3824.


Predicting account balance based on historical data

Given the strong correlation of the downward trend in each cycle, I offer the following that would offer an improved prediction over one linear regression, but not be significantly more complicated (but still would like to see what a precise "best estimator" would be):

Do a linear regression to determine the general upward trend. This would be the Income - Spend = Savings rate.
Do a least mean square linear fit for all the downward trend cycles. This would be the spend rate. (my temptation to do this easily would be detect the upward spike and superimpose the datasets, find the mean downward trend and do a linear regression to that.)
Determine the time of each upward spike to determine the best estimate of the frequency of occurrence of the upward spikes (clock recovery). For this the strongest tone in the FFT should reveal this.  This would be the pay dates.
Determine the mean of the magnitude of the upward spikes.  This would be the income per pay date.

From the information above you can create a future projection of continued upward spikes at the best estimated times along an overall upward trend that will provide a better estimate than one linear regression over the whole dataset.


Shannon entropy is a way of measuring the degree of unexpectedness or unpredictability of a random variable. For example rolling a die has higher entropy than flipping a coin because the die has more possible outcomes making it harder to predict. Same goes for a biased coin versus a fair coin.


How to detect underdamped harmonic signal?

First, de-noise the waveform using a technique suitable for non-stationary signals. Admittedly, the spectrum of the signal doesn't change too much over time. There are efficient iterative de-noising algorithms for this purpose, such as basis pursuit de-noising [1]. But even a simple low-pass filter may be sufficient depending on how bad the non-stationarity is.
Since your record is mostly oscillatory, the underdamped signal effectively begins when the successive peaks of your signal are lower in amplitude. With your de-noised smoothed signal, attempt to pick up local extrema. This can be done in real-time by looking at the first derivative or discrete difference. Once you have extracted the extrema, using the values of the extrema, you can construct the signal envelope. Empirically, you can decide that if the envelope is decaying, then it signifies the starting point for under-damped behavior.

Fourier and short-time Fourier transforms may not be as effective in this problem since you are looking for the exponential decay feature, which is a temporal signature.

[1] Chen, Scott Shaobing, David L. Donoho, and Michael A. Saunders. "Atomic decomposition by basis pursuit." SIAM journal on scientific computing 20.1 (1998): 33-61.


How to implement anisotropic diffusion using java?

I coded a modest Java implementation of the Perona & Malik anisotropic diffusion, based on the ImgLib2 library. You can find it here.

You will have to import and learn how to use the ImgLib2 library though. Here is an example using the anisotropic diffusion algorithm.

Otherwise, here is a 2D only implementation, by Vladimir Pilny and Jiri Janacek, made for the ImageJ software.

And here is another one for ImageJ as well, by Christopher Mei.

Best regards.


What methods can be used to remove ringing artifacts in the output of a digital filter?

Sticking with linear systems, removing the ringing is nearly the same as adding back some of the spectral content that your really steep transition filters removed.  Why use some crazy scheme to add back the stuff in the "softer" transitions that your hard-edged filters cut out?  Just use a more reasonable total filter response in the first place.

Going to non-linear systems, you could use some sort of AI pattern matching to determine what kind of sound waveforms would be perceived by a human a ringing, and just gate those waveforms out.  But that might just add weird sounding artifacts of its own (as well as also spreading spectral content outside of your really steep filters).


With no filtering or pulse shaping what you have is your signal with the default rectangular pulse shape
$$
\Pi(t) = \begin{cases}1 & \text{for} & -\frac T2\leq t\leq \frac T2\\0 & \text{otherwise}\end{cases}
$$
Where $T$ is the symbol duration. The Fourier transform of the such a pulse result in a sinc function of the form
$$
H(f) = T\operatorname{sinc}\big((f-f_c)T\big) \equiv T\frac{\sin\big(\pi (f - f_c)T\big)}{\pi \left(f-f_c\right)T}\tag{1}
$$
The zeros in Equation $(1)$ are at
\begin{align}
(f - f_c)T &= k\quad \big\vert \quad k\in \mathbb{Z}_{\neq 0}\tag{with $f\neq f_c$}\\
\implies 2(f - f_c) &= 2kR_s\tag{2}
\end{align}
The sinc function in Equation $(1)$ dies slowly resulting in infinite frequency content. Hence the choice of other practical pulse shaping filters such as the SRRC for bandwidth conservation (and ISI mitigation). In brief, having rectangular pulses one has remnant lobes outside the bandwidth of interest whilst with practical ones such as the SRRC you get high attenuation outside the main bandwidth.
For bandpass modulated signals (e.g. PSK) the required double-sided bandwidth is defined as
$$
W_{\rm DSB} = (1 + \alpha)R_s\quad,\qquad 0 \leq \alpha \leq 1
$$
Where $\alpha$ is the roll-off of the raised-cosine pulse, a measure of excess bandwidth, and $R_s$ the symbol rate.
For bandpass systems, the null-to-null bandwidth is defined as
$$f_2 − f_1\tag{3}$$
where $f_2$ and $f_1$ are the first nulls of the magnitude spectrum $|H(f)|$ above and below $\arg\max_\limits f\left(|H(f)|\right)$ respectively. So, from Equation $(3)$ and $(2)$ you have that
\begin{align}
\arg\max_\limits f\left(|H(f)|\right) &= f_c\\
\left(f_2 - f_1\right) = 2(f - f_c) &\equiv 2R_s\tag{with $k = 1$}\\
\end{align}


Inter-Symbol Interference (ISI) channels

Things like the good ole' GSM standard specify power delay profiles and a phase variation model. So, there you go, reference model delivered in the shape of a standard.

Other than that, if you're just looking for any ISI-inducing (i.e. frequency-selective) channel, why not use your favourite DSP framework to come up with a random FIR filter?

Also, a lot of frameworks do come with their own channel simulations; GNU Radio for example brings multiple channel models of varying realism in its gr-channel module.

I often find that students overestimate the complexity of actually measuring a channel impulse response to use it in their own channel simulations. Just send a PN, receive it at the same time, and use a correlation to find the impulse response; if your receiver and transmitter exhibit a frequency offset, you might need to try a couple of frequency shifts and use the one that gives the highest correlation coefficients.


If $x(t)$ is a deterministic signal, then you have $E[x(t)]=x(t)$. Furthermore, if $Y(t)$ is a random process you have $E[x(t)Y(t)]=x(t)E[Y(t)]$.


Can FFT information be used to obtain the DTFT?

Your question is a bit difficult to comprehend. It seems that you are requesting a closed form solution to the DTFT given the DFT of a sequence. This is of course possible. If there is no time-domain nor frequency-domain aliasing, then we can simply apply the formula for the IDFT on the samples (or equivalently perform the IFFT). Following this, we can apply the formula for DTFT on the time-domain samples. This may not be an answer in the simplified form you want, but it is a closed form expression for the DTFT and can be evaluated at any point. Furthermore, a very good computer algebra system (e.g., Maple, Mathematica, etc.) may be able to reduce it to the form that you have requested (meeting your requirement for no computation by hand).

The DFT is a frequency-domain sampled version of the DTFT. The FFT computes the DFT. Because of this, you could simply apply the Shannon-Nyquist Sampling Theorem and interpolate the points of the DFT using sinc functions,

$$\DeclareMathOperator{\sinc}{sinc}
\begin{align} X_{DTFT}(f) &= \sum_{k=-\infty}^{\infty} X_{DFT}[k] \sinc(f - k) \\
&= \sum_{k=0}^{N-1} X_{DFT}[k] \frac{\sin(\pi(f-k))}{N \sin\left(\tfrac{\pi}{N}(f-k)\right)}
\end{align}$$

where $\sinc(x) = \frac{\sin(\pi x)}{\pi x}$. This expression is equal to your simplified expression (i.e., it evaluates to the same value at every point). This expression assumes a unit length sampling interval (i.e., $T = 1$). This at least suggests a method for interpolating the DTFT from the DFT samples.

Looking for an exact simplified expression is difficult because different schools of thought will simplify it in different ways. Who is to say that your expression involving trigonometric functions is any better than one involving complex exponentials? For example, isn't the expression
$$ X_{DTFT}(f) = \frac{1}{4} + \frac{1}{4} \exp(-j 2 \pi f) + \frac{1}{4} \exp(-j 4 \pi f) + \frac{1}{4} \exp(-j 8 \pi f)$$
exactly what you are after? This is exactly what we would get if we used the approach in the first paragraph. Take the IFFT of the DFT, and then apply the formula for the DTFT. A computer program could do all of this and could evaluate this expression at any point $f$, requiring no computations by hand.


How to obtain round peaks instead of triangular peaks when smoothing an input?

Use two or three passes of moving average (can be pipelined).


FFT-like algorithm for fast DTFT computation?

According to Oppenheim and Schafer's "Discrete Time Signal Processing", the Goertzel algorithm will be more efficient than the FFT in computing an N point DFT if less than $2 Log_2 N $ DFT coefficients are needed.  So if you need to compute 15 frequency points, the Goertzel will be more efficient for total number of samples $N > 65536$.  Anything below that and the FFT will be more efficient than Goertzel as Marcus has implied.

As you probably discovered, one point in the DFT requires N phase rotations and N complex additions (assuming a complex input).  (As the DFT output is computed by simply multiplying your input signal through successive rotations and accumulating the result-- in this manner you are correlating to each output frequency of interest).

Each phase rotation requires 4 multiplications and 2 additions (if you have time to iterate with little multipliers, see the CORDIC algorithm as an alternative phase rotator that requires no multipliers).  The phase rotations that are +/-j are simplified, since a rotation by j is done by simply changing the sign of Q and swapping I and Q, but for a larger FFT this will be an insignificant number of the total computation.

So therefore for each point using the DFT, there are approximately 4N total multiplications and 4N total additions.  So for an N pt DFT that would be $4N^2$ of each - wow!

In comparison the FFT requires $2Nlog_2N$ multiplications and $2Nlog_2N$ additions. When considering a complete DFT, the FFT offers a dramatic reduction in required computations. (Pause here for some respectful bowing to Cooley-Tukey).

I would be interested in seeing how the Goertzel algorithm lands on the graph below but haven't muddled through those details yet.

Regarding your other challenge mentioned that your frequencies do not line up; for that I recommend windowing first to reduce the effect of distant frequencies on each tone you are measuring and then interpolate between the bins to more accurately estimate your frequency result for each tone. (A simple linear interpolation may be sufficient).  The windowing will reduce the frequency resolution of each bin but will also reduce the sidelobes from the other bins, so you would need to confirm that your frequency separation for your tones of interest are sufficiently larger than the equivalent bandwidth of your window+FFT. For a good summary of equivalent bandwidths vs window refer to fred harris' paper on windowing: https://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf


How to determine the power of a finite signal?

First of all, note that your formula for signal power is only valid for real-valued $x(t)$ that satisfy $x(t)=0$ for $t<0$.

A more general formula is

$$P_x=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^{T}|x(t)|^2dt\tag{1}$$

For signals with finite energy

$$E_x=\int_{-\infty}^{\infty}|x(t)|^2dt\lt\infty\tag{2}$$

the power $P_x$ defined by $(1)$ is clearly zero.