consider an linear time-invariant system (LTI) in continuous time.  the input-output relationship is:

$$ y(t) = \int\limits_{-\infty}^{\infty} h(u) x(t-u) \, \mathrm{d}u $$

and, in the frequency domain

$$ Y(f) = H(f) X(f) $$

where

$$ X(f) \triangleq \mathscr{F} \big\{x(t)\big\} = \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \, \mathrm{d}t  $$

and similarly for $H(f)$ and $Y(f)$.

now, suppose dimensionally $x(t)$ and $y(t)$ are the same species of animal, like they are both voltage or both position or both pure numbers.  then $X(f)$ and $Y(f)$ are also the same dimension of physical quantity.  then, dimensionally, what must $H(f)$ be?

$H(f)$, which is

$$ H(f) \triangleq \mathscr{F} \big\{h(t)\big\} = \int\limits_{-\infty}^{\infty} h(t) \, e^{-j 2 \pi f t} \, \mathrm{d}t  $$

must be dimensionless.  $H(f)$ must be a pure number with no units attached.  that is the only way it can multiply $X(f)$ and result in $Y(f)$ and the latter two have the same units or dimension.

now, in the integral above, what dimension must $h(t)$ have for the resulting integral to be the dimensionless $H(f)$?  consider the relationship between velocity and position:

$$ v(t) = \frac{\mathrm{d}r}{\mathrm{d}t} $$

or

$$ r(t) = r(t_0) + \int\limits_{t_0}^{t} v(u) \, \mathrm{d}u $$

$u$ is a dummy variable of integration and is time, just like $t$ is.
dimensionally, what does the little differential element $\mathrm{d}u$ have?  $v(u)$ is velocity, length divided by time, and it is multiplied by tiny time increments, $\mathrm{d}u$, to become tiny length increments.  those tiny length increments are added up to become the length between positions $r(t_0)$ and $r(t)$.

so in the integral $\int_{-\infty}^{\infty} h(t) \, e^{-j 2 \pi f t} \, \mathrm{d}t$ , what must $h(t)$ be, dimensionally, so that the result of the integral is the dimensionless $H(f)$?

and what is the dimension of $h[n]$?

so this is the first step, understanding what must happen dimensionally to $h(t)$ to scale it to $h[n]$.

perhaps the best way to see why the dimensional scaler should be the sampling period, $T$ is by viewing the convolution integral as a Riemann summation:

$$\begin{align}
 y(t) &= \int\limits_{-\infty}^{\infty} h(u) \, x(t-u) \, \mathrm{d}u \\
\\
  &\approx \sum\limits_{i=-\infty}^{\infty} h(i \Delta u) \, x(t-i \Delta u) \, \Delta u \\
\end{align}$$

$\Delta u$ is a vanishingly small period of time.  If the sample rate $\frac1T$ is sufficiently high or $T$ is sufficiently small that it can take the place of $\Delta u$, then the Riemann summation is exactly the same as the discrete-convolution summation:

$$\begin{align}
 y(t) &= \int\limits_{-\infty}^{\infty} h(u) \, x(t-u) \, \mathrm{d}u \\
\\
  &\approx \sum\limits_{i=-\infty}^{\infty} h(i \Delta u) \, x(t-i \Delta u) \, \Delta u \\
\\
 \to y(t) &= \sum\limits_{i=-\infty}^{\infty} h(i T) \, x(t - iT) \, T \\
\\
 \to y(nT) &= \sum\limits_{i=-\infty}^{\infty} h(i T) \, x(nT - iT) \, T \\
\\
  &= \sum\limits_{i=-\infty}^{\infty} T h(i T) \, x\big((n - i)T \big)  \\
\\
 \to y[n] &= \sum\limits_{i=-\infty}^{\infty} \Big( T h(i T) \Big) \, x[n-i]  \\
\\
  &= \sum\limits_{i=-\infty}^{\infty} h[i] \, x[n-i]  \\
\end{align}$$

there you can see (changing the dummy variable $i$ to $n$) that

$$h[n] = T \, h(nT)\ .$$


Is it okay to interpolate between overlapping areas obtained using a sliding window?

A first suggestion is to use an odd-sized sliding window (say $w=51$), so that you can really affect the black/white proportion to the center pixel. An even length can generate a half-pixel bias.
And you have a  window sliding by one pixel, therefore you have a black/white proportion for each pixel (except on the borders of the large image, of dimension $W\times W$). The resulting picture will be of size $(W-w+1)^2$, about $1950\times 1950$ here, which already seems to have a large definition for most displays. Hence, I cannot really understand the need to interpolate for visual purposes.
Nevertheless, two contiguous windows share $w\times (w-1)$ (sides) or $(w-1)\times (w-1)$ (diagonal) pixels in common, so more than 98% or 96 % respectively. Hence, the variations in the proportion should vary quite slowly. Then simple interpolation should work nicely. The problem is to give a meaning to the interpolated values. It should be close to having a window with fractional offsets over the big image, and counting partial pixels.


Can anyone explain how does CZT (Chirp Z Transform) really help in 'spectral zooming'?

The CZT allows for a fairly general evaluation of the Z transform - the more general evaluation path looks like a spiral, so it has a radial component step size as well an angular step size.For  spectral zooming, you're only using a subset of this. You're evaluating around the unit circle and only for a small set of frequencies.

The Zoom FFT can be implement in a more traditional manner:

Frequency shift by multiplying time domain signal by a complex exponential - baseband around the center frequency.
Low-pass filter and decimate - to reduce the sampling frequency.
Apply a (shorter) FFT

Note that the FFT in step 3 is shorter because of the reduced sampling rate. To get a frequency resolution of $\nabla f$ you need a time series of length $T=1/\nabla f$ and at a lower sampling frequency you can get that time length signal in fewer samples. This series of steps applies a series of operations to the signal.

Now another way of looking at those series of steps is to apply them to the filter rather than the signal. Now you're shifting a low-pass filter in frequency so it becomes a complex band-pass filter and the decimation step causes duplicates of your band-pass signal to appear around 0 Hz (and throughout the entire frequency spectrum). At this point you can apply a shorter FFT. To make the two series of steps exactly equivalent there are some restrictions on the sampling frequency - If you shift the filter to a center frequency of $f_c$ then you have to be careful of the decimation factor if you want the frequency content at $f_c$ to be mapped exactly to 0 Hz.

The CZT for Zoom FFT is essentially doing these same steps. If you examine the CZT from this point of view you can work the mathematical equivalence.


This appears similar to a classic least squared solution of an overdetermined equation that proceeds as follows:

Starting with:

$$ \mathbf{x} = \mathbf{A} \mathbf{s}  $$

$\mathbf{A}$ is not a square matrix if overdetermined (more equations than unknownns) so therefore an inverse does not exist. What you do then is multiply both sides by the transpose of $\mathbf{A}$ since $\mathbf{A^TA}$ is a square matrix. From that, assuming and inverse exists for$\mathbf{A^TA}$, you can solve for s:

$$ \mathbf{A^Tx} = \mathbf{A^TA} \mathbf{s}  $$

$$ \mathbf{s} =  (\mathbf{A^TA})^{-1}\mathbf{A^Tx} $$


Sampling the Wave Equation for String Vibration

Assuming I understand the question correctly, you need to sample the string's position at least twice as fast as the maximum frequency you anticipate the string can move or the highest frequency you care about. So if the highest overtones the string produces are 16kHz, you need to measure at at least 32kHz if you want to be able to reconstruct the original motion faithfully. This is the Nyquist-Shannon sampling theorem and you can probably find several hundred descriptions of it here and elsewhere on the internet.
If you know the lowest frequency the string can produce or that you care about, there are other tricks you can do (look up heterodyne) to reduce your sample rate requirements assuming you have the equipment to do so on the analog side.


Realization of a filter based on its transfer function

A transfer function is called realizable if it can be implemented by a causal and stable system. The given frequency response is continuous and doesn't have any impulses, so the corresponding system is stable.

The transfer function (as a function of $s$) is given by

$$H(s)=\frac{a}{b-s^2}e^{-st_0}\tag{1}$$

The given frequency response is obtained form $(1)$ by substituting $s=j\omega$. Since $b>0$, $H(s)$ has two real-valued poles at $\pm\sqrt{b}$, and, consequently, the corresponding impulse response is two-sided (i.e., non-causal).

Partial fraction expansion of $(1)$ gives

$$H(s)=\frac{a}{2\sqrt{b}}\left[\frac{1}{\sqrt{b}-s}+\frac{1}{\sqrt{b}+s}\right]e^{-st_0},\qquad -\sqrt{b}<\textrm{Re}\{s\}<\sqrt{b}\tag{2}$$

From $(2)$ it is straightforward to obtain the impulse response

$$\begin{align}h(t)=\mathcal{L}^{-1}\{H(s)\}&=\frac{a}{2\sqrt{b}}\left[e^{\sqrt{b}(t-t_0)}u(-(t-t_0))+e^{-\sqrt{b}(t-t_0)}u(t-t_0)\right]\\&=\frac{a}{2\sqrt{b}}e^{-\sqrt{b}|t-t_0|}\tag{3}\end{align}$$

No matter how large you choose $t_0$, the impulse response will never be causal. However, since the system is stable, the impulse response decays for $|t|\to\infty$, so you can approximate the actual impulse response well by sufficiently shifting it to the right (i.e., by choosing some large $t_0$), and truncating the non-causal part.

Also take a look at this related question and its answer.


number of samples and Amplitude accuracy,what's the solution to this problem

I think that what is happening to u is that the frequency which you are trying to measure is in between 2 of your FFT bins which will cause the power to spread between two different indexes of your result FFT.

Another possible solution , only if with higher resolution you are getting a lower value  , is that since as I understand this is a real signal with noise you may be summing more noise , on the first attempt, for each frequency of your FFT and as a result raising all the values of your FFT.

Without knowing the signal you are trying to process it is hard to know what is exactly the problem. also your FFT plots will help give a better answer.


What's the relation between the sampling rate and FFT size in OFDM

Yes, you can use different sample rates and corresponding FFT sizes to accomplish the same thing for OFDM (mod and demod). For example, 10 MHz LTE can be sampled with 15.36 Msps and modulated/demodulated with a 1024-point FFT for 15 kHz subcarrier spacing and  the 600 subcarriers around DC give the 9 MHz of utilized BW. Alternatively, you can use 30.72 Msps and 2048-point FFT to mod/demod and get the same 600 subcarriers.

What is effectively happening is that the basis vectors are resampled (by 2x in this case), but they remain orthogonal and at the same subcarrier spacing.


Alright, the way to think about the rectangular window (as well as all of the other windows) is that you begin with an infinite sequence of samples and you multiply the samples outside the window with 0 and the samples within the window with something non-zero (and usually, but not always, positive).

So the DFT (of which the FFT is fast method of doing the DFT) does periodically extend the $N$ samples passed to it.  But there are a couple of different ways of looking at it.  One way to look at it is with just the definitions of the DFT and iDFT:

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

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

The $x[n]$ that the DFT $X[k]$ is aware of is simply a periodic sequence with period of $N$. both $$ x[n+N]=x[n] $$ and $$ X[k+N]=X[k] $$ are true.

Now the alternative way to think about the DFT is that it is the sampling of the spectrum of the Discrete-Time Fourier Transform (DTFT):

$$ \hat{X}\left( e^{j \omega} \right) = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] \ e^{-j \omega n}  $$

where $ \hat{x}[n] =
\begin{cases}
0, & n < 0 \\
x[n], & \text{for }0 \le n < N \\
0, & N \le n
\end{cases} $

and this is the same as multiplying by a rectangular window.

$$ \hat{x}[n] \ = \ x[n] \ w[n]  $$

where $ w[n] =
\begin{cases}
0, & n < 0 \\
1, & \text{for }0 \le n < N \quad \quad \\
0, & N \le n
\end{cases} $  (length-$N$ rectangular window)

But here you cannot tell that you extracted $\hat{x}[n]$ from a periodic $x[n]$ or if those $N$ samples came from an infinitely long and non-repeating stream of samples.  Whatever it is, $\hat{x}[n]$ is non-repeating, infinite in length (with all those zeros padded), and you can apply the DTFT to it.

It doesn't have to be the rectangular window, it could be any length-$N$, finite-valued window for $w[n]$ and the math that follows is the same.

Now with that definition of $\hat{x}[n]$ and  $\hat{X}\left( e^{j \omega} \right)$, it's not hard to show that

$$ \hat{X}\left( e^{j \omega} \right) \Bigg|_{\omega = 2 \pi \frac{k}{N}} = \hat{X}\left( e^{j 2 \pi k/N} \right) = X[k] $$

Now, what happens in either domain, is that when you uniformly sample the function (say, in this case, the spectrum $\hat{X}\left( e^{j \omega} \right)$)  in one domain, it causes in the other domain, the signal ($\hat{x}[n]$) to be repeatedly shifted (by $N$ samples) and overlap-added:

$$ x[n] = \sum\limits_{m=-\infty}^{+\infty} \hat{x}[n+mN] $$

It is here where the discontinuity between $x[N-1]$ and $x[N]$ (which is the same as $x[0]$) becomes a concern, and tapered windowing to reduce or eliminate this discontinuity is perhaps indicated.


Random Signals - statistical properties are time dependant?

The answer to the question (a counterexample)

Properties of random processes will in general be time-dependent. They are not only when talking about stationary processes. Another related concept (not relevant here but that you might find interesting) is ergodicity. Most of us find it difficult to understand the difference between both concepts, but then you find this brilliant Dilip's answer to clear out all doubts.

Going back to your question: take for example the simple random process

$$X(t)=At$$

where $A$ is some random variable with mean $\mu_a$. The mean of the process would then be:

$$\mathbb{E}[X(t)]=\mathbb{E}[At]=t\mathbb{E}[A]=t\mu_a$$

It's easy to see that the mean depends on $t$, thus it's time-dependent. I hope this easy example shows why, in general, properties of random processes vary with time.

Additional info that might help

I think your confusion arises from thinking that to build a stochastic process you grab many realizations of a single random variable and put them in a sequence to represent a time-signal. That's not how a random process is built in general. If you do this, then the mean will not vary with time because the process will be stationary. But if you don't, it may or may not vary.

I'll use an example here to make myself clear. Suppose you have a random variable $Z$ that represents the result of throwing a fair 6-faced dice. Then $\mu_z = 3.5$. If you throw the dice 3 times and you put them together to form a discrete-time signal of 3 samples, then you will have the following random process:

$$X(n)=Z_1\delta(n) + Z_2\delta(n-1) + Z_3\delta(n-2)$$

where $Z_i$ represents the random variable corresponding to the $i$-th throw of the dice.

Note that each sample of $X(n)$ is a random variable. This is correct for any random process in general. For a given random process $X(n)$, $X(n_0)$ for any $n_0$ will be a random variable.

On the other hand, once you've thrown the dice 3 times, you are not dealing with anything random anymore: you have a signal. That is, a realization of the random process.

If you build a random process such as we have just built $X(n)$, then you are right: the mean of the process will be constant and equal to the mean of the random variable used for each sample (as in this case). In fact, all properties will be constant and equal to those of the variable. But that's just in that specific case, because all samples of the process will be i.i.d., and thus, as @Fat32 has explained in the comments, the random process will be independent of time ($t$). In our dice example:

$$\begin{align}
\mathbb{E}[X(n)]&=\mathbb{E}[Z_1\delta(n) + Z_2\delta(n-1) + Z_3\delta(n-2)]\\
&=\mathbb{E}[Z_1]\delta(n) + \mathbb{E}[Z_2]\delta(n-1) + \mathbb{E}[Z_3]\delta(n-2)] \\
& = \mu_z \qquad \text{for any} \  n\in\{0,1,2\}
\end{align}$$

But this will not always be the case with random processes, such as in the example I wrote at the top of the answer, in which every random variable obtained by evaluating the random process at different values of $t$ is different from all others, therefore preventing the process from being i.i.d. and making its statistical properties time-dependent.


The “O” in OFDM stands for orthogonal.  The IFFT takes a complex number (1 bin of the input) and turns it into samples of a sinusoid of a frequency that, over a certain length (of time), is orthogonal to (and thus, under ideal linear conditions, won’t interfere with) any other frequency subcarrier output by the IFFT.

If you just transmit the original bits, the impulse response of the channel can smear the bits together, interfering with the adjacent bits, and leave the decoder/demodulator the interesting job of unscrambling this mess.  Depending on the channel response and how it changes over time, this may or may not be more difficult to do.


How to interpret the signal in frequency domain

These questions will probably seem silly to you, sorry for that.

No question is "silly". My perception is that the questions are trivial and can be answered by looking things up in a textbook but if this is to get you out of a "difficult corner", here it goes:

  ...I was expecting that the real part will represent amplitudes of two sinusoidal signals and the imaginary part will represent a phase shift of them.

The real and imaginary part of the complex Discrete Fourier Transform are one number which carries many pieces of information. Specifically, the magnitude $|X[k]|$ of some $k$ complex sinusoid returns information about its amplitude.

$$|X[k]| = \sqrt{\mathcal{R}(X[k])^2 + \mathcal{I}(X[k])^2}$$

Or in more practical terms Y_amp = np.abs(np.fft.fft(y)) following the signals established in your example.

The "phase shift" is the angle $\angle X[k]$ of a given $k$ complex sinusoid.

$$\angle X[k] = \arctan\left(\frac{\mathcal{I}(X[k])}{\mathcal{R}(X[k])}\right)$$

Or in more practical terms Y_phs = np.angle(np.fft.fft(y)).

Where $\mathcal{I},\mathcal{R}$ denote the imaginary and real part of a complex sinusoid at discrete frequency $k$.

  However the amplitude 109 doesn't make much sense to me, nor does the phase shift of +/-128.

That's because of conventions used in the evaluation of the DFT integral. So, if you look at the documentation of the FFT, you will notice that there is a $\frac{1}{N}$ term (where $N$ is the length of the signal) that can be applied either on $x[n]$ before the transform, or on the $\sum$ at the inverse transform. In this case, the sum is evaluated without scaling in the forward transform and with scaling when performing the inverse transform. If you divide your amplitude spectrum values by 256, your strengths across the spectrum will be adjusted to the values you probably expect. Again, in practical terms that would be: Y_amp_normed = np.abs(np.fft.fft(y))/256.

And since we are on the topic of phase, sooner or later you are going to need unwrap to deal with phase wraping.

  How exactly can I figure out which frequency represents one peak or another?

$k$ denotes discrete frequency. To convert between discrete frequency and natural frequency use:

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

For more information about this you can have a quick look at this link

Hope this helps.


Below is a function which I wrote long back, when I needed to generate AWGN time-domain samples given Noise PSD in dBm/Hz.

AWGN_NOISE() : Generates Additive White Gaussian Noise of PSD power in dBm/Hz

AWGN has Gaussian PDF with 0 mean and $\sigma^{2} = N_{o}/2$

$NoisePSD_{dBm/Hz} = 10.log_{10}(\frac{N_o}{2.BW})$, Why?

Because, Output Noise Power(in dBm) : $$10.log_{10}(\frac{N_o}{2}) = NoisePSD_{dBm/Hz} + 10.log10(BW)$$

Hence,
$$\sigma = \sqrt{\frac{N_o}{2}} = \sqrt{BW.10^{\frac{NoisePSD_{dBm/Hz}}{10}}}$$

So, the python function would be as follows:

def awgn_noise(length, power, Bandwidth):
    sigma = np.sqrt(Bandwidth * 10**(power/10))
    noise = np.random.normal(0, sigma, length)
    return noise

Hope this answers the question.


Image Contrast Enhancement Threshold Selection

What you are looking for called Otsu Method / Otsu Thresholding.
It solves an optimization problem to set the optimal threshold between two modal distribution.


Transform coding gain of two channel Haar Wavelet transform (first stage)

I am not sure why you're trying to do this in MATLAB. This can be done quite easily with pen and papers.

Another way of seeing what you're doing here is the following:

Build vectors of size 2 from the sequence $X$
$$[X(0),~X(1)], [X(2),~X(3)], [X(4),~X(5)]\dots$$
Apply an orthogonal transform on each vector to obtain the corresponding transform coefficients $Y$.

The order 2 orthogonal transform you're using can mathematically be written as
$$\begin{bmatrix}Y_0 \\ Y_1\end{bmatrix} = \frac{1}{\sqrt{2}}\begin{bmatrix}1 & 1 \\ 1 & -1\end{bmatrix}\begin{bmatrix}X(n) \\ X(n+1)\end{bmatrix}$$
where $Y_0$ and $Y_1$ are the transform coefficients and $X(n)$ and $X(n+1)$ are successive (hence correlated) samples from the process to be compressed. From this, you can easily compute $\sigma^2_0$ (the variance of $Y_0$) and $\sigma^2_1$ (the variance of $Y_1$). Indeed, because $X$ is zero-mean (and so are the transform coefficients), you know that
$$\sigma^2_i = R_{Y_i}(0) = E\{Y_iY^*_i\}$$
for $i = 1, 2$.
Then, you simply need to plug in the expression of $Y_i$ as a function of $X(n)$ and $X(n+1)$ to find $\sigma^2_i$ as a function of $R_X(0)$ and $R_X(\pm1)$. At the end, you should obtain
$$\sigma^2_0 = R_X(0) + R_X(1) = 19.3,$$
$$\sigma^2_1 = R_X(0) - R_X(1) = 0.7.$$
You can clearly see from this the energy compaction in the first transform coefficient, which will contributes to the coding gain. You now have everything you need to compute the coding gain (note that $N$ in your expression for the coding gain is the number of transform coefficients, 2 in this case).


How do you reduce $H\left(e^{j\frac{\pi}{2}}\right)$ further according to a textbook solution

By plotting the unit circle you will find that $1-e^{-j\pi/2} = \sqrt{2}e^{j\pi/4}$.


Output of an FIR filter has a peak at the beginning of the signal

It's the initial transient at the filter output due to the implied zeros of the input signal for $n$ less than $0$.

Let me simulate the effect here, since you do not post the actual data I will model it as close as necessary, with the following Matlab/Octave code:

Fs = 250;                    % Sampling frequency in Hz.
Ts=1/Fs;                     % Sampling period in seconds.
t = [0:Ts:4];                % Simulation interval for the signals.
L = length(t);               % Simulation signal lengths in samples.

x = 8E4 + 1E4*exp(-t)  + 1E3*filter(fir1(43,0.03),1,randn(1,L));
% Note: x[n] is roughly modeling your input signal sampled at 250 Hz.
% I have not included those spikes, as they're irrelevant for this analysis

b = fir1(65,[10/125 20/125],'bandpass');   % a BPF with passband [10,20] Hz.
y = filter(b,1,x);           % Obtain the filtered output signal.

% Now, observe the outputs...
figure,subplot(2,1,1);
plot(t,x);title('Input signal x[n] of duration 4 seconds')
subplot(2,1,2)
plot(t,y);title('Output signal y[n] obtained from BPFing of the input');

As it can be seen below:

So the peak at the beginnig of the output is clearly observed.

It's intuitive to look at the bandpass filter's impulse response :

Looking at the BPF impulse response (especially the second subplot) reveals the fact that the transient is inline with the impulse response; their durations are the same and their shapes are very similar.

In fact, looking at the transient's duration (which is about 0.25 seconds) it matches with the number of samples the impulse response has:
$ N \times T_s = 0.25 s \rightarrow N = 0.25 \times 250 = 62.5 \text{ samples}$ This coincides with the BPF impulse response which was designed as a $65^{th}$ order filter having $66$ coeffcients.

So we have the clue that the peaking at the begining is indeed the impulse response (a transformed version of it) at the output.

The reason of this can be understood by looking at the convolution sum:
$$y[n] = \sum_{k=0}^{L-1} x[k]h[n-k] $$ which shall be evaluated for $n \geq 0$

now in the sum it can be seen that signal values for $k < 0$ are implied to be zero. Lets show this effect on the flip and drag graphical evaluation below:

Now first obervse that in its initial part, the signal $x[n]$ has an almost constant value. Then it can also be seen that until $n=65$ the overlap between the signal and impulse response is partial and therefore there won't be enough cancellation from the positive and negative samples of the impulse response to produce (apprx) zero output. But when $n \geq 65$ there is full overlapp between the impulse response and the signal (which is almost constant) therefore in the summation the positive and negative samples of the impulse response will cancel each other yielding a small, (zero-dc), output which represents the expected output after $n > 65$.


I wasn't familiar with that term in the context of signal processing. (Instead, I've seen the term being used in the context of the Riemann zeta function.) But I've found a document and this book where the term is used in a DSP context. The (obvious) definition is that trivial poles and zeros are the ones at the origin $z=0$ and at infinity $|z|=\infty$. They're called trivial because they don't affect the magnitude of the corresponding frequency response. Multiplying a given $\mathcal{Z}$-transform with $z$ (i.e. adding a pole at infinity and a zero at the origin) just advances the corresponding sequence by one sample, and multiplying by $z^{-1}$ (i.e., adding a pole at the origin and a zero at infinity) delays the corresponding sequence by one sample.

As an example, note that a causal FIR filter has as many poles as zeros, but all of them are at the origin $z=0$, i.e., they can't be used to shape the magnitude response. The same is true for an all-pole filter: all its zeros are at the origin, so they don't help to create a stop band (which would be one possible function of non-trivial zeros).

I think it's important to point out that the term trivial zero (or pole) is not a standard term used in DSP - as far as I know - but it appears to be an idiosyncratic use of a few authors.


Should I pass Kalman Filter absolute or offset-from-mean sensor values?

As mentioned in another answer, the variance of the measurement noise is a property of the model which is baked in to the Kalman filter.
I think you're running into a different problem, which is that your model is different from reality.  It sounds like you're modeling the system assuming that all the microphones are 'live'.  If that's the case, then the filter won't take a dead microphone as something to be ignored -- it'll take a dead microphone as a live microphone that's 'hearing' utterly valid silence.
Note that Kalman filters can do some apparently absurd things in the presence of a model mismatch.  Because they assume Gaussian distributions for everything, and because the PDF of a Gaussian drops off so profoundly as you stray from the mean, they really don't do a good job handling real-world processes that have low kurtosis, nor do they handle model mismatch very well.
You need to modify your filter, or sanitize its input.  It's probably best to modify the filter by augmenting it with an external "dead microphone" detector (and make sure you get it right -- you're basically automatically throwing out outliers, and there's known problems with doing that too enthusiastically).  Or, make sure a human goes through the data and makes sure the microphones are live.


A-weighted Digital Filter for 1/3rd Octave Band Analysis


Noise modeling and reduction in Video when we have just noisy video

Reduction of Noise in a video depends on what type of noise is present in the video.
The noise can be reduced by filters which are all present [There are software for this] and many type of filters such as median ,mean,wiener filter.
And The noise model can be build to add noise to the video there are many matlab code as well as software also there.
Gaussian ,salt and pepper,poisson and speckle noise model can be build and by adding the image sequence with the video you will get a noisy video.

imnoise(video,noisename);

this will give noisy video as output.

If you have any doubt ask me.


Lag in financial trading vs Group Delay in signal processing

i'm gonna work on this using the EE DSP notation i think is most conventional in the EE DSP field.  your exponentially-decaying 1-pole filter has this simple difference equation:

$$ y[n] = (1-p) \, x[n] \ + \ p \, y[n-1] $$

turns out (as we shall see below) that $p$ is the pole.  applying the Z transform gets you

$$ Y(z) = (1-p) \, X(z) \ + \ p \, z^{-1}Y(z) $$

from which we can solve for the transfer function:

$$ H(z) \triangleq \frac{Y(z)}{X(z)} = \frac{1-p}{1-p\,z^{-1}} = (1-p)\frac{z}{z-p} $$

so there is a zero at $z=0$ and a (real) pole at $z=p$.  the zero would go away if we took our input from $x[n-1]$ instead of $x[n]$ which would create a pole at $z=0$ and cancel the zero, but all that does is introduce another sample period of delay.  i don't think you wanna do that.

the (Kronecker) impulse response is

$$ \begin{align}
h[n] & = (1-p) p^n u[n] \\
& = (1-p) e^{\log(p) \, n} u[n] \\
& = (1-p) e^{n / \tau} u[n] \\
\end{align} $$

where u[n] is the Heaviside unit step function

$$ u[n] \triangleq
\begin{cases}
1, & \text{if }n \ge 0 \\
0, & \text{if }n < 0
\end{cases} $$

and $\tau$ is the "time constant" of the decaying exponential.  in order for the filter to be stable, the exponential must be decaying which means that $\tau > 0$, $\log(p) < 0$, and $0 < p < 1$

okay, so now the complex frequency response of $H(z)$ is

$$ \begin{align}
H(z)\bigg|_{z=e^{j\omega}} = H(e^{j\omega}) & = \frac{1-p}{1-p \, e^{-j\omega}} \\
& = \frac{1-p}{1-p \, \left(\cos(\omega) + j\sin(\omega) \right)} \\
\end{align} $$

the magnitude response is

$$ \begin{align}
|H(e^{j\omega})| & = \left|\frac{1-p}{1-p \, e^{-j\omega}} \right| \\
& = \left| \frac{1-p}{1-p \, \left(\cos(\omega) + j\sin(\omega) \right)} \right| \\
& = \frac{1-p}{\sqrt{ (1-p\cos(\omega))^2 + (p\sin(\omega))^2 }} \\
& = \frac{1-p}{\sqrt{ 1 + p^2 - 2 p \cos(\omega) }} \\
\end{align} $$

this is pretty familiar.  less oft observed is the phase response

$$ \begin{align}
\phi(\omega) \triangleq \arg \{ H(e^{j\omega}) \} & = \arg \left\{ \frac{1-p}{1-p \, e^{-j\omega}} \right\} \\
& = \arg \{ 1-p \} - \arg \{1-p \, e^{-j\omega} \} \\
& = 0 - \arg \{1-p (\cos(\omega) + j \sin(\omega)) \} \\
& = - \arctan \left( \frac{p \sin(\omega)}{1-p \cos(\omega)} \right) \\
\end{align} $$

the phase delay is

$$ \tau_\phi(\omega) \triangleq -\frac{\phi(\omega)}{\omega} $$

and the group delay is

$$ \tau_g(\omega) \triangleq -\frac{d\phi(\omega)}{d\omega} $$

gotta go.  i'll post this now and finish it later.  perhaps someone else wants to finish it.


There are two mistakes in your code/method. The first is the term $\sqrt{\Delta t}$ in your second formula; it should be replaced by $\Delta t$. The second is in the computation of the  power spectrum from the estimated auto-correlation. What you do is square the result of the FFT Y to obtain mY, but that's not correct. First of all, Y is complex-valued, and second, the Fourier transform of the auto-correlation is the power spectrum, so there's no need to square the outcome of the FFT, you just need to compute its magnitude.

In sum, replace the line

v*sqrt(sigma)*sqrt(dt).*normrnd(0,1);

by

v*sqrt(sigma)*dt.*normrnd(0,1);

and

mY = (Y).^2/N;

by

mY = abs(Y)/N;

This should give you the desired result, as shown in the figure below:


Multichannel input/output evaluation board

This is a little bit of a tricky question because it is homework and a direct answer would simply equate to doing that homework which is not what this website is about. Apart, from the purely technical bit of the right decisions, which you should do entirely on your own, I am interpreting this question as "There are so many options out there! Where do I start? I am lost!" and I am providing the following answer on a general framework that can be used in these cases.

Very briefly: Use a requirements matrix to rank your options.

The requirements matrix lists each of the requirements of the project. These are already available (Minimum sampling frequency, C compiler, number and type of I/O, real-time requirements, etc). Some of the requirements of course will need further unpicking. For example, the requirement "Capable of doing sound localisation" immediately leads to the question of "How accurate does this have to be?" which tightens the requirements of sampling frequency and data processing. Is floating point processing required or fixed point? Is the board fast enough to emulate floating point? And other similar considerations.

For each design option, add one column and provide a simple answer of how well does it cover that requirement. For example, suppose that you have an audio design project and your possible implementation options are Product 1, product 2, product 3, ... and each one has different ways of providing an audio output (Simple filtered Pulse Width Modulation, audio chipset, virtual sound card, etc).

At the very minimum, this matrix will help in selecting the one option that covers the most requirements for the design. More elaborate scoring scales are possible but also application dependant.

For more information, please see this link.

Hope this helps.


What is the Implementation of a Low, Mid or High Tone Dodge or Burn in GIMP?

It basically using what Photographers calls "Shadow Mask" / "Highlight Mask".

If the image is in the range [0, 1] "Highlight Mask" is actually the value of the pixel.
"Shadow Mask" is the inverse, namely 1 - Highlight Mask.

Use this mask as linear interpolator of the result with the originals.


Understading sampling theorem and aliasing

An ideal low pass filter has perfect attenuation above the cut off frequency.  So any aliasing distortion in that band will be filtered out.  In your example, this means any content above 3kHz.

For this application fmax can be made to extend up to the sample rate less the cutoff frequency of the low pass, or 21kHz.


How to construct ideal low pass filter in matlab?

The ideal lowpass filter is an infinitely long sinc function. It's Fourier transform is a rectangular shape as shown in your frequency spectrum diagram. In practice you have to window (truncate) it to a certain number of samples. The periodic width of the lobes of the sinc will correspond to the width of your frequency rectangle (lowpass cutoff frequency)


Technique to detect pitch after analysing only one pitch period?

in order to estimate the period of a quasiperiodic signal, normally you need an entire period plus a little more.  you need to compare a little snippet of the waveform with its counterpart in the adjacent period.  that's the fastest i expect a pitch detector to come up with anything.


Formant Analysis Frequency - What makes a peak value?

The spectra on the left are two line spectra of two different periodic functions, the top with fundamental frequency of 100 Hz and the bottom has $f_0$ of 200 Hz.  While they are not perfectly flat spectral lines, they are broadbanded with no apparent resonance.

It's a gross approximation, but I believe these two broadbanded periodic signals are meant to represent what comes out of one's voicebox in the throat.  It's the same broadbanded, except the $f_0$ depends on the pitch that one intends to sing (or speak) at.  If it were A above middle C (normally called "A4" or MIDI note 69) then $f_0$=440 Hz.

Now that broadbanded (same spectral envelope, but different $f_0$) driving signal on the left gets passed through an acoustic filter in the middle which is not dependent on $f_0$ but is dependent on the shape of the mouth, position of tongue, amount of nasal coupling, all which is dependent on the vowel one is singing.  In the frequency domain, the spectrum of the driving signal is multiplied by the frequency response to become the spectrum of the output.

The cause of the formant frequencies is due to the acoustic filtering of the vocal cavity that is one's head, mouth, lips, and nasel coupling.

A few of us have developed a pitch shifting algorithm that shifts pitch, but does not shift the formant frequencies.  I did this a quarter century ago and wrote a short paper about it.  You can sorta think of your voicebox as emitting a quasiperiodic string of impulses (or approximations to impulses) from the glottal closure and opening (with fundamental frequency $f_0$) and those impulses being convolved with the impulse response of the formant filter.  If the pitch shifting does not stretch or scrunch that impulse response, the formants remain unchanged.

But simple speeding up or slowing down the recording does scrunch or stretch the impulse response, which causes the formants to move.  That causes someone to sound like Alvin the Chipmunk (pitching up) or Darth Vader or Satan from hell (pitching down).


FIR filter reverse engineering from step response

The relation of the discrete-time step response, $s[n]$, and discrete-time impulse response, $h[n]$, is

$$ s[n] = \sum\limits^{n}_{m=-\infty} h[m] $$

and

$$ h[n] = s[n] - s[n-1] $$

from $h[n]$ you can compute the DFT (of a large finite segment of $h[n]$) and get the frequency response of this filter:

$$ H(e^{j \omega}) = \sum\limits_{n=-N/2}^{N/2-1} h[n] e^{-j \omega n} $$

remember, with a finite delay allowed for the filter (this is your 68 samples), you can model both your impulse response $h[n]$ and step response $s[n]$ as symmetric about $n=0$ and define the origin for $n$ there at the 68th sample.  doesn't matter if that makes the filter acausal.  don't sweat that.

also, keep in mind that the operation that turns the step response into the impulse response is a discrete-time differentiator, which means the information about DC is lost.  the DC component in your impulse response that is derived from your step response is the step size of the step response, $s[+68]-s[-68]$.  you can add to both any constant and that will not change the impulse response $h[n]$ at all.

oh, and expect 2×68 + 1 taps.


Here is the C++ code of the answer.

#include <algorithm>
#include <queue>
#include <numeric>
#include <vector>
#include <cmath>
const float N = 4;
float noise_th= 0.35;    //noise threshold
float calculate_SNR(queue<float> &signal_fifo)
{
    vector<float> noise_vect;
    vector<float> speech_vect;
    float SNR;
    float noise_std;
    float speech_std;
    auto fast_std = [](vector<float> input) //https://stackoverflow.com/questions/7616511/calculate-mean-and-standard-deviation-from-a-vector-of-samples-in-c-using-boos/12405793#12405793
    {
        float mean, acc;
        mean = accumulate(begin(input), end(input), 0.0)/input.size();
        acc = 0.0;
        for_each(begin(input), end(input), [&](const float d)
        {
            acc += (d - mean)*(d - mean);
        });
        return sqrt(acc/(input.size()-1));
    };

    while(!signal_fifo.empty()) //Seprate speech and noise samples(fifo has the ABS of input signal)
    {
        if(signal_fifo.front() < noise_th)
        {
            noise_vect.push_back(signal_fifo.front());
            signal_fifo.pop();
        }
        else
        {
            speech_vect.push_back(signal_fifo.front());
            signal_fifo.pop();
        }
    }

    noise_std = fast_std(noise_vect);
    speech_std = fast_std(speech_vect);
    SNR = 20*log10(speech_std/noise_std);
    noise_th = N*noise_std;     //Update threshold value
    return SNR;
}


Any meaning to the magnitude of a FFT when audio is recorded using a PC?

Actually, to say the either the sound values or the DFT magnitutes are unitless (or meaningless) is quite misleading.  The sound values actually represent changes in the sound pressure level at the microphone.  This pressure has physical units.  You can think of the microphone and A/D conversion as a "transducer" which converts the units into something else.  Ideally, this conversion is linear.  So the units in the sound file do have meaning, but most of the time it is not worthwhile to calibrate them to the actual physical units.

The DFT magnitudes are then in the same units, after you normalize by 1/N, where N is the sample count of your DFT.  The bit depth has to do with precision, not scale.

Your conclusion is still correct, but now better founded.  You can compare two signals on the resultant scale.

Ced


I took a look at your noise spectrum and it's pretty bad

You have strong lines at 20kHz and harmonics thereof. But also some spurs around  15 kHz, 18 kHz, 25 kHz and 44 kHz. The stuff at 20 you could potentially reduce with a comb filter, but your best bet (by far) is to clean up your measurement/recording setup.

There is no reason for this noise to be there, and you should try to identify the culprit(s) and eliminate them. Look at interconnects, other devices in the room especially lights or communication equipment, electric interference, power supplies, grounding issues etc.

This will be very hard to clean up, so having a good recording in the first place would be much preferable.


Impulse response convolution and normalization2

I think you're likely forgetting how the anti-derivative of $h(t)$ affects the gain of the convolution operation.  Recall that somewhere in your convolution integral, you'll be taking an integral of the form $\int We^{W\tau} d\tau$.  The Chain Rule requires the  $W$ in the exponent must appear as a $1/W$ factor after integrating.  This factor cancels the $W$ multiplier in $h(t)$ giving unity gain at DC.


GSM since it's first day has been digital (Fat32 is incorrect). That means you can't send analog audio. Luckily, analog audio had already been phased out in the 80s, in telecom.

In digital telephony, audio was digitally sampled at 8 kHz, 8 bits per sample, for a 64 kbit/s signal. If you has an ISDN phone, you got this quality at home, but it certainly was used inside and between the central phone switches.

GSM built on that 64 kbit/s idea, but used audio compression. The first codec used was GSM Full Rate, which achieved 5x compression for a 13 kbit/s net data stream.

GSM furthermore divides every frequency band into 8 slots, which means that phones can only transmit short 5 ms bursts and then have to wait for 7 other bursts from other phone users. Band spacing is a mere 200 kHz; a typical carrier gets 5Mhz blocks of spectrum which means 25 bands and 200 slots total.


Calculate the Inverse DTFT of the DTFT Derivative in Terms of $ x \left[ n \right] $

This is pretty straight forward using the definition of the Discrete Time Fourier Transform (DTFT).

The definition of the DTFT:

$$ X \left( {e}^{j \omega} \right) = \sum_{m = -\infty}^{\infty} x \left[ m \right] {e}^{-j \omega m} $$

Differentiating with respect to $\omega$:

$$\begin{align*}
\frac{d}{d \omega} X \left( {e}^{j \omega} \right) & =  \sum_{m = -\infty}^{\infty} \frac{d}{d \omega} x \left[ m \right] {e}^{-j \omega m} \\
& = \sum_{m = -\infty}^{\infty} \left( -j m \right) x \left[ m \right] {e}^{-j \omega m} \\
& = \frac{1}{j} \sum_{m = -\infty}^{\infty} m x \left[ m \right] {e}^{-j \omega m} \\
\end{align*}$$

Now, just plug it into the Inverse DTFT:

$$\begin{align*}
\frac{1}{2 \pi} \int_{0}^{2 \pi} \frac{1}{j} \sum_{m = -\infty}^{\infty} m x \left[ m \right] {e}^{-j \omega m} {e}^{j \omega n} d \omega & = \frac{1}{j} \sum_{m = -\infty}^{\infty} m x \left[ m \right] \frac{1}{2 \pi} \int_{0}^{2 \pi} {e}^{-j \omega m} {e}^{j \omega n} d \omega \\
& = \frac{1}{j} \sum_{m = -\infty}^{\infty} m x \left[ m \right] \delta \left[ n - m \right] \\
& = \frac{n}{j} x \left[ n \right]
\end{align*}$$

This implies:

$$ j \frac{d}{d \omega} X \left( {e}^{j \omega} \right) = j \frac{n}{j} x \left[ n \right] = n x \left[ n \right] $$


Establishing Phrase Frequency in Audio Files

This is not a simple problem.  Almost certainly any "home-brewed" system you create looking for "frequency patterns" is going to have horrible performance, unless it develops complexity comparable to full-fledged speech recognizer.  The reason for this is that recognizing speech, while it seems like a simple task to our ears, is in fact quite complex, and requires large amounts of context to have good accuracy.  What I'm trying to say is that if you just try to recognize the phonemes from the audio locally, you're going to get gibberish. You need to take into account the neighboring phonemes, the vocabulary of the language being spoken, and the large scale structure and content of the surrounding sentences in order to get anywhere near human level performance.

At the end of the day if you try to write something which "recognizes frequency patterns of speech" from scratch, you're going to be reinventing the wheel, and a pretty dang big wheel at that.

For this reason, the best thing for you to do would be to take an existing speech recognition system, and modify its output to suit your needs.  I would suggest using CMU Sphinx to decode the audio into an N-Best list, and then search this list for the text you are looking for.  CMU Sphinx has good recognition performance and is highly configurable.


Spatial multiplexing MIMO, antenna critical spacing, antenna isolation

What you ultimately want is that the channel coefficients are uncorrelated. MIMO (and diversity techniques in general) depend on the fact that it's unlikely that many different channels will all be bad at the same time. When the channels are correlated, the performance of your receiver will go down dramatically, because correlation means that when one channel is bad, the others are likely to be bad too.
Regarding spacing: if the antennas are too close together, then they see essentially the same channel. So, some spacing is needed; the more distance between antennas, the better. The $\lambda/2$ antenna spacing is a rule of thumb that may or may not work in practice.
Regarding isolation: It turns out that when antennas are close together, they become "mutually coupled". The voltage produced by each antenna depends not only on the received signal, but on the voltages produced by the other antennas. Mutual coupling is reduced with increased spacing, or it can be attenuated using "decoupling" circuits.
Both of these effects result in an increase in channel correlation, and a decrease in performance. Most textbooks mention only spacing, and not coupling, even though this effect has been known in the RF community for a long time. See for example "Mutual coupling effects on the capacity of multielement antenna systems", by T. Svantesson et al.
To answer your question: in an actual system implementation, both spacing and coupling must be taken into account. However, keep in mind that the ultimate objective is to have decorrelated channels; antenna spacing and decoupling are just two techniques to achieve that.


Difference between frame rate and sampling rate?

One method uses np.linspace() and the other one np.arrange() to create the time domain vector.

I suggest reading the documentation

https://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html
https://docs.scipy.org/doc/numpy/reference/generated/numpy.linspace.html

The first method has a small error, so the frequency will not be 100% correct.

Update

Common error with using linspace: Let's say you want three samples that are 1 apart. linspace(0,3,3) will give you [0 1.5 3], so the spacing is bigger than you want. The correct way is to adjust the upper limit by subtracting 1 interval (or 1/samplerate), so it should be linspace(0,2,3) or

linspace(0,5-1/framerate,framerate*5)


The problem is the assumed knowledge that the receiver needs to have. In your coding scheme you assume that the receiver knows that you transmit exactly $1000$ symbols. If the receiver didn't know that, there is no way to distinguish the following two cases:

$2$ tails at certain positions inside one block of $1000$ symbols
$2$ tails at the same positions as for case 1. but in two different blocks

So your coding scheme is incomplete and that's why you end up with an average bit rate less than what we would expect from the entropy of the source.
Also note that your scheme, even though incomplete, will exceed the minimum possible bit rate for larger block lengths. That limit of the block length beyond which the required rate becomes larger than the minimum possible rate can be computed as the smallest integer value of $N$ satisfying
$$p\lceil\log_2(N)\rceil>-p\log_2(p)-(1-p)\log_2(1-p)\tag{1}$$
where $p$ denotes the probability of a tail.
Apart from the above, you would also need to assign a codeword to the case that there is no tail inside a given block. Of course, for this likely case it would be wise to choose a short codeword.


After wasting forever on this problem, I found that it was just a question of window function. Doh!

Window functions that are not zero-ended produce these vertical lines. For example, Rectangular Window, Hamming Window, Gaussian Window (with low sigma) produce these lines, and Barlett Window, Hanning Window, Blackmann Window, Welch Window, don't.

I had tried other window functions but maybe during my tests, I made a coding mistake so I didn't find the fix. Moral of the story : do tests carefully!


The multiplication of the two unit step sequences $u[k]\cdot u[-n+k-1]$ is only non-zero if both sequences are non-zero. This means that the condition $k\ge 0$ as well as the condition $k\ge n+1$ must be satisfied. So you have two cases: for $n<=-1$ you have to evaluate the sum with the lower limit $k=0$, and for $n>-1$, you have to evaluate the sum with a lower limit $k=n+1$.


Music/Speech Classification: How are frame based features used in a feature vector?

Let us say you have a set of songs you want to classify by genre.

There are two approaches if the classification algorithm you use works on unidimensional vectors (eg: neural networks, SVM...):

Either you have one feature vector per song; and your feature vector contains statistics (mean, standard deviation, etc.) summarizing all the MFCCs computed at the frame level. In this case, your training set contains one vector per song; and during classification, the classification algorithm is run once per song.
Or you have one feature vector per frame. In this case, your training set is made of frames; and the classification algorithm is run once per frame. Voting is used to select the label for the whole songs. It is common to weight the votes to give more important to the middle of the signal (for example, for musical genre classification, the song might start with a slow/soft intro which might bias the classification).

A particularity of both solutions (the so-called "bag of frames" representation) is that they are invariant to permutations in the sequence of MFCC vectors - they do not take into account the temporal structure of the signal. Another approach is to use machine learning algorithms or models (such as HMMs) which can intrinsically deal with time-series - in which case your representation is the sequence of feature vectors.

The approach you describe in which MFCC vectors for each frame are concatenated to yield a gigantic feature vector is not practical: the large feature vectors would make the learning task computationally extensive, it would yield feature vectors of different dimensions for signals of different durations; and two signals identical modulo a delay would have very different feature vectors.


You don't. Only circles can be rotated by an arbitrary degree and still have the same shape. A 41x41 image is an axis-aligned square, and only 90 degree rotations keep the shape.

And so rotating over 60 degrees means that the width and height needed to fit the square have increased by sqrt(3/2)..


Why does an fftshift in the time domain give a different result compared to the frequency domain?

There are a couple issues with your code. First, let's talk about why you can do what you're trying to do.

What you're trying to take advantage of is the duality between a multiplication in the time domain and convolution in the frequency domain. The fftshift() function simply swaps the two halves of a vector (which is assumed to be the result of a call to fft()). This operation is equivalent to a circular rotation of half of the vector length in the frequency domain.

One can accomplish the same operation by multiplying by a complex exponential function in the time domain before the FFT, as you noted. To get an $M$-sample shift, you need to multiply by a complex exponential that looks like:

$$
s[n] = e^{\frac{j2\pi nM}{N} }
 $$

For the case where $M = \frac{N}{2}$, this simplifies to:

$$
s[n] = e^{j\pi n} = [1, -1, 1, -1, \ldots]
 $$

This illustrates problem #1 with your script: it multiplies by the inverse of the above: $-s[n]$. You need to invert the sense of the if statement in the for loop to ensure that the vector s has the appropriate value.

Problem #2 stems from another simple MATLAB mistake: all of the vectors in your simulation are length 50001. So, the time-domain rotation that you're doing by multiplying by the $s[n]$ illustrated above actually results in an effective frequency domain shift of 25000.5 bins (instead of the 25000-bin shift implemented by fftshift()). The mismatch in results is due to this extra half-bin shift.

You can fix this problem in two ways:

Change the time-domain multiplication to multiply by an actual complex exponential vector at the appropriate frequency (i.e. set $M=25000$ and $N=50001$ in the first equation). The resulting signal at the input to the FFT will be complex.
Change your simulation to use even-length vectors (i.e. chop off the last element of the t vector that you initialize at the beginning).

I tweaked your simulation to fix the generation of s and truncate the vector length to 50000. It is as follows:

%1khz
f=1000;
t=0:.0001:5;
t=t(1:length(t)-1);
x=sin(2*pi*f*t);
%Wave
subplot(4,1,1);
plot(x)
%FFT
subplot(4,1,2)
ffted=fft(x);
plot(real(ffted));
%FFT Shifted
subplot(4,1,3)
shifted=fftshift(ffted);
plot(real(shifted));
%FFT Shifted by preshifting
subplot(4,1,4)
s=zeros(1,length(x));%In languages like C++ this is computationally inexpesnive
for i=1:length(x)
    if (~bitget(i,1))%odd evens, works for real signal
        s(i)=-1;
    else
        s(i)=1;
    end
end
chirpshifted=fft(s.*x);
plot(real(chirpshifted));

I then compare the results of the two methods:

>> diff = abs(chirpshifted-shifted);
>> sum(diff)
ans =

   3.3406e-21

As you can see, the results agree as perfectly as you would expect from floating-point math.


Confused on application of Anti-Aliasing Filter

Math and Physics are not the same thing. They can be very close but infinity is hard to realize.  Going from analog to digital, a more accurate description is an average of the waveform over a small time window and the average is in turn quantized. Filters don’t have infinite stop bands but 80 dB is usually more than good enough.  There is always a small error, but one can bound it.


Basically, I see a plot of a 2D function of discretized scale and translation parameters. Instead of a smooth 2D surface, it looks like 1D plots of coefficients  at all scales $s_n$, put behind each other along each location $\tau_n$ on the translation axis. And each 1D plot is colored in a level-set fashion: the  "vertical" coloring below the 1D curve is related to the amplitude, potentially with color cycling: red for low amplitudes, then yellow, green, cyan, blue, magenta, red (for high amplitudes). This is apparently an instance of a waterfall plot:

  curves are staggered both across the screen and vertically, with
  'nearer' curves masking the ones behind

with amplitude coloring. So:

Coefficients (absolute value) gives you the height of the 1D curve (top) and the coloring below.
You don't have approximations and details with CWT. This is not DWT. "Only" low-scale to high-scale "detail" coefficients. There is not father wavelet or scaling function in that case.
Yes, without special settings, standard wavelet codes compute coefficients at each sample.

Alternatively, you can draw ribbon plots, with a Matlab code for the color (ribboncoloredZ.m):


According to MATLAB documentation on wavedec,

  The length of each filter is equal to $2N$. If n = length(s), the
  signals $F$ and $G$ are of length $n + 2N −1$ and the coefficients
  $cA_1$ and $cD_1$ are of length

  $$\text{floor}\left( \frac{n-1}{2} \right) + N$$

Here, $n=16$ is the length of your signal, and $N=2$ is the Daubechies number.

Putting all that together, your detail coefficients at the second level should be of length

$$\text{floor}\left( \frac{16-1}{2} \right) + 2 = 7+2 = 9.$$

At the second level, you coefficients should be of length

$$\text{floor}\left( \frac{9-1}{2} \right) + 2 = 4+2 = 6.$$

If you're wondering why this must be the case, imagine the filtering-decimation procedure. Scaling and wavelet function for $\text{db}m$ wavelets are of length $2m$. When you convolve length $n$ signal with length $2m$ signal, you get back a signal of length $l_0=n+2m-1$. If you take every second sample of this resulting signal, you get back something of length $l_1 = \frac{n+2m-1}{2} = \frac{l_0}{2}$. Of course, if $l_0$ is odd ($n$ is even), then we cannot split the signal exactly into two parts. With clever mathematics that I will not go into here, you can show that this last non-paired coefficient is redundant anyway (carries no information that we don't yet know), so you can just omit it. Therefore, we will always have the resulting decimated signal of length

$$l_1 = \text{floor}\left( \frac{l_0}{2} \right) = \text{floor}\left( \frac{n+2N-1}{2} \right) = \text{floor}\left( \frac{n-1}{2} \right)+N.$$


Pole zero plot, normalizing frequency response plot?

For the frequency axis, you're going to have to normalize by $\pi$ to plot digital frequencies in $[0,\pi]$.

For the magnitude axis, normalize by $\mathrm{max}\left(|H\left(e^{i\omega}\right)|\right)$.

I'm surprised you weren't asked to plot magnitude response in dB (a logarithmic axis).


Why duplicating the signal length alters magnitude in FFT spectrum?

This is due to Parsevals Theorem. The energy before and after DFT is equal (or depending on the normalization factor, at least proportional).

So, if you double the input energy by duplicating the input, you also get a different output energy. In your special case, you have a two-times repetition of the input, this corresponds to a discrete spectrum, where only every second bin is non-zero. Hence, the value on these bins needs to increase to fulfill the energy constraint.

Edit - More explanation on the discrete spectrum:
For the continuous-time Fourier Transform, every periodic time domain function of period $T$ has a discrete spectrum, where the spectrum consists only of discrete lines (Dirac impulses) with distance $1/T$.

A similar property holds for the DFT: if you input a periodic time function, you get out a discrete spectrum (i.e. it only every $n$th frequency bin is non-zero). Namely, if your time-domain input function has period of $N/L$, your spectrum is discrete, i.e. every $L$th bin is non-zero ($N$ is the DFT length).

In your case, your input function has period $N/2$, i.e. only every 2nd frequency bin is non-zero. Now, the double input energy (you have input 2 periods instead of one period) needs to be represented on these bins, leading to higher amplitudes compared to a single period.


Generation of noise floor file for estimating signal-to-noise ratio of a wave file

One method that is reasonable for white noise is to use a spectrogram. For each spectrum, compute a quantile value of the power bins and use that as the estimated noise floor. The quantile should be less than the maximum spectral occupancy of the speech. For example if the speech never fills more than 50% of the spectrum then the median (quantile 0.5) might be reasonable. If the speech fills more than half the spectrum then use a smaller quantile.


Cross-correlation and auto-correlation of an image using $\tt fft$

For 1D signals, the autocorrelation of a rectangular pulse will be a triangular pulse.
For your image, the autocorrelation will be a triangular pulse in both directions.
The image below shows:

The original image
The autocorrelation without zero padding
The autocorrelation with zero padding
The autocorrelation with zero padding, viewed as a mesh.

For your example image, the zero padding is not required because the support of the non-zero pixels is way less than the image size. For non-zero pixels covering the whole image, you'll need to zero pad when using the FFT.
Your image does not appear to show the gradual change. I suspect you will need to normalize the pixel values manually. It's not clear to me that 'CDataMapping','scaled' does the right thing (it may, I'm just not familiar with it).

what does the output value in each pixel mean for the auto and the cross-correlation?

As calculated, not much.  It would mean a little more if the mean was subtracted from the image before calculation. As it is, the  most you are seeing is the auto-correlation of the mean pixel value (128?) over a limited range.
If the mean were subtracted, then it would tell you about how each pixel influences the pixels around it: if the autocorrelation values are close to zero, then that means "not much". If the autocorrelation values are high, then it means "greatly".

R Code Below
# 32079

par(mfrow=c(2,2))
image <- array(0,c(20,20))
image[8:13,8:13] <- 1
image(image, col= grey(seq(0, 1, length = 256)))
title("Original")

fftimage <- fft(image)
congfft <- Conj(fftimage)
ans1 <- fftimage * congfft
ans2 <- Re(fft(ans1, inverse=TRUE))/20/20
image(ans2, col= grey(seq(0, 1, length = 256)))
title("Without zero padding")

imagePadded <- array(0,c(40,40))
imagePadded[11:30,11:30] <- image
fftimagePadded <- fft(imagePadded)
congfftPadded <- Conj(fftimagePadded)
ans1a <- fftimagePadded * congfftPadded
ans2a <- Re(fft(ans1a, inverse=TRUE))/40/40
image(ans2a/max(ans2a), col= grey(seq(0, 1, length = 256)))
title("WITH zero padding")

persp(ans2a, theta=30)
title('Mesh plot WITH zero padding')


Parameters of standard audio channel filters

Your question is not very specific, nor the answer will be. For example what do you mean by 'standard filter parameters'? I assume that you want to make given speech signals sound like ones being recorded with a specific microphone. Solutions would be

Hard way. You must obtain the impulse response of your microphone, for example by using Maximum Length Sequence or Sweep Sine method. In theory such measurement should be performed in anechoic environment (free space outside, well damped room, anechoic chamber), otherwise you will add even more reverberation. When you have your impulse response, simply convolve it with your signal and you will get very good approximation of how your microphone modifies the sound. Obviously speaker characteristics are also part of this impulse response, so you might want to use one with most linear amplitude response. In theory speaker response can be deconvolved to some extend by using reference microphone.
Easy way. Measure the frequency response of your microphone by using signals described above. You also measure the frequency response of the reference (maximally flat) microphone. Just subtract two amplitude characteristics and you get how your microphone distorts signals. Having that, use it as the input for the Parks–McClellan algorithm and design FIR filter that resembles this response in a best way.
Dummy way. Take the frequency characteristic provided by manufacturer in the manual, trust him, sample it, and use as the input for algorithm mentioned above. In MATLAB that would be firpm function.


Modeling audio effects, reverse engineering

Now is there a general method to analyze and determine the exact function of that effect

No, there can't be. Effects can (and will be) arbitrary, non-linear, memory-affected mappings.

You will have to make a model of what the effect does and then look for an estimator for the parameters of that model.

For example, if you radically restrict yourselft to LTI systems (e.g. linear filters), then you can represent the effect by its impulse response, which is easy to observe using  an impulse.

But, as said, the effects that are specific to expensive DAW will at least have memory and often quite a bit of nonlinear behaviour. Your model becomes very complex then, and there's no guarantee any model matches what your effect actually does!

Really, effects are arbitrarily complex computer programs, if you look at it this way. You want to figure out what a program does by observing its output. Can you tell me the source code of an encryption program when I give you the plaintext and key and ciphertext? How many sets of inputs and outputs would you need to be able to do that in general?

  cool audio effect that works well on guitar, sadly it exists only in the DAW environment

Don't underestimate the computational load of effects; many will simply be too hard for anything that's not as powerful as a PC. Modern effect boxes that can digitally emulate a wide range of analog distortions as well as apply purely digital effects use pretty beefy specialized digital signal processors, and/or FPGAs to even be able to compute these things.

Also, and that's not even solvable by throwing more compute power at the problem: A lot of post-processing algorithms simply aren't causal. That simply makes it impossible to apply live.

  would using Audio Neural Style Transfer be the best option

"Best" is never defined unless you give an objective measure for goodness. However, no:

Style transfer is style transfer. You want to estimate the properties of a system; that's a different problem.

But neural networks/ machine learning in general:

Yes! The central statement of neural networks is the universal approximation theorem, which says that given enough nonlinearly weighted linear projections, you can approximate any function.

I'm currently working with a student on training neural networks to learn how the nonlinear distortion of a high-power amplifier looks, so that he could use that model to train a predistorting network to counter these effects. It's working well, but the model we're using now is memoryless and relatively benign in nonlinear distortion. Learning an arbitrary guitar effect would, on our high-end GPUs, take months to learn (to the degree that satisfies our requirements).

Whenever you can, start with a known model rather then letting your machine learn from scratch. It shows that really, it's way quicker converging if you can offer a preprocessed version (that might not perfectly match the desired function). And you have a really hard problem: While my student has the luxury of not having to assume memory, you necessarily have to consider previous input. That makes the dimensionality explode.


Is there a difference measure that is insensitive to shifts?

Time scattering by design yields timeshift-invariant representations, with control over amount of desired invariance. An advantage is ability to compare within a frequency range of interest rather than in entirety.
It's the continuous wavelet transform, followed by modulus and lowpass filtering; nice lecture.

Distance with T=2048: 0.022988
Distance with T=1024: 0.172749
Distance with T=512:  0.415810

Code
Can only run from a dev branch, library's renovating for next release.
import numpy as np
import scipy.signal
from kymatio import Scattering1D
from kymatio.visuals import plot, imshow

def _l2(x):
    return np.sqrt(np.sum(np.abs(x)**2))

def l2(x0, x1):
    """Squared relative distance measure."""
    return _l2(x1 - x0) / _l2(x0)

#%% configure scattering & signal ############################################
N = 4096
Q = 16
J = 9
width = N//8
shift = N//8
# set amounts of invariance
T0, T1 = N//4, N//8
freq_fracs = (4, 8, 16, 32)

#%% make signal & shifted ####################################################
window = scipy.signal.tukey(width, alpha=0.5)
window = np.pad(window, (N - width) // 2)
t  = np.linspace(0, 1, N, endpoint=False)
x  = np.sum([np.cos(2*np.pi * N/ff * t) for ff in freq_fracs], axis=0
            ) * window
xs = np.roll(x, shift)

#%% visualize signals ########################################################
plot(x,  title="x, original",    show=1)
plot(xs, title="x, timeshifted", show=1)

#%% make scattering objects ##################################################
max_order = 1  # 2 yields better results but requires understanding; see docs
kw = dict(J=J, Q=Q, shape=N, average=1, out_type="array", pad_mode="zero",
          max_pad_factor=1, frontend='numpy', max_order=1)
ts0 = Scattering1D(T=T0, **kw)
ts1 = Scattering1D(T=T1, **kw)
#%% scatter
out0_x  = ts0.scattering(x)
out0_xs = ts0.scattering(xs)
out1_x  = ts1.scattering(x)
out1_xs = ts1.scattering(xs)

#%% compute & print distances ################################################
l2_0 = l2(out0_x, out0_xs)
l2_1 = l2(out1_x, out1_xs)
print("Distance with T={}: {:.6f}".format(T0, l2_0))
print("Distance with T={}: {:.6f}".format(T1, l2_1))

#%% visualize coeffs
imshow(out0_x,  title="scattering(x), original    | T=%s" % T0, abs=1)
imshow(out0_xs, title="scattering(x), timeshifted | T=%s" % T0, abs=1)
imshow(out1_x,  title="scattering(x), original    | T=%s" % T1, abs=1)
imshow(out1_xs, title="scattering(x), timeshifted | T=%s" % T1, abs=1)


Why the plots of spline and cubic interpolation are exactly same?

The best way to understand cubic splines (under whatever name or how specified) is to derive and code them yourself.
What you are looking for is a parametric equation with these properties:
$$
\begin{aligned}
f(0) &= y[n]    = y_0\\
f(1) &= y[n+1]  = y_1\\
f(2) &= y[n+2]  = y_2\\
f(3) &= y[n+2]  = y_3\\
\end{aligned}
$$
Your candidate function is a cubic polynomial, of course.
$$ f(t) = a t^3 + b t^2 + c t + d $$
Plug 'em in.  You get four equations with four unknowns ($a,b,c,d$)
Once you have solved the equation, use it to interpolate between the middle two points.
Extra credit:
Do the vector equivalent in N-space with the conditions given as two points in space and the velocity vectors at each point with one unit of time to get there.
Then code your solutions.  Print the results so you can see some precision.  And compare your answer to Matlab's.
Then post your results here.


Find out when an audio file is actually recorded?

Baring information content that contains verifiable historical landmarks, you can't date data .  Anything can be altered in a file, including the recorded content, so the data itself is no more reliable than a time stamp in a file header.  The normal assumption about data files is that they can be faithfully reproduced over and over without degradation, so there is no time information in the data itself.  You may be able to infer rough time estimates about the age of a recording based on certain qualities in the recording (assuming that the data has been unaltered).  For example a digital reproduction of an analog recoding that was originally mastered on a particular media has qualities that may be discernible upon analysis.  If you can identify qualities that indicate the type of equipment used to capture the audio and type of media the audio was captured on, you may be able to infer that the recording is from a certain time period.


Recommended C/C++ General Purpose DSP Libraries

Maybe a bit late, but since others might land here like me.

The following are good signal processing libraries/frameworks:
  * https://www.gnuradio.org/ - lots of classic signal processing + cascading options
  * https://root.cern.ch/ - advanced statistical signal processing
Both provide generic facade, fall-back implementations as well as performance optimised sub-libraries (ie. fftw, atlas, LAPACK, gsl, VOLK, ...).

Together with OpenCV, these probably cover 90++% of use-cases.


Highpass filter design that is better than the Moving Difference filter

[ 0.25      -0.6035         0.25       0.1035] seems to work.
Here is how I did that: Since the "breakeven point" is specified at $pi/2$, we know that $|K(0)|^2 = 0$, $|K(\pi/2)|^2 = 0.5$ , and $|K(\pi)|^2 = 1$ . The magnitude of the spectrum is thus given at 4 different points. DC and Nyquist are real and we have a conjugate complex pair at $\omega = pi/2$. We can use this to populate a 4 point frequency domain vector and perform an inverse FFT to get a 4 tab impulse response.
The only thing we don't know is the phase at $\omega = pi/2$. A highpass filter generally has a falling phase and the phase should be zero at Nyquist. Hence we should interpret the phase at DC to be 180 degrees and so it's reasonable to guess the phase at $\omega = pi/2$ half way in between and, $\phi = -pi/2$ seems like a reasonable first guess.
y = ifft([0 .707i 1 -.707i].');

You can play around with the phase and see if it can be made better (whatever that means), but given the amplitude constraints, any 4-tab solution must more or less look this way in the frequency domain.
EDIT
Turns out a phase of 75 degrees gives optimum performance in a sense that it both maximizes the gain at $\omega =3pi/4$ and minimizes the gain at $\omega =\pi/4$


Does Forward Error Correction (FEC) affect upsampling factor?

Is the redundant data introduced during forward error correction (FEC) counted as upsampling?

Yes.

BUT, under certain conditions.

The whole point of FEC is to be able to detect and correct an error instead of simply detecting it and asking for re-transmission.

To achieve this, we encode the source symbols in a predictable and reversible way, adding some redundant information that, upon reception, helps us infer the symbol that was transmitted even if what was received was incomplete or "damaged" (but always within reason).

This does not come for free. The code rate expresses the amount of redundancy inserted in the data stream. For example a convolutional 1/3 code adds two redundant code bits for every source bit. Consequently, we now have to send 3 bits down the same line, instead of 1. Therefore, our effective data rate is lowered.

One of the simplest ways to "encode" the source data is to repeat them a number of times. Instead of sending 1, send 11111 (a 1/5 code).

On the receiving side, you decode by majority decoding. If most of the symbols received were 1, it is inferred that what was intended to be transmitted was 1 (and conversely for 0).

Upsampling is composed of two steps: The actual upsampling followed by interpolation which is required to "make up" the values of those in-between samples whose value we don't know because it was never captured.

So:

if you were to take a data stream, upsample it and use nearest neighbor interpolation (effectively repeating each sample a number of times) and...
if you were then to send this augmented data stream down a communication channel, and...
if you were to now use a majority decoder on the other side of the channel to infer which one was the level that was actually sent...

then, yes, THAT could count as upsampling.

Otherwise, FEC and upsampling / interpolation are only somehow related as both operations can modify the data rate and content of the signal in a predictable and reversible way. But they do not overlap so much as to count as the same thing.

Hope this helps.


How to choose a windowing function when performing a Chirp $\mathcal Z$-transform?

When we add the window function to the definition of the DFT we get

$$ X(k) = \sum_{n=0}^{N-1} w[n]x[n] e^{-j \omega n} $$
with $\omega_k= \frac{2\pi k}{N}, k\in \mathbb{Z} , $ $w[n]$  the window function and $x[n]$ the data. The multiplication of the data $x[n]$ with the window function $w[n]$ in the time domain corresponds to a convolution in the frequency domain.

The CZT and the FFT are basically just different ways to compute this sum.  Each algorithm has its advantages and disadvantages. The FFT is fast, the CZT is more flexible. If you use both algorithms  to compute the same coefficients $X(k)$ you get the same results. So, in this case, the window choice considerations should be the same for CZT and FFT.

So far I can answer only have of the original question: If the CZT is used to compute the same frequency bins as delivered by the FFT, there is no difference with respect to window selection.

Note1: Are you aware of the issues regarding frequency resolution of different DFT algorithms, as discussed e.g. in this thread ? The frequency resolution topic is the same for "CZT vs. FFT" as for "zero padding of FFT" or for "Goertzel vs. FFT".

Note2: This was the basic message of my original post.

If  the CZT is used to compute the DFT on a finer frequency resolution, things are more complicated. In this thread the effect on the DFT is discussed. To make it short, the CZT interpolates the DFT. This interpolation makes window selection more involved.

If (as intended by the OP) the CZT is used to interpolate the DFT at a finer frequency resolution, the underlying maths remain the same, however the different frequency resolution must be kept in mind.

This is illustrated in a very simple example. The first figure shows the absolute values of a rectangular window. The FFT (without zero padding or similar things) computes the DFT right at frequencies where the window is zero, here marked in orange.

With finer frequency resolution this changes, as shown in the next figure.
 This may (as in the example) introduce additional leakage.


What is conceptually happening in this modulation function?

My question is what is significant of this two lines

The code receives a buffer of bytes but the modulation is based on a stream of bits. The outer loop loops over all bytes in the buffer. The
 for (int j = 0; j < 8; j++) {

This loops over all 8 bits in the current byte
int bit = (buf[i] >> j) & 0x01

This tests whether the j'th bit in the current byte is 1 on 0, where j = 0 produces the LSb and j = 7 produces the MSb.
If your current byte is 0b11001001, then the variable "bit" will be 1, 0, 0, 1, 0, 0, 1, 1 in each of the 8 loop iterations. Depending whether the current bit is 1 or 0, the modulated signal is generated in phase or out of phase.


Estimating confidence intervals for the mean value

To the extent you can assume your process is a white stationary ergodic process, the variance of the mean is the variance of the process divided by N where N is the number of samples.

Assuming that you can estimate the variance from your sample using

$$\sigma_x = \frac{1}{N-1}\sum_{i=1}^N(x_i-\mu_x)^2$$

NOTE: This is the unbiased estimator since the variance is being estimated from a sample of N elements using a mean that is estimated from the sample itself. (see Zwillinger, D. (Ed.). CRC Standard Mathematical Tables and Formulae. Boca Raton, FL: CRC Press, 1995.). Otherwise if the mean was known (which it isn't) there would be a $\frac{1}{N}$ in front of the summation above.

You should then be able to derive confidence intervals on the mean itself with the variance of the mean as $\frac{\sigma_x}{N}$:.

This is seen in the expression for the mean and variance of $Y$ where $Y$ is the average of $X$ over $N$ samples, as follows where $\mu_x$ is the mean of x and $\sigma_x$ is the variance of x:

\begin{align}
E(Y) &= E\left[\frac{1}{N}\left(X_1+X_2+ \ldots +X_N\right)\right] = \frac{1}{N}E\left[(X_1+X_2+ \ldots +X_N)\right]\\
\textrm{Var}(Y) &= \textrm{Var}\left[\frac{1}{N}\left(X_1+X_2+ \ldots +X_N\right)\right] = \frac{1}{N^2}\textrm{Var}\left(X_1+X_2+ \ldots+X_N\right)
\end{align}

\begin{align}
E\left[(X_1+X_2+ \ldots +X_N)\right] &= E[X_1]+E[X_2]+ \ldots E[X_N] = N \mu_x\\
\textrm{Var}(X_1+X_2+ \ldots +X_N) &= \textrm{Var}(X_1)+\textrm{Var}(X_2) + \ldots \textrm{Var}(X_N) = N\sigma_x^2
\end{align}

Therefore
\begin{align}
E(Y) &= \frac{1}{N}E[(X_1+X_2+ \ldots +X_N)] =\frac{1}{N}N\mu_x =  \mu_x\quad\text{and}\\
\textrm{Var}(Y) &= \frac{1}{N^2}\textrm{Var}(X_1+X_2+ \ldots +X_N) = \frac{1}{N^2}N\sigma^2 = \frac{\sigma^2}{N}
\end{align}

How to make use of the Autocorrelation Function

Note that this reduction by N is valid as long as the sequence is white. Once samples are correlated, there will be no further reduction in the variance of the estimate. You can therefore make use of your autocorrelation to determine the number of samples that will reduce the variance. I do not have an exact calculation for this, but to provide a rough order magnitude I would use the number of samples that results in the normalized autocorrelation dropping to below 0.5. For example, in your data the total number of samples is 30,000. If the autocorrelion immediately dropped below 0.5 after just one sample, then all samples are independent and your variance estimate would be the $\sigma_x/30,000$. However if the autocorrelation does not drop below 0.5 until 100 samples, then the variance estimate would be $\sigma_x/300$.


autocorrelation A to K conversion for LPC speech analysis?

After applying a Hamming window to your data I get very similar autocorrelation values:

R =
   1.6765e+00
   1.5828e+00
   1.3697e+00
   1.0688e+00
   7.4186e-01
   4.6660e-01
   2.5672e-01
   1.2121e-01
   4.9779e-02
   1.6453e-02
   3.1701e-03

You can check your LPC filter coefficients (not the reflection coefficients) by computing the matrix-vector product $\mathbf{R}\mathbf{a}$ where $\mathbf{a}$ is the vector of LPC coefficients and $\mathbf{R}$ is a Toeplitz matrix with the autocorrelation values:

$$\mathbf{R}=\begin{bmatrix}R(0)&R(1)&\ldots && R(p)\\
R(1)&R(0)&R(1)&\ldots&R(p-1)\\
\vdots& & \ddots & &\vdots\\
R(p) & R(p-1) & \ldots && R(0)\end{bmatrix}$$

This product must result in a vector with the first element being non-zero, and all the other elements must be zero (up to numerical accuracy).

Using the Octave function [lpc,v,k]=levinson(R(1:11)); with the autocorrelation function given above, I get the following LPC filter coefficients:

lpc =

   1.000000
  -1.370890
  -0.386146
   0.813691
   0.992372
  -0.949576
  -0.642884
   0.427407
   0.515102
  -0.310107
  -0.039544

And these are the corresponding reflection coefficients:

k =

  -0.9440538
   0.6826730
   0.2730461
  -0.1825513
  -0.4886926
   0.1882461
   0.3690086
  -0.0056205
  -0.3648883
  -0.0395440

So we get very similar values from which you can conclude that your code for computing the reflection coefficients from the autocorrelation function is correct.


How to select the length of a smoothing moving average filter?

The magnitude frequency response $|Y(\omega)|$ of moving average of length $N$ can be approximated by:

$$|Y(\omega)| \approx |\hat Y(\omega)| = \left|\frac{2\sin(\frac{N\omega}{2})}{N\omega}\right|$$

Figure 1. Horizontal axis: $N\omega$, vertical axis: $|\hat Y(\omega)|$.

The frequency $\omega$ is in radians, meaning that $\omega=2\pi$ corresponds to the sampling frequency. You should choose $N$ such that the noise is attenuated as much as possible while keeping a flat enough frequency response over the signal's bandwidth so that the signal's waveform is not distorted too much. What is flat enough depends on what you are going to do with the signal. A good first bet is to have the -3 dB cutoff frequency at the upper bandlimit of your signal so that no frequency of the signal is attenuated by more than 3 dB. The -3 dB cutoff frequency $\omega_c$ of moving average can be approximated by:

$$\omega_c \approx \hat\omega_c = \frac{2.783}{N}$$

These approximations were derived here.


How Is the Gaussian Kernel Related to the Euclidean Distance of the Neighborhoods in the Non Local Means Denoising Algorithm?

So we have 2 neighborhoods (Which are Matrices) $ v \left( {\mathcal{N}}_{i} \right) $ and $ v \left( {\mathcal{N}}_{j} \right) $. How could we calculate the distance between them?
One option would be the Euclidean Distance:
$$ {\left\| v \left( {\mathcal{N}}_{i} \right) - v \left( {\mathcal{N}}_{j} \right) \right\|}_{2}^{2} = \sum_{k} {\left( {v \left( {\mathcal{N}}_{i} \right)}_{k} - {v \left( {\mathcal{N}}_{j} \right)}_{k} \right)}^{2} $$
Yet this gives each pixel in the neighborhood window the same weight. The writes of the Non Local Means Denoising Algorithm thought it would be better to give the pixels near the center of the window higher weight.
They did it easily using by multiplying the difference window by a Gaussian Function and then summing. Something like:
$$ {\left\| v \left( {\mathcal{N}}_{i} \right) - v \left( {\mathcal{N}}_{j} \right) \right\|}_{2, a}^{2} = \sum_{k} {w}_{k} \left( a \right) {\left( {v \left( {\mathcal{N}}_{i} \right)}_{k} - {v \left( {\mathcal{N}}_{j} \right)}_{k} \right)}^{2} $$
Where $ {w}_{k} \left( a \right) $ is a weight based on a 2D Gaussian function with its center aligned to the center of the neighborhood window and its Standard Deviation is given by $ a $.
The Weight Function
In the classic Non Local Means implementation the Gaussian functions is used as weighing.
Assuming the $ v \left( \cdot \right) $ operator is the Vectorization operator (Column Wise) then for a $ 31 \times 31 $ window the following MATLAB code will define $ {w}_{k} \left( \sigma \right) $:
windowRadius    = 15;
kernelStd       = 3.5;

vX = -windowRadius:windowRadius;
vK = exp(-((vX .^ 2) / (2 * kernelStd * kernelStd)));
vK = vK ./ sum(vK);
mW = vK.' * vK;
vW = mW(:);

figure();
surface(vX, vX, mW);

The weight function is given by:

Indeed centered pixels will have higher weight.


Analog butterworth to digital - bilinear transform

look at MATLAB bilinear().  you are specifying fs=1 in your call to it.

if you're gonna do any digital filtering, sometime, somewhere you need to commit to a sampling frequency and you haven't yet.


Dummy DSP student request for IDTFT Calculation and Integration Process

This is using Euler's Identity for a sine which should get you past the point where you are stuck:

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


Is there any way to make a sampling rate reduction in the frequency domain?

Yes, you can do what you're looking for. Think about what happens in the frequency domain when you downsample a signal in the time domain: the spectrum "folds over" on top of itself and sums, the well-known aliasing effect. If your signal is already sufficiently bandlimited, then you can obtain the downsampled time-domain signal without computing the samples that you would subsequently throw out.

Given a signal $x[n], n = 0, 1, \ldots , N-1 $ that has DFT $X[k], k = 0, 1, \ldots , N-1$, a new signal $x_d[n]$ that is $x[n]$ decimated by a factor $D$ can be generated by:

$$
X_d[k'] = \sum_{d=0}^{D-1} X[k' + d\frac{N}{D}], k' = 0, 1, \ldots , \frac{N}{D} - 1
$$

$$
x_d[n] = IDFT(X_d[k'])
$$

The inverse DFT (and subsequently $x_d[n]$) is of length $\frac{N}{D}$. This method is sometimes referred to as "stack and add".


Relation between original points and 1st/2nd derivative points

Numerical signals are discrete and finite. So you are likely to hit issues either on end-points or discrete locations. I won't treat end-points here.

If you use only one operator, many signal processing folks care relatively little about a shift, as long as it is known, and constant. Here, the two-point derivative yields a $0.5$-point shift. You can consider that the derivative is valid at some mid-point. And you can interpolate it an integer locations if needed. Smoothing first can help limit a too-high fluctuations of differences.

A classical option is to choose an odd-length derivative operator. The three-point derivative writes:

$$ d_3[k] = (p[k+1]-p[k-1])/2$$

and gives you an information at the center of the schemes. Indeed, it coincides with the average of the left- and the right two-points skewed differences:  $ d_2[k] = p[k]-p[k-1]$ and $ d_2[k+1] = p[k+1]-p[k]$, since:

$$d_3[k] = (d_2[k] +d_2[k+1])/2$$

Many domains uses such schemes, and name them differently. To better find resources, some uses other terms, like generic ones, like "Finite difference coefficients", or specific ones, like "five-point stencil" or $n$-point gradient/Laplacian. One example is:

$$ d_5[k] = (p[k+2]- 8p[k+1] + 8p[k-1]-p[k-2])/12$$

The above thinking extends to higher order derivatives. One can further constraints methods using data properties, knowledge of noise, additional penalties, etc. To start with,  some literature on numerical derivatives:

CSERD Numerical Differentiation
Finite difference coefficients


pyrtlsdr: Why is my output different than universal radio hacker?

Your signal recording clearly shows that you have long streaks of 1.0 – that probably means you're clipping.

Your signal is thus broken. Make a new recording with less gain.


Does Scrambling increase speed of signal transmission overall? And is Scrambling a redundant process?

Scrambler is a randomizer with a purpose to encode transmitted data so that there will be no long sequences of only zeros or only ones in transmitted data. It is done to support proper timing recovery, when receiver waits for bit transitions to synchronize. Scrambling is usually done for channels with self-synchronized line codes, e.g. ISDN u-interface and Ethernet.
Generally, scrambler doesn't need to change the amount of transmitted bits, hence it doesn't change the speed of data transmission.
You can compare a scrambler to error-detection and correction schemes, which add redundant bits to support its functions. You can similarly compare it to 8b/10b encoding, when 2 bits are added to every 8 bits of data. As a scrambler doesn't add additional information to the data, it is not a redundant process in that sense.
Also, probably irrelevant in most cases, but it can be considered redundant, when used for a signal that inherently doesn't need scrambling. For example, UART uses short transactions and receiver is capable to get all bits only based on knowledge of a baud rate, so use of a scrambler would be redundant.


I'd say it's extremely difficult to obtain the sound pressure from a WAV file.

Consider what happens to the sound signal between the microphone and the WAV file. Pressure is applied to the mike. It is converted to an electrical signal, with power that depends on the properties of the mike transducer. Then it is attenuated on its way to the sound card, where it is amplified, filtered, sampled, quantized, encoded, and further massaged by the sound card driver and the OS sound system.

(Furthermore, if the signal is ultrasonic, the frequency response of your sound-card front-end may very well come into play too).

Then, in software, the samples are further processed to fit a particular WAV format (the details of which are amply documented in Wikipedia and other sources).

One possible way to proceed is to obtain experimental readings with your particular set-up: generate a set of known, controlled sound pressures on your microphone, and correlate them with the numbers you get in the WAV file.


why bandpass filtering (butterworth) centers the signal

High pass filtering means that low-frequencies, including DC, are removed. So if your signal has a DC offset (is not zero mean), then this will be removed by a highpass filter or a bandpass filter that does not include 0 in the passband.


How Is the Formula for the Wiener Deconvolution Derived?

The Wiener Filter can also be derived by another (Easier) way.

Let's assume the following model:

$$ y = h \ast x + n $$

Namely the data is a result of a linear combination (Convolution) of $ x $ with Additive Noise.

If we assume the noise model is Gaussian and our data is also formed by a Gaussian distribution then we should try to minimize (MAP Estimator):

$$ \hat{x} = \frac{1}{2 {\sigma}_{n}^{2}} \left\| h \ast x - y \right\|_{2}^{2} + \frac{1}{ 2 {\sigma}_{x}^{2} } \left\| x \right\|_{2}^{2} $$

Since the convolution is Linear Operation and assuming we in finite dimension then it can be rewritten as:

$$ \hat{x} = \frac{1}{2} \left\| H x - y \right\|_{2}^{2} + \frac{{\sigma}_{n}^{2}}{ 2 {\sigma}_{x}^{2} } \left\| x \right\|_{2}^{2} $$

This is just a Tikhonov Regulaized Least Squares with the solution given by:

$$ \hat{x} = \left( {H}^{T} H + \frac{{\sigma}_{n}^{2}}{{\sigma}_{x}^{2}} I \right)^{-1} {H}^{T} y $$

Now, if we use the discrete form of the circular convolution then $ H $ is a circulant matrix which means all can be solved in the Frequency Domain which results in the same equation as above.


Scale and Rotation invariant color based feature descriptors

There is absolutely no difference in the way these algorithms operate over grayscale and RGB images and colour SIFT algorithms simply take into account the additional information contained in colour.

These algorithms operate on typical images that are represented as matrices of (usually) integer numbers. Each number represents the intensity, or brightness, of light coming in to the imaging equipment from a particular direction. Colour information is (usually) represented with more than one information channels due to the way light works. These additional channels, also look like matrices. The raw form of a typical RGB optical photography image is three grayscale images packed together, OR, one image with 24 bit colour depth (8 bit for every colour channel).

There are however different ways by which you can visualise the RGB colour space. One of them is the Hue Saturation Value that returns one single number representing colour (Hue) another one representing the strength of colour (Saturation) and another one representing the brightness (Value).

The point in colour SIFT is to diminish the effect that colour has on tracking features but do it in the same way across all different colours (colour invariant). This can be done, either by pre-processing the input image (e.g. transforming it to a different colour space) so that the effect of colour is diminished and THEN sending it to a typical SIFT, or by incorporating colour invariant features to the standard keypoint selection process employed by SIFT (thus, cSIFT).

Besides SIFT, colour invariant metrics are also required in image indexing where they have been in use for longer, so you might get some additional information from "colour invariant metrics in content based image retrieval" (or this) type of searches.

Hope this helps.


Digital filter design of time series for specified frequency response function

This is the classical problem of filter design:

  You've got a frequency response that you want, how to implement it using a FIR?

I'm not going to lay out all the theory here, because it can easily be found by looking for Fourier approximation methods for FIR design or windowing methods, but the idea is easy enough:

Take your frequency response. Sample it at many points.
Take the inverse FT. You get something long, which, when convolved with a signal, would have the sampled frequency response (because of the convolution theory of the Fourier Transform).
Because the result is a very long filter and exhibits undesirable fringe behaviour, apply a window function to cut off some of the coefficients. The result is your filter taps.

There's iterative methods of increasing the match between desired and actual response given a limited amount of freedom by a limited number of taps; but this would exceed the scope of the question.


Why did a digital HPF increase the gain of high-frequency signal?

The following coefficients $a$ and $b$ implement a second order digital DC notch filter :

b = 0.9900   -1.9800    0.9900

a = 1.0000   -1.9800    0.9801

Whose frequency response is (y-axis in dB):

Given the raw-input signal $x[n]$ as a vector of samples in $x$ in matlab, you would apply the following line to filter it:

y = filter(b,a,x);

Note that the filter is not linear phase, hence you will have some phase distorions especially about DC frequencies... if you need a strictly linear phase DC blocking filter please indicate.

Given the raw data with the following plot:

The resulting filtered signal plot is the following:

It can be seen that initially there is a transient starting from a value of $0.428$ and takes some samples to stabilize to remove the DC. This transient is a direct consequence of long impulse response associated with the IIR filter.

Shorter impulse responses would provide faster response, but more steady state distortions.

Assuming that a shorter filter would do better, here is another set of DC notch coefficients with shorter impulse response:

b = 0.9100   -1.8200    0.9100
a = 1.0000   -1.8200    0.8281

In which case the output becomes:

Clearly an improvement in dynamic response, but keep in mind that the steadystate response now gives more distortions.

You can even further get a quicker response with the following set:

b = 0.7000   -1.4000    0.7000
a = 1.0000   -1.4000    0.4900

whose result is:

Much faster response, but there could be untolerable distortions on the samples, depending on the application requirements.


How to choose the number of codewords in bag of words model

It is a trade-off between recall and precision. If you have only one codeword in your codebook then everything will map to it and you will have high recall. If you have 20,000 codewords in your codebook then you will have high precision and low recall. You can use a validation set to try different codebook sizes and see which one works best for your application.

If you have a smaller size codebook I don't think that will help you get a sparse feature. It will just mean that more features map onto the same small set of codewords, increasing the recall.


MRI reconstruction using windowing based apodization

You say:

I have a 128 point one dimensional k-space samples...
The hanning window is the same size as the k-space vector (256)...

Make sure that you have the appropriate sizes in your algorithm.

Next, I convolved the k-space signal with the hanning window...

Windowing is applied in $k$-space by multiplication - you simply have to multiply your window and the original $k$-space line. In image space, this is a convolution, hence you could take the Fourier Transform of the window function and convolve the image with it.

Another view:
Apodization is applied to suppresses certain frequencies. For example, the highest frequencies are simply skipped by setting them to 0 due to the multiplication of the window function and the $k$-space data. If you were to convolve the two, you would "spill" the information of the higher frequencies all over $k$-space - which certainly is not what you want to achieve by apodization.