Explanation a flat fading channel in case MIMO OFDM

Frequency-flat fading essentially means that the channel's transfer function is (approximately) constant in your frequency band of interest, which means you can treat it as if it would not depend on frequency at all. This means two things:

Your MIMO transmission model is simply $\mathbf y(t) = \mathbf H \mathbf x(t) + \mathbf w(t)$ with a transmit signal $\mathbf x(t)$, a MIMO channel matrix $\mathbf H$, a received signal $\mathbf y(t)$ and your additive noise $\mathbf w(t)$. Note that $\mathbf x(t)$ needs to be narrow-band for this to make sense.
Your wouldn't use OFDM in this case. It just doesn't make sense: Every subcarrier would see exactly the same channel and you would gain nothing from the OFDM pre/post processing.

As for modelling the elements of your channel matrix, you can use statistical models that assume some random distribution that the elements are drawn from but be wary of the assumptions those are derived under. For instance, if you assume dense multipath from all directions, you end up with a Rayleigh distribution (resembling a case where you have rich scattering and no line of sight component). For cases with line of sight you might want to consider a Rician distribution. This is fine for not too high frequencies. For mmWave bands, this would be a bad model, since there the multipath is typically quite sparse and your channel matrix tends to be low rank.


The common definition of SNR is the power of the wanted signal divided by the noise power. Suppose you have obtained the wanted and the noise signal as arrays, calculation of the SNR in Matlab before noise reduction can be done like this:

snr_before = mean( signal .^ 2 ) / mean( noise .^ 2 );
snr_before_db = 10 * log10( snr_before ) % in dB

After noise reduction, the residual noise can be calculated as the difference of the wanted signal and the actual signal. Calculation of SNR is then straightforward:

residual_noise = signal - noise_reduced_signal;
snr_after = mean( signal .^ 2 ) / mean( residual_noise .^ 2 );
snr_after_db = 10 * log10( snr_after )


Estimating the complexity by the calculation of the number of real multiplications in a general polynomial function

The direct computation should be straightforward. Since I suspect the exercice has other purposes (more advanced algorithmic complexity), there additionally are a couple of tricks to reduce the complexity of power evaluations.

First, if you have to evaluate all powers $y^0,y^1,\,\cdots,y^Q$, you can store the power-of-two powers, and use the binary expansion of a power $q$. It can be called exponentiation by squaring. For instance, to compute $y^7$, you need one multiply for $y^2 = y\times y$, one for  $y^4=y^2\times y^2$ and two more for $y \times y^2 \times y^4$, so four multiplies instead of six by direct computation. And the $y^2$ and $y^4$  can be stored for the next powers.

Second, a classical polynomial evaluation is "Horner's scheme":

$$a_0 x^n + a_1 x^{n-1} + \cdots + a_n =  \left(\left(\left( a_0x+a_1\right) x \right)\ldots a_{n-1}\right)x+ a_n$$
which can be adapted to your case.


You can get pretty good results with an approximate linear chirp spread out over the period:

$$s(x) = \sum\limits_{n=1}^{N} \cos\left(2\pi \left(x n+ \frac{n^2}{2\times N}\right)\right), \quad N=10$$

$$m \approx 4.28$$

Group delay is the derivative of the phase with respect to frequency. For a linear (frequency as function of time) chirp, the phase should be a quadratic function of frequency. With discrete frequencies we don't get exactly that, but all frequencies seem to roughly find their slot in time if we sample the phase function. I wonder if a ping pong chirp would be possible here and if it would give a flatter envelope.


Video Compression - When would an average of the previous and next I or P frame be useful

Before elaborating anything, i would urge you to go through my answer here that simplifies some of the common confusions prevailing about motion vectors and estimations.

Your essential question is - when or how does encoder decides to put B frame or bidirectional motion vectors. I will break your questions into few as follows:

1. When does encoder decides to use B frames?

Almost always the number of B frames between I or P is always fixed. This is usually 2 or 3 pictures. The distance between I frame (called GoP size) may keep varying, but number of B frames will always remain fixed. From the point of view of MPEG, this is not really a restriction, but as one of your other question discusses, if you keep changing B frames, there will be serious problems of rhythms between the transmission and receptions.

Also, even before estimating whether a given frame will be a good choice to act as a B frame or a new I or P frame, you have to actually decide the forward reference and past reference for from where the prediction. Iterating over whether given frame will be perfect to work as P against frame already gone or should you choose some distance frame as P and current as B is too much overhead to get encoding done, and in the process you are fluctuating delay heavily.

So in essence, typically encoders decides number of B frames as 1,2 or maximum 3 and corresponding pictures are always fixed to be B frames.

2. So then what if i have a huge black/uncorrelated picture between last I and future P?

Ok - now that you have chosen to see that given frame is B, but almost nothing matches of that frame between either past or future; or sometimes only either side of the frame? This doesn't really put encoder in trouble though.

In MPEG, every frame is divided in terms of macro blocks. Each macroblcok has motion prediction from forward or reference blocks from respective past or future pictures.

I frame can ONLY have Intra blocks - which has no motion compensation, so it doesn't need any other frame ready for decoding itself.
P frame can have INTRA as well as FORWARD-PREDICTED macro-blocks. So even in case of P frame, when current block doesn't match with anyone in the past frame, P frame can simply use INTRA block.
B frame allows all type of macroblock including INTRA (as in I), FORWARD-PREDICTED mb (where motion vector depends only on past picture) , BACKWARD-PREDICTED mb (where motion vector depends only on future picture) and BIDIRECTIONAL macroblock which has two motion vectors from past and future and they are averaged out to form current block.

Hence, in case, the B picture or part of it happens to be very independent it is quite possible to put INTRA blocks for those MBs, and if it matches only either side of the picture encoder can use only FORWARD or BACKWARD prediction either way.

3. So does all MB's of a given frame have to be of same type?

Not at all, suppose a new small object just entered in a frame which is encoded as B frames, only part of the frame needs INTRA Macroblocks but majority of others can predict their motion from other blocks.

Understand, that all MPEG cares about is, if encoder can find a block which is almost as good as a current one, use to minimize the error, if not bad luck and encode this as a fresh new information. Even this ability to find itself varies from encoder to encoder.

4. So how does encoder make so many choices?

Usually, for the B pictures, encoder starts a usual business by searching the most appropriate blocks from past and followed by future picture.

a. if both frame produce some reasonable motion vector, this is most preferred to minimize the energy and book the block as BIDIRECTIONAL MB.

b. if either past frame or future frame doesn't produce motion vector - i.e. the minimum energy (SAD) is higher than a useful number, then correspondingly only 1 motion vector is encoded making MB FORWARD or BACKWARD prediction only.

c. if none of the reference frames have good enough motion vector - the MB is encoded as INTRA.

5. So why don't I simply encode everything as INTRA block and save from all the hassel?

You can absolutely do that. It will still be perfectly legal MPEG stream. Just that INTRA coding will take more symbols to encode and hence it will have to be quantized badly to fit under a given bit budget. Putting everything as INTRA has significant loss of efficiency because typically video always has huge temporal redundancy that INTRA blocks don't exploit; but they do give random access freedom.


ISI Nyquist criterion and spectral efficiency

I don't have Sklar's text but I guess he is referring to the baseband case. For bandpass signals the required bandwidth is doubled. So in the bandpass case you get a maximum spectral efficiency of 1 bit/s/Hz for BPSK, and $\log M$ bits/s/Hz for $M$-ary QAM.

Just as a side note, spectral efficiency is measured in bits/s/Hz (and not in bits/Hz as you seem to suggest).


How to fix frequency drift when producing an arbitrary length signal by repeatedly looping over a single-period waveform?

My suggestion is to use a Numerically Controlled Oscillator (NCO). A typical NCO effectively has one cycle of a sine wave in memory and can play back that sine wave as a continuous waveform with very high frequency resolution. But this needn't be a sine wave: you could have any waveform in memory and adjust the rate at which it plays back in the same fashion. That said, I believe these posts will be helpful where I provide more details on how an NCO works:

Numerically Controlled Oscillator (NCO) for phasor implementation?

Problem in understanding DDFS (Direct Digital Frequency Synthesizer.)

As a simplest explanation the approach is to count through the the waveform at a fixed rate (the sampling rate) with an increment that is proportional to the output frequency desired. The higher the precision of your counter (accumulator), the higher precision in frequency can be achieved. If you count with an increment if 1, you will get the lowest non-zero frequency out, and it will take $2^n$ counts to get through the entire waveform in memory once where n is the precision of the accumulator. Only the MSB’s of the accumulator are needed to address the memory based on how big the memory is. If you count by 2 that frequency will twice as fast, etc. So ultimately the frequency resolution is

$$f_\Delta=\frac{f_{clk}}{2^n}$$

Where $f_{clk}$ is the sample rate and n is the number of bits in the accumulator.

And the output frequency assuming the memory contains exactly one cycle of your waveform will be

$$f_{out}= Nf_\Delta$$

Where N is the count increment.

So for example with a 10 MHz clock (update rate) and a 32 bit accumulator which are both typical and practical values, you can generate output frequencies in the 0 to 5MHz first Nyquist zone with 2.33 mHz resolution!  (If combining with an D/A converter the practical upper limit would be approximately 40 MHz to allow for image filtering and an inverse Sinc filter would be considered to compensate for the Sinc reconstruction droop).

The links above go into full details of design parameters including signal quality (SNR and spur levels) and considerations on how many address bits should go to memory etc.


OFDM sensitivity to sampling time offset (with half of time interval between samples)

Yes, fine timing synchronization is still very important for OFDM. As you know, an FFT is typically used to demodulate an OFDM signal once it's been properly synchronized and stripped of cyclic prefix. To see what would happen if you aren't finely synchronized in time, recall the time-shifting property of the DFT:

$$
x[n-l] \Leftrightarrow X[k] e^{j\frac{2\pi kl}{N}}
$$

That is, an arbitrary timing offset $l$ (which need not be an integer for this formula to be valid) causes a linear phase shift in each output bin that is proportional to its distance from zero frequency.

So, your subcarriers near the center of the band will not be affected that much, but as you move outward, your demodulated signals will experience an increasing amount of phase rotation. That phase rotation is related to the unknown timing error $l$, and this value will likely vary slightly from symbol period to symbol period (assuming you don't have perfect symbol rate synchronization either). This all adds up to degraded performance.

The solution? You need fine symbol synchronization in OFDM, just like with any other digital modulation format.


How to implement Viola-Jones algorithm in Matlab?

Somebody did it for you:

Viola-Jones Face Detection for Matlab
Viola Jones Object Detection

and I guess you can find a few others. A nice description, in pseudo-code, can be found in An Analysis of the Viola-Jones Face Detection Algorithm, IPOL, 2014, which you can follow to code your own.


Try flattening your 2D array along the Y or frequency axis instead of along the X or time axis before doing the peak detect.


Bounds of higher order derivative of bandlimited signals

What simulations did you try? I think what you mention is right.
If $\mathcal{F}[f(t)] = F(\omega)$, according to the inverse Fourier transform, we have
$$
f(t) = \frac{1}{2\pi}\int_{-\infty}^\infty F(\omega)e^{j\omega t}d\omega
$$
$$
f'(t) = \frac{d}{dt}\Big(\frac{1}{2\pi}\int_{-\infty}^\infty F(\omega)e^{j\omega t}d\omega\Big)= \frac{1}{2\pi}\int_{-\infty}^\infty j\omega F(\omega)e^{j\omega t}d\omega
$$
So the Fourier transform of derivative should be
$$\mathcal{F}[f'(t)] = j\omega F(\omega)$$
and the Fourier transform of $n$-th order derivative is
$$
\mathcal{F}[f^{(n)}(t)] = (j\omega)^nF(\omega)
$$


Finding values of DTFT without explicitly computing

You just need to know the formulas for the DTFT and its inverse:

$$X(e^{j\omega})=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\tag{1}$$

and

$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\omega})e^{jn\omega}d\omega\tag{2}$$

From $(1)$ you see that your answer for $(a)$ is correct. Also the answer for $(d)$ follows immediately, if you realize that $e^{-jn\pi}=(-1)^n$.

For $(b)$ it's important to see that the coefficients are asymmetrical. Looking at $(1)$ you see that you can pair the terms for indices $\pm 1$, $\pm 2$, etc., which allows you to rewrite the sum $(1)$ as a sum of weighted sines time $j$ (imaginary unit), because $e^{jx}-e^{-jx}=2j\sin(x)$. So you'll have a purely imaginary expression. Computing the phase should then be easy.

$(c)$ is easily solved using formula $(2)$ (which value of $n$ do you need to plug into $(2)$?). And finally, for $(e)$ you're right that you need Parseval's theorem.


Digital receivers for the reception of non pulse shaped signals

Since you are talking about digital communication you are going to have a symbols.
When u send a symbol $s(t)$ you are going to send it with a symbol time $T_{sym}$
so we can say that any symbol that is sent equals to $s(t)\Pi(\dfrac{t}{2T_{sym}})$ where $\Pi$ is the rectangular function.
So basically any discrete time signal has a rectangular shaping pulse.

E.g.

In BPSK you will be sending one of the following symbols:
$$
    S_1(t)=1
$$
$$
    S_0(t)=-1
$$
Since I am sending the symbols for a known time $T$

I am practically multiplying them with a rectangular window with the length of $T$

so for the data series $a$ the signal $S(t)$ will be :
$$
   S(t) =  \sum_{k=0}^\infty S_{a_k}(t)\Pi(\dfrac{t-kT_{sym}}{ 2T_{sym} })
$$

And the matched filter in the receiver will be the average over $T_{sym}$ (or integration over $T_{sym}$) which is otherwise known as a rect.

Which is the same way you would do with other shaping pulses


In general, the data rate in bits per second is the symbol rate times the number of bits per symbol; independent of the modulation used.

The B in both BPSK and BFSK stands for binary, so there is 1 bit per symbol, or 2 possible symbols (2 possible phases or frequencies depending on the case).

So the data rate in both cases is the same.


You're thinking too complicated.

For a digital signal to be periodic, it simply has to repeat every $P\in\mathbb N$ samples.

Now, you have a continuous signal which is $\frac{2\pi}{\omega_0}$-periodic.

Obviously, if a period happens to be exactly "worth" a natural number $N\in\mathbb N$ sample times, the resulting sample from sampling it will repeat with a period of $N$.

More generally, if $M\in\mathbb N$ periods of the continuous signal are as long as $N$ sampling times, the digital signal will also be periodic.

So, as a general rule, one period of the input signal must be worth $\frac{N}{M}$ sampling times, or:

$$\begin{align*}
\frac{2\pi}{\omega_0} &\overset != \frac NM \frac {1}{f_\text{sampl}}&\quad N,M\in \mathbb N\\
&= R\frac {1}{f_\text{sampl}} &\quad R\in \mathbb Q\text,\\
\text{which means that}\\
\frac{f_\text{sample}}{\omega_0} &= \frac {R}{2\pi}\text;\\
\text{with}\\
f_{sample}&=\frac{\omega_s}{2\pi}\\
\text{it follows that}\\
\frac{\omega_s}{\omega_0} &= R\text,
\end{align*}$$

i.e. the relationship between signal frequency and sampling rate needs to be rational.

The fact that you have a zero-order interpolator afterwards doesn't change this consideration -- the interpolator's output will be periodic iff its input is.


How to Apply Least Squares Estimation for Sparse Coefficient Estimation?

I am not entirely sure what matlab's LASSO routine does so I started with Ordinary Least Squares (OLS) and worked backwards.  From an OLS perspective X1 as you have it won't work.  You've got a regressor that is as all ones, but your parameter inputs for you example data (r) doesn't contain an offset.  Essentially your model doesn't fit you're data well, which is why you get a weird answer.

Second, you have repeated regressors which is going to confuse things even more.  I've rewritten your code to shift your regressors to line up the three inputs that would result from MA using p = 3.

You can perform regular old OLS on the data X1 and y to get the exact input r.  If I use matlab's LASSO function I can get the proper components of r but they are off by a scalar.

    % Code for LASSO estimation technique for
    %MA system, L = 3 is the order,

    %Generate input
    x = -5:.1:5;
    xm = circshift(x,[0,1]);
    xm2 = circshift(x,[0,2]);

    r = 0.5*[1.0    1.0   0.0]'; % L elements of the channel coefficients

    %Data preparation into regressors
    X1 = [ x' xm' xm2'];

    y = X1*r;
    y1 = filter(r,1,x); % Generate the MA model, so I can plot and check.

    r_hat = X1'*X1\(X1'*y);

    [r_hat_lasso, FitInfo] = lasso(X1, y, 'alpha', 1, 'Lambda', 1, 'Standardize', 1);

Output:

    r_hat =
    0.5000
    0.5000
    0.0000

    r_hat_lasso =
    0.3235
    0.3231
         0


Fourier Transform of Morlet wavelet Function?

If you approximate the Fourier transform

$$X(f)=\mathcal F(x)(f)=\int_{-\infty}^\infty x(t)\,e^{-2\pi j\,ft}\,dt$$

by the discrete Fourier transformation for by sampling on the time segment $[-T,T]$ as

$$X(f_n)\approx \sum_{k=-N}^{N-1} x(k\tau)\,e^{-2\pi j\,f_nk\tau}\,\tau=s[n]\,\tau$$

with $T=N\tau$, $f_n=n/(N\tau)=n/N*f_s=n/T$, $n=-N,...,N-1$, $s$ the result of FFT on the sampled $x$ sequence after shift, ...

then you have to multiply the result of the FFT method by $\tau=T/N=1/f_s$.


Faster R-CNN and 4-step Alternating Training

Since no one answered this, I will post my attempt after further research.
Here is my understanding after having thought about it more. Lets say we are using VGG-16 as the backbone.
There are three fundamental components, the backbone, the RPN heads (objectness and RPN regression), and the ROI heads (classification and regression). Using the RPN heads, we can sample positive and negative ROIs.
Step 1) Train a region proposal network (backbone + RPN head) with the VGG-16 backbone pre-trained on imagenet with the head popped off. Some layers of the backbone will be trainable.
Step 2) Using the object proposals (300 or so per image) generated from the network obtained in step 1) (the backbone trained with the rpn), train a second network with another separate VGG-16 backbone pre-trained on imagenet, and a network head as in fast R-CNN (we will call ROI head). Then as in fast R-CNN we train for object classification of a region, as well as regression to learn the displacement from the proposal to a bounding box around the object (if the highest activation in the softmax does not correspond to a "background" class, i.e. no object).
Step 3) Use the backbone from step 2) for RPN training, which will be fixed (= not trainable), with the RPN head trained from step 1), and fine tune the RPN head.
Step 4) Finally keeping the backbone from step 2) fixed, and the RPN heads from step 3) fixed, fine tune the ROI Heads from step 2) to train bounding box regression and classification once more on the region proposals generated from
the backbone and RPN head.
Edit: Original answer had many errors so I came back to correct them, also original question doesn't really make sense now that I read it (I was new to this topic). Please just read this answer.


Correct comparison of FFT noise with FFT signal

It's interesting that there are no slam-dunk answers here - I think that's good, it means that there is no right answer and what works is going to hinge on the specifics of my implementation. Here is what I did:

I switched to measuring the power spectrum normalized for the Hanning window. This seems to give me a more stable center frequency and reduces the noise when there is no signal
I measured the noise on a per-bin basis and recorded the power spectrum and then averaged the output over a series of periodigrams. Since this reduces the gaussian noise I then multiplied by sqrt(N) where N was the number of periodigrams that I collected. I realise that this will amplify the non-random noise, but it also gives me a reasonable noise floor for the gaussian noise
I compared the noise power with the signal noise, converted to dB and allowed this to be configured by the user. 25dB works well in my case. I special-cased the situation where the noise was < 1 and made this the floor

This seems to work pretty well.


Why look at power spectral density for stochastic processes?

Because a stochastic process itself doesn't have a Fourier transform.

That's really all there is to it.

You can only transform signals (i.e. functions over a body isomorphic to $\mathbb R$, for example, functions of time). You can't transform a random variable whose individual realizations are such functions!


How to find phase shift and process I/Q channel

The way to do it is to correlate the received signal with the signal you sent out- likely a chirp or something.  You correlate it with the signal from antenna 1, which will give you a time (location of the cross correlation peak) and phase (phase of the complex value of the peak).  You do the same with the signal from antenna 2, getting a time and phase.  If all is right the time should be the same (or, depending on your time resolution, very, very close).  The phase difference is simply the difference of the phases of the cross correlation peaks from the two antennas.


In standard 1D signal processing, notch filters are filters that try to remove one specific frequency, and pass all others.

I say "try to" because this is not usually achievable, but it comes close.

Steve W. Smiths's DSP Guide has a nice plot (Figure 19.6b below) that shows the frequency response of such filters.

In image processing, it will be the same thing, but instead of a 1D filter it'll be a 2D filter probably comprised of two 1D notch filters --- one in each image direction.


It seems like you’ve by and large got it, but I did notice a couple things. You’ll probably want to include support for windowing and the corresponding correction factor as shown by APs article on FFT noise.  Rectangular windows aren’t great for RMS type measurements (at least in my opinion) for all the side lobe type explanations that are common with windows.  Second, I’m pretty sure you need to sum from 1 to N, not N/2.  You may have realized that because the input signal is real, that there is redundant data that doesn’t need to be calculated twice. That is true, with the exception of the 0 and N/2 (if N is even) bins.  Granted this probably wouldn’t make the output too wrong, but it would be off a bit.
Edit:  you know what, AP may have designed those coefficients so that windowing either isn’t necessary, or would do more harm than good.  You may still want to try it, since it’s easy enough to support, but it may not help.  Also, I don’t see how this would be more computationally efficient than doing the calculations directly in the time domain.


Find cut-off frequency of a filter given its transfer function

First you discretize you equation (note: I assumed that your $x(t-1)$ means a delay of one samples and not 1 second, but you can adapt easily if that's not the case) :

$$ x'[n] = x[n] - \alpha \times x[n-1] $$

Where
$$ x[n] = x(n\times T_s)$$

Then you use the $\mathcal{Z}$ transform :
$$X'[z] = X[z] - \alpha \times X[z]\times z^{-1}$$

And then you get the transfer function in the Z domain :
$$ \frac{X'[z]}{X[z]} = 1 - \alpha \times z^{-1}$$

This is a causal FIR filter, called a 1-zero filter with its zero at $\alpha$. Your cutoff frequency will depend on your sample frequency.

Illustration :

Code for the illustration :

figure,
hold on
for alpha=0.1:0.1:1
    h = freqz([1 0], [1 alpha]);
    f = 0:1/length(h):1-1/length(h);
    semilogy(f, 20*log(abs(h)));
    hold on
    grid on
end
legend('0.1', '0.2', '0.3', '0.4', '0.5', '0.6', '0.7', '0.8', '0.9', '1');
title('normalized frequency response')


Spatial Aliasing - Wrap Around F-K Spectra

Aliasing and frequency wrap is a consequence of violation of Nyquist-Shannon sampling theorem which states that a continuous signal must be discretely sampled at least twice the frequency of the highest frequency in the signal.Hence we need to briefly go into the mathematics of the same.

Let x(t) be a continuous signal, y(n) is discrete signal where y(n)= x(nT). Please note that x(t) will have a continuous time Fourier transform (CTFT) while y(n) will have a Discrete Time Fourier Transform (DTFT)

We can construct a mathematical model of sampling using the continuous Dirac delta function p(t) :
$ p(t)= \sum_{k=-\infty}^{k= \infty}\delta(t-KT) \forall t \in Real $

Let  w(t) =x(t)p(t).We can show that CTFT of continuous function w(t) is DTFT of y(n). Writing the frequency domain representation of w(t)

$W(\omega) = 1/2\pi \  X(\omega) * P(\omega) = \int_{-\infty}^{\infty} X(\Omega) P(\Omega -\omega)d\Omega \ \ \  $ , where * is convolution

now CTFT of p(t)  can be written as

$P(\omega) = 2\pi /T \sum_{-\infty}^{\infty} \delta(\omega -k\ 2\pi /T) $

thus $W(\omega) = 1/2\pi \int_{-\infty}^{\infty} X(\Omega)2\pi /T \sum_{-\infty}^{\infty} \delta(\omega -k\ 2\pi /T) d\Omega \\
=   1/T \sum_{-\infty}^{\infty} \int_{-\infty}^{\infty} X(\Omega)  \delta(\omega -k\ 2\pi /T) d\Omega \\
=1/T \sum_{-\infty}^{\infty}X(\omega -k\ 2\pi /T) \   \ \ \ \ \ \ \ \     ...using\  \ shifting\  property  $ ,

Thus we can say that $Y(\omega) = 1/T \sum_{-\infty}^{\infty} X((\omega -2\pi k)/T) $

So we can say that DTFT of y(n): $Y(\omega) $, is a shifted and repeated version of CTFT x(t) i.e. $X(\omega)$. The DTFT is the sum of the CTFT and its copies shifted by multiples of 2π/T.This is shown in the figure bleow. The frequency axis is normalized to −π/T < ω < π/T.

If X(ω) = 0 outside the range −π/T < ω < π/T, i.e. x(t) has no frequency greater than nyquist, then the copies will not overlap the range −π < ω < π and there is no problem as seen in figure above.

However if X has non-zero frequency components higher than π/T (fs/2 or nyquist). There will be overlap causing a wrap around in f-k domain. Refer below:

Notice that in the sampled signal, the frequencies in the vicinity of π are distorted by the overlapping of frequency components above and below π/T in the original signal causing wrap in F-k domain. Now in often spatial sampling is low as in case of Seismic or MRI due to constrains on number of detectors hence there is a aliasing in wavenumber. This spatial aliasing appears as wrap around in Frequency wavenumber Spectra.

Source of Images :R. G. Lyons: Understanding Digital Signal Processing (2nd Edition)

Reference : EECS, University of California Berkeley


Let me try to explain what happens to the PSD (power spectral density) of a discrete-time Gaussian white noise $x[n]$ of power $\sigma_x^2$ after it's either expanded by $L$ (insertion of L-1 zeros between consequitive samples of $x[n]$), or compressed by $M$ (aka decimation, selecting every M-th sample of $x[n]$).

First, let's put a convenient framework, by defining required quantities and their assumptions.

A discrete-time random process (RP) is constructed by a set of random variables (RV) and represented by $\{X[n,s)\}$ where integer $n$ denotes time-indexing for each RV and $s$ denotes outcome of a random experiment associated with those RVs or the RP. Customarily, an outcome $s$ selects a complete discrete-time waveform $x[n]$, which is an instance of the RP, from the ensamble set of it. The ensamble set includes all possible discrete-time waveforms. For notational simplicity, the sequence $x[n]$ is referred to as the random process, unless stated otherwise.

Probabilistic analysis of the RP $x[n]$ proceeds by defining various joint PDFs (probability density/distrbution functions) between random variables of the RP. Among those joint PDFs, (and moments associated) two of them become most useful and important for practical engineering applications. They're the mean and the Auto-Correlation Sequence (ACS) of the process, as defined below for a real $x[n]$ :

1- Mean of a RP is: $\mu_x[n] = \mathcal{E}\{x[n]\} $

2- ACS of a RP is: $\phi_{xx}[n,m] = \mathcal{E}\{x[n]x[n+m]\}$

As can be seen, in general the mean and ACS of a RP depend on time index $n$. The ACS depending also on the lag m. A very important simplification for the analysis of RP is made by defining what's known as a wide sense stationary (WSS) RP, whose mean and ACS are defined to be independent of time; i.e.,

1- Mean $\mu_x[n]  = \mu_x = \mathcal{E}\{x[n]\} $

2- ACS $\phi_{xx}[n,m] = \phi_{xx}[m] = \mathcal{E}\{x[n]x[n+m]\}$

Now, for given a WSS random process $x[n]$, its ACS $\phi_{xx}[m]$ holds a very important position. The power spectral density (PSD), $S_{xx}(e^{\omega})$, of a WSS RP is defined to be the discrete-time Fourier transform (DTFT) of its ACS:

$$ S_{xx}(e^{\omega}) = DTFT \{ \phi_{xx}[m] \} = \sum_{m=-\infty}^{\infty} \phi_{xx}[m] e^{-j\omega m} $$

At this point, it's clear that the PSD of a RP exists iff the RP is WSS and further that its ACS is stable so that its DTFT converges. I assume  @DilipSarwate may comment if a (major) ill-statement was presented up to here.

we are now ready to tackle your problem !

The compression by integer M (which you refer to as downsampling without filtering) operation is:

$$ x[n] \longrightarrow \boxed{ M \downarrow} \longrightarrow y[n] ; ~~~y=x[Mn] $$

Given an i.i.d (independent, identically distributed) zero mean,  WSS, (Gaussian) white noise RP $x[n]$ with power $\sigma_x^2$ it can be shown that its ACS and PSD are:

$$\phi_{xx}[m] = \sigma_{x}^2 \delta[m]$$
$$S_{xx}(e^{j\omega}) = \sigma_x^2 ~~~, ~~~\text{ for all } \omega $$

Then what is the PSD associated with the compressed output $y[n]$ ? To show this, first we must find its ACS and see that the RP $y[n]$ is indeed WSS.

$$\phi_{yy}[n,m] = \mathcal{E}\{y[n]y[n+m]\} = \mathcal{E}\{x[Mn]x[M(n+m)]\} = \mathcal{E}\{x[Mn]x[Mn+ Mm)]\} = \phi_{xx}[Mm] = \sigma_x^2 \delta[Mm] = \sigma_x^2\delta[m]$$

Hence we see (with a slightly bold step of treating $Mn$ the same as $n$ alone, which's justified by the WSS and whiteness of the input $x[n]$) that the ACS associated with the compressed signal $y[n]$ is identical to the ACS of the input $x[n]$. Hence their PSD's are also the same:

$$S_{yy}(e^{j\omega}) = \sigma_x^2 ~~~,~~~\text{ for all } \omega  $$

Next we apply the same for the expansion (you refer upsampling) process which can be shown to be as:
$$ x[n] \longrightarrow \boxed{ L \uparrow} \longrightarrow y[n] = \begin{cases} { x[n/L] ,~~~n=rL \\ 0 ,\text{ otherwise} }\end{cases} $$

Let's try to compute the ACS of the expanded output $y[n]$.
$$\phi_{yy}[n,m] = \mathcal{E}\{y[n]y[n+m]\} = \begin{cases} { \sigma_x^2 ,~~~ m=0 ~~~\text{ and } n=rL \\ 0 ,~~  m \neq 0~~~\text{ or }  n \neq rL }\end{cases} $$

When computing $\phi_{yy}[n,m]$, I assumed that the stuffed zeros were deterministic random variables with a pdf of $f_X(x)=\delta(x)$, which can only take a value of $0$ with a probability of one. They are not independent in themselves but independent from nonzero samples. Furthermore, since their value is always zero, when multiplied inside an expectation operator they always yield a zero result.

Now there is a problem! Eventhough the input $x[n]$ is WSS, the output is not. As its ACS $\phi_{yy}[n,m]$ depends not only on the lag $m$, but also on the exact time instance $n$. Hence its ACS is not of the form that's suitable for the DTFT expression. So at this point we can say that the PSD associated  with $y[n]$ does not exist in the form of $$S_{yy}(e^{j \omega}) = DTFT\{ \phi_{yy}[m] \} $$

May be a parametric Fourier transform or a Short-Time Fourier transform (STFT) or a 2D Fourier transform can be utilized but I'm not very sure about it.

Finally, we can state the following observation that since the expansion and compression operators are not LTI systems (they are time varying), then we are not surprised to see that their outputs are not WSS even if their inputs are so, which would not be possible with LTI systems.


How to calculate the discrete power of a discrete signal?

The instantaneous power is simply $p(k) = x^2(k)$. The signal energy is given by $W_x = \sum_{k=0}^{k=N} x^2(k) $.


How do you choose the optimal number of taps in a windowed image scaling filter?

Images quantized to some bit depth have error due to the quantization. Your filter also introduces error, compared to using an infinitely large ideal filter. Using larger and larger filters, eventually you will arrive at a knee in the total error vs filter size diagram where those two errors are equal. After that the quantization error dominates so it doesn't pay off much to improve the filter further.

A commonly used error metric is peak signal-to-noise ratio.


We did the same a couple of years back with what we called "echotimer", but not form FMCW, but for OFDM radar: basically, if you have a very large target that's very close and at 0 doppler, that's your crosstalk, and you can just subtract the distance of that.

Same applies here: with TX and RX hopefully sharing the LO, you're looking for a zero-velocity target that's close. You determine your range based on that, it represents your RX-TX time offset, and you just subtract.


would the linearity of a fourier property apply?

Of course, why not?

  How would you compute the power and the phase of the sum?

Computing each term in complex numbers, summing and converting to power/phase representation

  a step by step breakdown would greatly be appreciated

\begin{equation}
\begin{split}
C(t)&=\left(1.42 \exp^{-1.192t}- 12.44 \exp^{-1.192t} +11.02 \right) u(t) \\
 & = 1.42 \exp^{-1.192t} u(t) - 12.44 \exp^{-1.192t} u(t) +11.02 u(t)
\end{split}
\end{equation}

\begin{equation}
\mathcal{F}(C(t)) = 1.42 \mathcal{F}\left( \exp^{-1.192t} u(t)\right) - 12.44 \mathcal{F}\left(\exp^{-1.192t} u(t)\right) + 11.02\mathcal{F}\left( u(t)\right)
\end{equation}

Can you go from here? And why aren't you summing your exponential terms?


Does taking the Hilbert transform extract the envelope of audio signals?

That really depends on your definition of "enevelope" and what you need it for.

The Hilbert transform calculates the "analytic" signal, i.e. it calculates a matching imaginary part to a real signal by shifting the phase by 90 degrees in the frequency domain. It's reputation of calculating the "envelope" comes mainly from communication technology. It works very well with a narrow band signals, like amplitude or frequency or phase modulated signals. It's based on the simple fact that the analytic signal of a sine wave is a complex exponentional, the magnitude of which is 1. So the hilbert transform (or to be precise calculating the magnitude of the analytic signal) will eliminate the carrier quite well and leaves you with the modulation signal, or "envelope".

For broad band music, this is mostly not the case and the analytic signal isn't particularly helpful. A cymbal crash is basically a noise burst. The hilbert transform of this is just another noise burst and so is the magnitude of the analytic signal. You can't eliminate a carrier, if there isn't one to start with.

IF you are looking for something that roughly gives you "perceived loudness as a function of time", there many different ways to do this from simple lossy peak detectors all the way to sophisticated perceptual algorithms including cochlea models, binaural & temporal masking, nerve excitation models etc. It really depends on your application.


Is my algorithm about blob tracking for vehicle following correct?

You might want to consider tracking using MeanShift & Cam Shift:

http://opencv-python-tutroals.readthedocs.org/en/latest/py_tutorials/py_video/py_meanshift/py_meanshift.html#meanshift


What is the difference between photon shot noise and salt and pepper noise?

photon shot noise follows a Poisson Distribution, which is essentially a shifted Normal Distribution. It is an electrical noise usually caused by charge surges. Salt and Pepper Noise is more due to digital conversion, and is usually modelled using a Fat-tailed distribution (which is essentially a shifted Normal Distribution). What to note is for any unaltered image, the order of magnitude of each type of noise, with photon shot noise usually much lower than salt and pepper noise (unless specifically created)

True both noise can be treated using median filter, and often people use salt and pepper noise to represent shot noise.


Is there a noise cancellation technique for over-lapping spectra?

First of all the Wiener filter does not remove the noise but reduce it for WSS signals.

It does this based on the relationship between the power spectral density of the clean signal $x[n]$ and imposed noise $v[n]$, both assumed as WSS. The frequency response of a noncausal IIR Wiener filter is given as

$$
H(\omega) = \frac{ P_x(\omega) }{ P_x(\omega) + P_v(\omega) }
$$

where $P_x$ and $P_v$ are PSDs of clean signal and uncorelated noise imposed on it. As can be seen, the filter spectrum involves all the bands including where noise and signal spectrums are overlapping too, and that's one of the benefits of using such an advanced tool.

Note that this is not a non-destructive noise reduction and it wont necessarily preserve the waveform of the clean signal, but its PSD is improved; made closer to clean signal's PSD.


If the GPS signal is centered (carrier) at 4.309 MHz, and this is sampled at
5.714 MHz, then there would be an alias copy at 5.714-4.309 = 1.405 MHz.  To
acquire this signal, you would first "down-convert" the signal to baseband, or
0 frequency by multiplying by $exp(-j\omega t)$ where $\omega$ is
$2\pi(1.405e6)$ Hz. The reason to use an exponential and complex downconverter
such as this process instead of just multiplying with a sine, is that it will
facilitate fine adjustments to carrier frequency and phase as you do your
acquisition (which you can do manually to get a feel for how a carrier tracking
system would work).  After you multiply the signal as such to move it to 0
frequency, then correlate with the code for your specific SV of interest.

When correlating to a known code, the magnitude of the correlation will be
proportional to the residual frequency offset of your carrier as a Sinc
function with the first nulls located at $F= \frac1T$ where $T$ is the length
in time of your correlation. (For example, if you choose to correlate over a
full PRN which is 1 ms, the null to null width will be ± 1kHz, meaning you
should be able to see an SV even if it is still ±500 Hz off. For this reason if
I was doing a manual search with a full length PRN, I would step my frequency
in steps of 500 Hz and do a correlation in each "frequency bin" until I stepped
over the full possible Doppler range for my system. For your case to give a
specific example, I would do a first correlation after multiplying by -1.405
MHz ($e^{-j2\pi1.405e6t}$), and know that I would see and SV if it was within
500 Hz (the null extends to 1kHz offset, but beyond 500 Hz and the correlated
signal would start to compete with noise). If I didn't see anything for an SV
that I knew was present, I would then step my frequency up 1 kHz (to -1.406 MHz)
and try there, and I would continue this over all possible Doppler offsets
(typically up to ±5 kHz). If I knew my signal was strong, I could correlate over
1/2 a sequence to double the "window" in which I would see a signal in Doppler
offset. But this would reduce the processing gain by 3 dB but useful for an
initial acquisition search.

Note, if needed see my code generator available here: https://www.mathworks.com/matlabcentral/fileexchange/14670-gps-c-a-code-generator  which includes the ability to resample the code to any arbitrary frequency (in your case 5.7144 MHz).

And if you really want to get fancy (at the expense of resources associated with a 2D FFT), see my implementation for "one-shot" code delay and carrier phase here: https://www.mathworks.com/matlabcentral/fileexchange/26035-joint-frequency-and-delay-correlation   This shows how to determine code delay and frequency offset from one block of data in one matrix operation, using the cross-correlation theorem to relate the fft to a circular correlation;

$\text{xcorr}= \text{ifft}\left(\text{fft}\left(a\right)\text{fft}\left(b\right)^*\right)$

combined with the property that a circularly rotated FFT of a time domain signal is the same as the time domain signal in a different Doppler bin. (So once one FFT is computed for the CDMA PRN, you have the FFT's for the PRN at all other Doppler bins just by rotating that one FFT!).   The code provides a similar plot as shown below which is the correlation for an SV vs code and frequency offset.


How to approximate the transfer function of compensator using root locus?

K is the gain parameter for which you will be at any particular point on the root locus. As you increase K, the closed loop poles will start from the open loop poles and move toward the open-loop zeros as depicted in the root locus showing all possible locations of closed loop poles versus gain K.
This is because the general form for the closed loop and open loop gain in a general loop equation is related as:
$$G_{CL}(s) = \frac{G_F(s)}{1+K G_{OL}(s)}$$

Where $G_{CL}(s)$ is the closed loop gain, $G_{F}(s)$ is the forward gain from input to output and $G_{OL}(s)$ is the open loop gain with the real gain $K$ pulled out as a independent variable that we will sweep when creating the root locus plot.
The denominator of the closed loop gain is called the "characteristic equation" and given properness of realizable causal systems (properness means the number of zeros is always less than or equal to the number of poles, not including poles and zeros at infinity), we can determine all loop characteristics (stability, dynamics, etc) from the characteristic equation alone. To really see what is occurring in the root locus plot, consider this equation expressed in terms of poles and zeros of the open loop gain $G_{OL}(s)$ as follows:
$$1+KG_{OL}(s) = 1 + K \frac{G_z(s)}{G_p(s)}$$
Where $G_z(s) = (s-z_1)(s-z_2)...$ is the polynomial of open loop zeros, and $G_p(s) = (s-p_1)(s-p_2)...$ is the polynomial of open loop poles. This is rearranged to be:
$$= \frac{G_p(s)+KG_z(s)}{G_p(s)}$$
Since this equation is in the denominator of the closed loop gain $G_{CL}(s)$ then the roots of the numerator (the zeros of this equation) will be the poles in the closed loop gain!
Knowing that, consider when $K$ is varied from very large to very small: When $K$ is very small then $G_p(s)$ dominates the equation $G_p(s)+KG_z(s)$ and therefore dominates the poles of the closed loop gain. Specifically in this case of very small $K$, the poles of the open loop gain will be the poles of the closed loop gain.
Similarly, when $K$ is very large then $G_z(s)$ dominates the equation $G_p(s)+KG_z(s)$ and therefore dominates the poles of the closed loop gain. Specifically in this case of very large $K$, the zeros of the open loop gain will be the poles of the closed loop gain.
This is why the poles start at the open loop poles and move toward the open-loop zeros in the root locus.


Objective Noise Classification for speech enhancement

Environmental Sound Classification might be suitable for what you are asking for. Datasets include Urbansound8k, ESC-50 and Audio set. Pretrained classifier like VGGish are available, which cover the 527 classes in Audioset.


I am looking for an analytic description of a continuous-time Butteworth High-pass filter in the time domain (=impulse response)

$\frac{s^n + b_{n-1}s^{n-1}+\cdots+b_0}{s^n + a_{n-1}s^{n-1}+\cdots+a_0} = 1 + \frac{(b_{n-1}-a_{n-1})s^{n-1}+\cdots+b_0 - a_0}{s^n + a_{n-1}s^{n-1}+\cdots+a_0}$, and the inverse Laplace transform of 1 is $\delta(t)$.
Is that enough to get on with?


Distinguish harmonic sounds from non-harmonic

A standard difference between harmonic and non-harmonic sounds is that harmonic sounds tend to be periodic (in the time domain) in nature whereas non-harmonic sounds don't show any periodic properties since they're completely random.

One initial method you could try is simply performing autocorrelation in the time domain to see whether that portion of your signal is periodic. This could help you differentiate between the different sounds. There are definitely more complicated ways to test this. You could go to the extent of training a machine learning model based on certain features in the audio signal. But maybe you could just give the autocorrelation a try and see how well it performs for your application.


How to deal with "weird" phase plots in bode diagram when designing a controller

I'm a bit puzzled with your Bode plot, here's what I get with your transfer function and Matlab 2019b. My phase starts at almost -180 degrees while yours start at 0 degree. Is it possible that you've made a mistake somewhere in your analysis?

Your transfer function is a bit complex, you could simplify it by removing the fast poles and zeroes. For example, you have 2 fast poles at -1714 and -516 and a fast zero at -1088. You can remove the 2 fast poles and the fast zero to simplify your transfer function. This will make your analysis easier.
Secondly, I think the rule-of-thumb formulae you use for the cross-over frequency is only valid when your open loop system has a low-pass characteristic. Your system has a band-pass characteristic.


Discrepancy when calculating LTI system output using inverse z-Transform

As Matt L pointed out, your initial condition is wrong. y[0] is NOT 0, it's 4. This is obvious from the difference equation. Once you use y[0] = 4, both results come out to be the same.


Average identical sweeps or lowpass filter sweeps first?

If your are using a linear time invariant lowpass filter (such as a Finite Impulse Response (FIR) or an Infinite Impulse Response (IIR) filter), and have negligible rounding errors in each operations then doing the lowpass filtering before or after the averaging will yield the same results.
Note that from an implementation perspective, doing the averaging first would likely require fewer computational resources.


At it's most fundamental, the DFT is about fitting a set of basis functions to a given set of sampled data.  The basis functions are all sinusoidal functions, expressed as the complex exponential with a purely imaginary exponent.  Using the most common scaling convention each basis function, without its scaling coefficient, is:
$$ g_k[n] \triangleq \tfrac1N e^{j \omega_k n} $$
For the DFT:
$$ \omega_k = \frac{2 \pi k}{N} $$
and the basis functions add up as
$$\begin{align}
 x[n] &= \sum\limits_{k=0}^{N-1} X[k] \ g_k[n]\\
\\
  &= \sum\limits_{k=0}^{N-1} X[k] \ \tfrac1N e^{j \omega_k n} \\
\\
  &=  \tfrac1N \sum\limits_{k=0}^{N-1} X[k] \ e^{j \frac{2 \pi k}{N} n} \\
\end{align}$$
It's real easy to solve for the coefficients:
$$\begin{align}
 X[k] &= \sum\limits_{n=0}^{N-1} x[n] \ e^{ -j \frac{2 \pi k}{N} n} \\
\end{align}$$


How to detect Frequency Sweep Direction when Channel is unknown

See below for a python matched filter implementation. I generate an up-chirp and down-chirp. Then, subsequently I filter this signal by two matched filters, each filter is matched to one signal shape. As you see, there are two very strong peaks in the matched filtering. The reason why there are two peaks instead of (what you expect) slowly growing and falling triangles is that the chirps create an orthogonal basis, i.e. two time-shifted chirps are orthogonal, i.e. the MF output is zero. (The chirps are only orthogonal for infinite-time, which we dont have here, so there is still some small non-zero MF output).

import scipy

Fs = 10000
t = np.arange(0, 2, 1/Fs)

f = 500*(2-t)

signal = np.sin(2*np.pi*f*t)  # generate an down- and upchirp

f_ax, t_ax, Sxx = scipy.signal.spectrogram(signal, Fs)
plt.figure(figsize=(8,8))

plt.subplot(211)
plt.pcolormesh(t_ax, f_ax, Sxx)

plt.subplot(212)
plt.plot(t, np.real(signal))

chirpDown = signal[:int(Fs)]  # extract both signal shapes from the signal
chirpUp = signal[int(Fs):]

rx = signal + 0*np.random.randn(*signal.shape)

# perform the matched filtering (this means convolving with the time-reversed sequence)
t_c = np.arange(0, 3, 1/Fs)[:-1]
plt.plot(t_c, np.convolve(chirpDown[::-1], rx))
plt.plot(t_c, np.convolve(chirpUp[::-1], rx));

So, in order to detect those peaks, you can select the maximum of the MF output for each symbol duration (in case you are not perfectly synchronized) and compare which MF had higher output. Then, you reset the maximum value and search the next symbol period. Alternatively, you can convolve the filter output with some triangle or Gaussian, to get wider peaks.  For example, you can do the following to get an output similar to your posted figure:

triang = 1-abs(np.arange(-1,1,1/Fs))
# perform the matched filtering
t_c = np.arange(0, 3, 1/Fs)[:-1]
mf1 = np.convolve(chirpDown[::-1], rx)
mf2 = np.convolve(chirpUp[::-1], rx)
plt.plot(t_c, mf1)
plt.plot(t_c, mf2)

plt.figure()
plt.plot(np.convolve(mf1-mf2, triang))


2D Digital Signal Processing - Applying High Pass on and Image

First things first: If you use fftshift, you have to undo it before doing the ifft. Try

af=fftshift(fft2(k));
aflp=af.*lp;
aflpi=ifft2(ifftshift(aflp))

The result should improve, especially this checkerboard-overlay should go away.


How to proceess Image with curved corners and edges to give sharp Corners?

here's how houghLinesP would look like. Not a solution yet, but maybe you can work with that:

    int main()
    {

        cv::Mat image = cv::imread("SharpenCorners.jpg");

        cv::imshow("image", image);

        cv::Mat gray;
        cv::cvtColor(image,gray,CV_BGR2GRAY);

        cv::imshow("gray", gray);

        cv::Mat output; image.copyTo(output);

        // gradient:
        cv::Mat gx, gy;
        cv::Sobel(gray, gx, CV_32F, 1, 0, 3 );
        cv::Sobel(gray, gy, CV_32F, 0, 1, 3 );
        cv::Mat gm;
        cv::magnitude(gx,gy,gm);

        //cv::Mat g_thres = gray < 100;
        cv::Mat g_thres = gm > 100;

        //int lineThres = 5;    // 1st image
        //int lineThres = 10;// 2nd image
        //int lineThres = 25;// 3rd image
        int lineThres = 50;// 4th image

        std::vector<cv::Vec4i> lines;
        cv::HoughLinesP( g_thres, lines, 1, CV_PI/(4*180.0), lineThres, 10, 10 );

        for( size_t i = 0; i < lines.size(); i++ )
        {
            cv::line( output, cv::Point(lines[i][0], lines[i][1]),
                    cv::Point(lines[i][2], lines[i][3]), cv::Scalar(155,255,155), 2, 8 );
        }

        cv::imwrite("houghpSharpen4.png", output);
        cv::imshow("g thres", g_thres);
        cv::namedWindow("output"); cv::imshow("output", output);
        cv::waitKey(-1);

        return 0;
    }

with these results for given thresholds:

1st:

2nd:

3rd:

4th:

Additional idea:

Giving your own comment, that you might try to approximate each contour with rectangles, you could start with houghLinesP results (with a pretty high threshold to get at least all "long" line elements) and try to create the rectangles based on those lines


How to calculate the delay generated by a 2nd order all-pass filter

The phase of a second-order all-pass system with poles at $z=re^{j\theta}$ and $z=re^{-j\theta}$ is $[1]$

$$\phi(\omega)=-2\omega-2\arctan\left[\frac{r\sin(\omega-\theta)}{1-r\cos(\omega-\theta)}\right]-2\arctan\left[\frac{r\sin(\omega+\theta)}{1-r\cos(\omega+\theta)}\right]\tag{1}$$

The poles are the roots of the denominator polynomial $z^2+a_1z+a_2$. For the given values of $a_1$ and $a_2$ you get $r=0.93264$ and $\theta=0.040271$.

The delay of a sinusoid with frequency $\omega_0$ when passing through the all-pass filter is given by the phase delay

$$\tau_p=-\frac{\phi(\omega_0)}{\omega_0}\tag{2}$$

Evaluating $(2)$ using $(1)$ with $\omega_0=2\pi\frac{262.5}{44100}$ gives a phase delay of $42.7$ samples, which is pretty close to the value you've determined experimentally.

$[1]$ Oppenheim, Schafer, Buck, Discrete-Time Signal Processing, 2nd ed., p. 275.


How is Weighted Thresholded Histogram Equalization different from Gamma Correction Image enhancement

This Method can very well pass for being Gamma Correction Technique as We are Increasing the Weight of the Less probable Intensities in the Image  by the Exponential Factor r So that the Less Clear Details may be Visible upon the increase in their Weight thus Providing a Enhanced Image , On the Similiar Note the Gamma Correction Function Does this by using Power Law Transformation.


Long term spectral flatness measure voice activity detection

Or you can try and implement the adaptive threshold method:
This is taken from the EURASIP Paper which i have implemented.

UnlikeM and R, a fixed threshold would lose its efficiency
when facing varying acoustic environments. Therefore, it
is more suitable to design an adaptive threshold [22]. From
Equations 2 to 4, we can conclude that (R +M−1) frame
(0.39 s for fixed R = 30 and M = 10) information is
needed to acquire the first LSFM feature value. In our
implementations, the initial 1.39 s of the test signal x(n) is
always assumed to be non-speech. From this 1.39 s of x(n),
100 realizations of LN can be collected and saved to ψINL.
The threshold is initialized to be
THRINL = min(ψINL). (8)
To update the threshold at the mth frame, we used two
buffers ψS + N and ψN. ψS + N stores the LSFMmeasures
of the last 100 long windowending at themth frame which
was decided as containing speech; similarly, ψN stores
the LSFM measures of the last 100 long window ending
at the mth frame which was decided as including nonspeech
information only. The adaptive threshold for the
mth frame is then updated as:
THR(m) = λ×min(ψS+N)+(1−λ)×max(ψN), (9)
where λ is the parameter of the convex combination. We
experimentally found that λ = 0.55 results in the maximum
accuracy rate in VAD decisions over the TIMIT
training set.


Conceptually you should split this into two steps

Up-sample by a factor of 3, i.e. x = [a 0 0 b 0 0 ...]. This results in a 3 times periodic repetition of the spectrum.
Convolve with a rectangular pulse of length three, i.e. h = [1 1 1]. This crates time stretched sequence you want. Convolution in time is multiplication in frequency, your three copies of the original spectrum get multiplied with a $sinc$ function.


Approach to Determine Short-Term Delay Between Two Signals

Unless I misunderstood completely, you are looking for the lag between the two waveforms red and green over smaller windowed intervals of that function as the window is moved over the sequence? In this case it can simply be done with a moving window using the following function in Matlab:

[corr, lag] = xcorr(a,b)

This will return the cross-correlation as corr, and the lag variable is the delay between a and b associated with each cross correlation.

Given you want to do this as a moving window, that can be achieved by setting a window size and index, and then index the window through your data:

Note this may be pseudo code below as I haven't tested the actual code written, but should give you the idea, as shown is a moving window of size K indexing 1 sample at a time, and will compute the relative delay between the waveforms within the window view.

x= green-mean(green);    % remove mean
x=x(:)';                 % ensure row vector
y = red-mean(red);
y=y(:)';
[~,n]=size(x);

for shift = 1:n-K    % ensure x and y same size
 xwin= x(shift:shift+K); % K is the window size
 ywin= y(shift:shift+K);
 [corr,lag]= xcorr(xwin,ywin);
 disp(shift)= lag(corr==max(corr));
end

If this is on track to what you are looking for, I would also consider normalizing the correlation by the standard deviation of the two waveforms within each window, thus giving you the correlation coefficient. This could allow for the option of optimally weighting the results for estimation of overall delay (or ignoring regions where the match was poor etc).


Different results for seperable filtering v/s direct filtering on 2-D image in MATLAB

Turns out that when doing filtering by columns, i.e. for the following section

for j = 1:size(I,2)
   I_seperable(:,j) = imfilter(double(I2(:,j)),g1,'circular','same','conv');
end

we need to make sure the filter is aligned along columns. Changing g1 to g1' makes the difference almost zero as seen below.

ans =

   1.9782e-11


Intuition for how a Fourier transform handles signals where the frequency varies with time?

Each frequency burst is equivalent to a rectangular window on an infinite sinusoid.  A rectangular window in the time domain is the same as circular convolution with a Sinc function in the frequency domain.  Thus, you end up with a spike with side humps (the spike for the sinusoid, and the humps due to the convolution with the transform of the rectangular window) for each frequency.

Then, since the FFT operation is linear, the 3 sine bursts or, equivalently, the 3 spiky hump groups, just sum additively in both the time and the frequency time domain, respectively.

You are correct in that the phase information is important for reconstruction.  The FFT phase results includes the information about the location of each (non-overlapping in your case) rectangular window for each sine burst.  Twist (or untwist) each spiky hump group's phase in the frequency domain, and you will move the burst's location (the rectangular window) in the time domain.


Discrete-time system: getting the true frequency from normalized frequency

The maximum frequency of a discrete-time signal is half the sampling frequency. A signal with maximum frequency in discrete time is an alternating sequence, and since its period is two sample intervals, its frequency is $f_s/2$. All frequencies in the discrete domain are normalized by that maximum frequency.

This implies that a discrete-time system's behavior can be described independently of the sampling frequency. E.g., if you have a low pass filter, the cut-off frequency is fixed in terms of relative frequency. The actual cut-off frequency in Hertz depends on the sampling frequency.

Often you will encounter the normalized frequency in radians:

$$\omega=2\pi\frac{f}{f_s}\tag{1}$$

where $f_s$ is the sampling frequency, and $f$ is the frequency in Hertz. Note that the frequency response of a discrete-time system is $2\pi$-periodic with respect to $\omega$, i.e., with respect to $f$ it is periodic with period $f_s$.


Why does downsampling stretch a signals frequency response and upsampling shrink and create images of a signals frequency response?

I try to think of it this way: In the time domain, when downsampling occurs, the signal gets compressed; while on upsampling, the signal gets stretched.

Then, from the Fourier transform we know that time stretching means frequency compression, and vice-versa.

It may not be a rigorous answer, but hope this helps!


Is there a strategy for discrete control of a system with dynamics near sample rate?

Depending on your plant you might be able to do something. If the plant is unstable, but dominantly linear, you should be able to get a very accurate model using exact discretization (zero-order hold). Either you will be able to stabilize your plant using PID, or PID will be of too low order, in which case you need to switch to a control law matching you plant order. The H-infinity synthesis framework (robust control) will do that for you. Here you can also include uncertainty weights to limit control law gain, which might help esp. if your system is not linear. Depending on the uncertainty you might not be able to find a feasible control law for stabilization, however...

If your plant is stable, or stabilizable, you can further improve performance using feed-forward control. Here the sampling time should not be of much concern, and you can use something like zero phase error tracking.


Posting a standard is actually a good point.
After reading several definitions i've figured out that i have to keep the symbol (and thus the phase) if I want to transmit 1 and reverse it if i want to transmit zero.
(see: http://www.arrl.org/psk31-spec, PSK-63 has just higher speed, about 63 baud.)
if(ctr==(11025/63)){
    ctr=0;
    c=message[subsym]-'0';
    subsym++;
    if(subsym>=strlen(message)){
        subsym=0;
    }
    if(!c) // here instead of just assigning bit=symbol
    sym_in = !sym_in;
}


Are SOS structures better for linear phase FIR filters?

well, i see little value in factoring an FIR into quadratics and implementing them in cascade.  it won't be computationally cheaper and it's not better from a quantization noise POV if your FIR has access to an accumulator that is double width.

perhaps it will help with limiting coefficient range so that it's less likely that you have non-zero coefficients (at the tails of the FIR) with values of, say, $10^{-8}$ the size of other coefficients (in the middle).  there are other, more efficient ways of dealing with that problem (like computing the FIR from the outside in, doing the tails first and the middle last).


Understanding SURF Features Calculation Process

What's SURF?

In order to correctly understand what is going on, you also need to be familiar with SIFT: SURF is basically an approximation of SIFT. Now, the real question becomes: what's SIFT?.

SIFT is both a keypoint detector and a keypoint descriptor.
In the detector part, SIFT is essentially a multi-scale variant of classical corner detectors such as the Harris corner, and that has the ability to auto-tune the scale. Then, given a location and a patch size (derived from the scale), it can compute the descriptor part.

SIFT is very good at matching locally affine pieces of images, but it has one drawback: it is expensive (i.e., long) to compute.
A large amount of time is spent in computing the Gaussian scale-space (in the detector part), then in computing histograms of the gradient direction (for the descriptor part).

Both SIFT and SURF can be seen as difference of Gaussians with automatic scale (i.e., Gaussian sizes) selection. This, you construct first a scale-space where the input image is filtered at different scales. The scale-space can be seen as a pyramid, where two consecutive images are related by a scale change (i.e., the size of the Gaussian low-pass fiéter has changed), and scales are then grouped by octaves (i.e., a big change in the size of the Gaussian filter).

In SIFT, this is done by repeatedly filtering the input with a Gaussian of fixed width until the scale of the next octave is reached.
In SURF, you do not suffer any runtime penalty from the size of the Gaussian filter thanks to the use of the integral image trick. Thus, you compute directly the image filtered at each scale (without using the result at the previous scale).

The approximation part

Since computing the Gaussian scale-space and the histograms of the gradient direction is long, it is a good idea (chosen by the authors of SURF) to replace these computations by fast approximations.

The authors remarked that small Gaussians (like the ones used in SIFT) could be well approximated by square integrals (also known as box blur).
These rectangle averages have the nice property to be very fast to obtain thanks to the integral image trick.

Furthermore, the Gaussian scale-space is actually not used per se, but to approximate a Laplacian of Gaussians (you can find this in the SIFT paper).
Thus, you do not need just Gaussian-blurred images, but derivatives and differences of them. So, you just push a bit further the idea of approximating a Gaussian by a box: first derive a Gaussian as many times as needed, then approximate each lobe by a box of the correct size. You will eventually end up with a set of Haar features.

Increment by 2

This is just an implementation artifact, as you have guessed. The goal is to have a central pixel.
The feature descriptor is computed with respect to the center of the image patch to be described.

Middle region

When going from a black ray to a white ray, you have something like $\sum_{\text{all pix in column}} \partial x = A$. Then, going from white to black, you have the opposite sum: $\sum_{\text{all pix in column}} \partial x = -A$. Thus, you have a small $\sum \partial x$ for the window, but a higher sum of the magnitudes.

Magic number

The first scale is obtained by applying a blur with $\sigma = 1.2$ (or 1.4 in some papers). This is because a natural (real) sharp image can be considered as being the result of the convolution of an ideal (without aliasing) image with a blur kernel of width $\sigma = 1.2$. I can't really remember where it comes from, but it was also explicitly studied in Guoshen Yu's work on A-SIFT, so you might check this page.


It really depends on which application you have...the algorithms for speech recognition are different from speech indentification for example. Based on what you need, just look for an open source library for Java. Getting into the theory is definitely desired but as a first step you can just look for a software block that has inputs/outputs, later you can drill into it.

Here is some info on speech recog: https://stackoverflow.com/questions/3018418/open-source-simple-speech-recognition-in-c-in-windows


Why does the Fourier series not include phase information?

(The following is the answer to the original version of the question asking about the Fourier Transform. But exactly the same argument also applies to Fourier series).

They are. Assume that $F(\omega)$ is the Fourier transform of the signal $f(t)$:

$$F(\omega)=\int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt$$

In general, $F(\omega)$ is a complex-valued function and can be written in polar form

$$F(\omega)=|F(\omega)|e^{j\phi(\omega)}$$

where $\phi(\omega)$ is the phase of $F(\omega)$.
The signal $f(t)$ can be written in terms of $F(\omega)$ by virtue of the inverse Fourier transform:

$$f(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{j\omega t}d\omega=
\frac{1}{2\pi}\int_{-\infty}^{\infty}|F(\omega)|e^{j(\omega t+\phi(\omega))}d\omega\tag{1}$$

From (1) it can be seen that $f(t)$ is represented by complex sinusoids with varying amplitudes $|F(\omega)|$, phases $\phi(\omega)$, and frequencies $\omega$.


Why is In-phase Quadrature sampling not used to record and store digital audio?

It is because the audio signals are real and already at baseband. In contrast radio frequency signals are often represented as complex numbers once they are brought back to baseband. Real signals can be represented as a single stream of real numbers, while for complex numbers two streams of real numbers are required to represent them (as in $I+jQ$).

When the upper sideband of a radio waveform at a given carrier frequency does not match the lower sideband of the waveform in magnitude and opposite phase, then the waveform when brought to baseband must be represented with a complex signal (because at baseband the positive and negative frequencies will not be complex conjugate symmetric, which is the requirement for any real signal). A simple example of this is Quadrature Phase Shift Keying and Quadrature Amplitude Modulation where samples are mapped to locations on the complex (IQ) plane as part of the modulation.

It is helpful to have a full grasp of a complex frequency represented as $e^{j\omega t}$ vs $cos(\omega t)$ and what happens when you translate either in frequency to a real carrier: in the case of the exponential only a positive frequency exists so when translated using a real carrier frequency will be a single real tone at a higher frequency than the carrier (upper side band) and similarly $e^{-j\omega t}$ will be a single real frequency at a lower frequency than the carrier (lower side band). The sinusoid in contrast contains both exponentials (see Euler’s identify) and therefore will have both sidebands when translated. This shows how the complex form can be used to occupy half the bandwidth which is a motivation for doing so (in fact the higher modulation orders mentioned offer better spectral efficiency and are complex waveforms in order to do this.)

Similarly the audio signal could be made complex if it was to be converted to a carrier for the same reason (for example single-sided AM is one variant of an analog modulation that can be used with audio in complex form.)

But ultimately an audio signal that is going to one speaker is represented by a single stream of real signals so as explained earlier is not complex — but between the microphone and the speaker we can choose any way to encode the data including complex signals.


literature related to Spectral resolution

I would recommend
Array Signal Processing - Johnson, Dudgeon. This book covers classical spectral estimation, Minimum Variance Distortionless Response (MVDR), Linear Prediction, and subspace methods (e.g. MUSIC and ESPIRIT). It provides examples of the resolution capabilities between these methods. The are quite a number of good references in this book if you require additional detail.

MUSIC and ESPIRIT type algorithms are dicussed in most advanced DSP texts:
Advanced Digital Signal Processing, Proakis, Rader, Ling, Nikias
Discrete Random Signals and Statistical Signal Processing, Therrien

For parametric methods:
Spectral Analysis - A modern Perspective (Kay, Marple) Proc. IEEE Vol 6, No 11, 1981
Modern Spectral Estimation: Theory and Application, Steven Kay
Digital Spectral Analysis - S.L. Marple
SPECTRAL ANALYSIS OF SIGNALS - P Stoica and R Moses  PDF is available here

The books by Marple and Kay tend to focus on algorithms i.e. implementations of various MA, AR, and ARMA approaches (e.g. fast Lattice filter implementations) rather than the performance of a basic approach

Compressive Sensing / Sparse Reconstruction techniques can also be used. There are too many references here to list. There are a few books by these authors Michael Elad, Yonina Eldar, and Holger Rauhut. These techniques are often quoted as "having super resolution properties"

The algorithms are often hard to compare directly because they all have twiddle factors which affect their performance.  For example
 1. For MA, AR, ARMA, what order of model to use?
 2. For Compressive techniques - how closely to space the template vectors? What level of regularization to use? How many frequencies are there? How many samples do we have?


Why do we read that Hilbert transform can be used for envelope detection?

Using a Hilbert transform with a sufficiently narrow-band signal, the low pass filter may not be needed even for continuous envelope estimation.


How to get Harris corners that are present in watershed segmented Bitwise image(road part)?

You have a mask of the road were pixels of the road have value 1 and non road pixels as value 0.
Also you have a list of Harris corners were each corner as a coordinate Xi,Yi

For each Xi,Yi i corners
If mask[Xi,Yi] == 1
keep corner


As you have observed, instantaneous signal amplitude is an easy to understand, one-dimensional quantity, whose averaged absolute value is loosely associated with perceived loudness, and which can be easily altered by multiplying your signal by a constant.

Frequency is a very different affair. If we lived in a world where the audio signals were made of sustained pure tones, it would make sense to ask "what is the frequency of this sound?". But things are more complicated:

Stationary speech or music signals can be described as a sum of pure tones. Thus, it makes more sense to ask "what is the signal amplitude at this particular frequency?". This is why we look at spectra - graphs showing signal amplitude as a function of frequency. Note that frequency is not the plotted  quantity, it is the $x$ axis!
Speech or music signals are not stationary: some components of the signal decay over time or are modulated, some appears... Thus, it makes even more sense to ask "what is the signal amplitude at this particular frequency and particular time frame?". This is why we commonly see for audio analysis a representation called "spectrogram" - showing the measured signal energy as a function of time and frequency (again, note that frequency is the axis, not the measured quantity).
In one of your comments you say "as well as be able to read the frequency data in each sample". Frequency is a rate of change, and to observe a change, we need a collection of samples - it makes no sense at all to do any frequency analysis on a single sample, and if you use Fourier analysis (which is just one of the ways of estimating a sum of sine waves from a signal), you'll need more and more samples to get a more and more precise resolution on the frequency axis. But remember that sounds change over time, and that performing a Fourier Transform discards timing information (the magnitude spectrum of two simultaneous 100 Hz and 300 Hz tones is the same as the spectrum of a 100 Hz tone followed in time by a 300 Hz tone). Thus, if you care about the temporal position of events, you will need to perform your analysis on slices of the signal, computing a FFT on overlapping blocks of the input signal - this is called the short-term Fourier transform (STFT). You'll have to deal with resolution problems here - the smaller your windows, the more accuracy you have on the time axis, the less accuracy you have on the frequency axis. The bigger the analysis windows, the less accuracy you have on the time axis, the more accuracy you have on the frequency axis. Again, and again, note that frequency is not what you estimate - your STFT is just answering the question "What is the signal energy in this particular frequency range and this particular time frame?".

I hope these will make you realize that there's no such thing as "the frequency" of an audio signal. Frequency is an axis, "where we look at things", not a measured quantity, "what we measure".

Note that there are some one-dimensional quantities homogeneous to frequencies (and thus expressed in Hz), that can be used to characterize stationary (stable) sounds. For generic audio or music, these quantities will vary over time:

Fundamental frequency, is the inverse of the period of the sound. It relates to the perceived musical pitch (note) of the sound. This quantity is undefined for aperiodic sounds. Pitch-shifting algorithms will alter this quantity, to make a recording sound as if it was played/sung higher or lower on the musical scale.
Roll-off frequency, is the frequency under which a given fraction (say 95%) of the energy of a sound is contained. This quantity can be affected by low-pass / high-pass filtering, to make a recording sound "brighter" or "duller". It is one of the many dimensions of the perceived timbre of sound.

Maybe some of these are the quantities you really care about?

One final word about the FFT - it is the hammer in a signal processing version of "if all you have is a hammer, everything looks like a nail". It is an extremely useful tool, but it is not the best fit for many audio signal processing/analysis problems, and more often than not it is used only as an intermediary processing or computational step. In particular, for some of the questions you might be interested in:

Estimating the fundamental frequency might an involve a FFT, but the spectrum is not the right thing to look at, because of the "missing f0" phenomenon (a signal might be perceived as being a 263Hz C3 note, without having energy at this frequency). Don't make the mistake of trying to estimate the musical pitch of a signal by looking at the peak in its spectrum.
Computing a FFT, messing with the numbers, and performing the inverse FFT is a poor way of filtering a signal. This should be done with digital filters, available in scipy.signal.lfilter. You'll find formulae to compute IIR filter coefficients for a variety of musically useful second-order filters here.


Pitch shifting with bins creating large feedback

...when I do this, a large amount of noise relative to the size of "n" is generated (spectral leakage maybe?) The smaller the size n, the larger the amount of noise

No. This is not "spectral leakage". The flutter that you hear comes from the fact that you are pitch shifting each signal frame without taking into account the phase shifting caused by the pitch shifting of the earlier blocks.

To put it in a different way: When you work with an audio recording, the "zero" for time is the first sample and that is the starting point when evaluating phase. By pitch shifting block $n$, you alter its phase (relative to the first captured sample) by some $d\theta_n$ but when you pitch shift block $n+1$, you are not taking into account the $d\theta_n$ you just applied. Therefore, subsequent blocks are out of phase.

What you hear at play-back are those phase discontinuities.

  How do I implement the pitch shifting on a bin by bin basis while minimizing the noise/leakage?

Exactly how you are going to take this into account depends on the algorithm used. If the algorithm is a linear phase algorithm, then you might get away easily by correcting the phase for a constant term. Another thing you can do is to process the audio in overlapping blocks (for example, one block from 0 to 99 and the next from 60 to 159, rather than 0..99 followed by 100..199), still as if they were independent of each other. The final result is produced by averaging over the overlapping region. This does not correct the "problem" but it softens the phase discontinuities. By tuning the overlap (and possibly using a low pass filter too), you can supress those "jumps" to the point that they are not too noticeable. In a real-time setting, the overlap can be implemented as preserving a proportion of the currently playing output buffer which you average with the right segment of the next processed frame as the final step.

But also, saying this, you might be able to find an "online" version of a phase-vocoder with all these considerations implemented.

Hope this helps.


IEEE 802.11 is not related to 5G NR, so consideration for WLAN networks following 802.11ay don't apply to cellular infrastructure.

I don't know what "ad-hoc" modi 5G NR brings; I don't think it's any.

So, 5G doesn't do any collision avoidance, because it's not a problem it has: it's a cellular standard that defines that the basestation / cellular infrastructure defines who sends when, where. The mobile network operators own the spectrum, so there's by definition no uncooperative users.

Generally, whereas interference mitigation is in fact a non-negligible topic for coordinated microwave networks, with these deci- to centimeter waves still taking many non-optical paths and receiver sensitivity being high enough for that to create "unwantedly far reach", you can basically ignore that for cellular mmWave applications, where you can be pretty sure that you, as the owner of the used spectrum, will have a relatively easy time making sure that your own femtocell base stations don't lead to mutual interference: The propagation is pretty much optical, so, the potential numbers of base stations a single handset might reach is pretty limited by lines of sight and a pretty small radius.


What does "A Target with Zero Doppler" Mean?

In Pulse Doppler Radar, the Doppler frequency that a target's return signal exhibits can be used to calculate the Radial Velocity of target.  That means, the return frequency of a pulse will be different than the frequency of the pulse that was sent out depending on how fast the target is moving toward the source of the pulse.  If the target's "Radial Velocity" is zero, zero doppler will result, i.e. the resulting doppler caused by the target's Radial Velocity will be 0 Hz.

Wikipedia shows the calculation of Doppler to be:

$$
F_{Doppler} = \frac{2 * F_{Transmit} * V_{Range} }{c}
$$

$c$ being the Speed of light.

So, if the target's Radial Velocity is zero (not moving towards or away from the observer), $F_{Doppler}$ will also be zero.

https://en.wikipedia.org/wiki/Pulse-Doppler_radar


Struggling to find a good definiton for Gabor-Descriptor

In this context, $j$ is the imaginary unit: $j = \sqrt{-1}$. The complex exponential function represents what they call an "oriented sinusoidal grating." The fact that $j$ is the imaginary unit in this case can be verified by noticing the relationship between $(1)$ and $(2)$ in Section 2; this is an example of the frequency-shifting property of the Fourier transform. That is, multiplying the Gaussian function $h(x,y)$ by the complex exponential function results in a shift in the Fourier transform of the result along the axis that the sinusoid runs across (the $x$ axis in your example).


Upsample / Downsample a Signal from 1920.93Hz to 1920Hz

Assuming the signal was reasonably band-limited before sampling, you can just use a high quality interpolation, such as a Sinc interpolation kernel, to calculate new equally spaced points at any lower density.

For small frequency deltas, even a long interpolation kernel or high degree polynomial interpolation can be a lot more efficient than naively upsampling and then downsampling.


why is it so hard to simulate up and down position on headphones?

Front/back is actually very hard too, at least for stationary simulation without head tracking.

The main reason is simple: Left/Right is done by looking at the differences between the two ear signals, i.e. interaural level differences and interaural time differences. If a source is located to the left, the sound will arrive earlier on the left ear and will be louder (depending on frequency).

In the median plane (front, up, back down), the ear signals are almost identical. The only information available are pinna cues, shoulder reflection, interaural changes due to micro rotations of the head, small asymmetries, contextual information, etc. . That's a lot less robust so the localization acuity in general is a lot worse. It's also harder to model since it depends a lot on the individual human and one-size-fits-all doesn't work well.


How to perform a Rubberband-Correction on spectroscopic data?

This can be easily done in R or Python. There are well-tested functions available, so you don't have to worry about any boundaries or nuances. Moreover, both are free and popular among scientists.

Solution for R

There is a special package to handle spectral data, called hyperSpec. The rubberband baseline correction is already implemented there (function spc.rubberband). All details are highlighted in the documentation. The usage is pretty simple:

require(hyperSpec)
spc <- read.spe("paracetamol.SPE")
baseline <- spc.rubberband(spc)

corrected <- spc - baseline

Solution for Python

There is no (to the best of my knowledge) out-of-the box solution for python, but you can use scipy.spatial.ConvexHull function to find indices of all points that form a convex hull around your spectrum. Suppose that the spectrum is contained in x and y arrays:

import numpy as np
from scipy.spatial import ConvexHull

def rubberband(x, y):
    # Find the convex hull
    v = ConvexHull(np.array(zip(x, y))).vertices

Array v contains indices of the vertex points, arranged in the CCW direction, e.g. [892, 125, 93, 0, 4, 89, 701, 1023]. We have to extract part where v is ascending, e.g. 0–1023.

    # Rotate convex hull vertices until they start from the lowest one
    v = np.roll(v, -v.argmin())
    # Leave only the ascending part
    v = v[:v.argmax()]

    # Create baseline using linear interpolation between vertices
    return np.interp(x, x[v], y[v])

Now the baseline is corrected like this:

y = y - rubberband(x, y)


Let's give it a shot.

We know that during sampling, we obtained a discretized version, $x(nT_s)$ of a continuous time signal $x(t)$ sampled every $T_s$ seconds. I will denote their corresponding Fourier Transforms as $X_{DT}(e^{j\omega})$ and $X_{CT}(j\Omega)$, respectively, according to the majority of the bibliography. Indices denote discrete time (DT) and continuous-time (CT) Fourier Transforms.

The Discrete time Fourier Transform (DTFT) of the sampled signal is given by $$X_{DT}(e^{j\omega}) = \frac{1}{T_s}\sum_{l=-\infty}^{+\infty}X_{CT}\Big(j\Big(\frac{\omega}{T_s} - l\frac{2\pi l}{T_s}\Big)\Big)$$

Assuming that there is no aliasing (like in your case), over one period this turns into $$X_{DT}(e^{j\omega}) = \frac{1}{T_s}X_{CT}\Big(j\frac{\omega}{T_s}\Big)$$

Due to the $2\pi$-periodic nature of the DTFT and the fact that we're dealing with real-valued signals, we can rewrite this as $$X_{DT}(e^{j\omega}) = \left\{\begin{array}{ll}
\frac{1}{T_s}X_{CT}\Big(j\frac{\omega}{T_s}\Big), & 0 \leq \omega < \pi \\
\frac{1}{T_s}X_{CT}\Big(j\frac{\omega-2\pi}{T_s}\Big), & \pi \leq \omega < 2\pi \\
\end{array}\right.$$

The Discrete Fourier Transform (DFT) is an $N-$point uniformly sampled version of the DTFT, so $$X_{DFT}[k] = X_{DT}(e^{j\omega})\Big|_{\omega = 2\pi k/N}$$ for $0 \leq k\leq N-1$.

Combining the last two equations we have $$X_{DFT}[k] = \left\{\begin{array}{ll}
\frac{1}{T_s}X_{CT}\Big(j\frac{2\pi k}{NT_s}\Big), & 0 \leq k < N/2 \\
\frac{1}{T_s}X_{CT}\Big(j\frac{2\pi (k-N)}{NT_s}\Big), & N/2 \leq k \leq N-1 \\
\end{array}\right.$$

Except for the constant $N$ in your results - that I have failed to extract - I think what I wrote will help you.


Resample an array in real-time (application: pitch-bend)

Audio-optimized resampling is done by using anti-alias filtering and interpolation kernels with good stop-band attenuation and low pass-band ripple/distortion.  On current PC and mobile processors (the kind that can easily compute hundreds of transcendental functions per sample time), these filter and interpolation kernels can not only be precomputed and cached in poly-phase lookup tables, but for certain window types, (re)computed in real-time per sample.

Thus, resampling can be done in real-time by computing the needed interpolation kernels and interpolating each sample as needed at the required sample delay from the last for the current sample rate.

Note that changing the sample rate too quickly may result in FM modulation skirts that may need additional filtering to reduce aliasing artifacts.


From transfer function to differential equation

Given a transfer function
$$G_v(s) = \frac{k_v}{1 + sT} \tag{1}$$
the corresponding LCCDE, with $y(t)$ being the solution, and $x(t)$ being the input, will be
$$ T ~\dot{y}(t) + y(t) = k_v ~x(t) \tag{2} $$
Your formulation replaces $x(t)$ with a unit-step $u(t)$, and $y(t)$ with $x(t)$, yielding
$$ T ~\dot{x}(t) + x(t) = k_v ~u(t) \tag{3} $$
or equivalently
$$ \dot{x}(t) + \frac{1}{T} x(t) = \frac{k_v}{T} ~u(t) \tag{4} $$
Eq.4 represents a typical first order, constant coefficient, linear, ordinary differential equation (abbr LCCDE) whose solution procedure is as follows:
First, find the homogeneous solution to the Eq.4 with RHS being zero, as
$$ x_h(t) = A e^{-t/T} \tag{5} $$
where $A$ is an unknown coefficient to be determined by the initial condition $x(0)$.
Second, find the particular solution of Eq.4 with the RHS being $u(t)$. For this you can use the method of undetermined coefficients, by assuming a solution of the type
$$ x_p(t) = B u(t) \tag{6}$$
an letting $t \to \infty$, will yield in Eq.4 as
$$ 0 + \frac{B}{T} = \frac{k_v}{T} \tag{7} $$
where you find $B$ as
$$ B = k_v \tag{8}$$
Then you find the complete solution by adding the homogeneous part to the particular part, and solving for the unknown constant $A$ from the initial condition $x(0)$
$$ x(t) = x_h(t) + x_p(t) = A e^{-t/T} + k_v ~~~~,~~~~ t >  0\tag{9} $$
The initial condition $x(0)$ is given as $0$ in the problem, hence
$$ x(0) = 0 = A + k_v \implies   A = -k_v \tag{10} $$
And the overall solution of the initial value problem in Eq.4 with $x(0)=0$ becomes:
$$ x(t) = -k_v e^{-t/T} + k_v = k_v (1 - e^{-t/T}) u(t) \tag{11} $$
Note that we look for a causal solution, unless otherwise stated, hence assume $ t>0$.


It would be more accurate to say that the MFCC "capture" or "represent" the spectral envelope - in the sense that two signals with a similar spectral envelope will have a similar sequence of MFCC.

The MFCC coefficients themselves are not supposed to "look" like the spectral envelope when plotted. After all the spectral envelope is a "big" vector (as long as the size of an analysis window), while the MFCC are designed to be short (a handful of coefficients). The spectral envelope is highly redundant (this is a very smooth function, consecutive values are correlated), while the MFCC are, by design, very little correlated.

One could visualize the spectral envelope from the MFCC by inverting all steps in the process except the first one (note that you probably omitted one last step: truncate the vector and keep first 13 or 20 coefficients). The reverse operations would thus be: pad with zeros, compute inverse DCT, take exponential, and compute a matrix product with the matrix containing the triangular bin responses. This will yield a smoothed version of the spectrum that will show what the MFCC "see" of the spectrum.

Finally, unless "example.wav" is very short, there is something wrong in your code. Audio signals are analyzed in short overlapping frames rather than wholly. It's likely that you might instead want to look at spectrograms; and at the MFCC as a 2D matrix (imagesc).


How to convert between function forms of single frequency signals


Myth: DTFT is Sinc-interpolated DFT.

Problem with the above statement: Sinc is not $2\pi$-Periodic function, but all DTFTs are.

Correct Answer:

Theoretical, Continuous-$\omega$ $2\pi$-Periodic DTFT can be obtained by continuous Lagrangian-interpolation of the DFT Samples. So that the values at $\omega = 2\pi k/N$ will be the DFT Samples $X[k]$ for $k=0,1,...,N-1$ and the Interpolation-function's zero-crossings are at $2\pi k/N$.

In other words, DTFT will take same values at roots of unity as DFT Samples, but it will be a smooth interpolation of the DFT at other values of digital frequency $\omega$.

Mathematically, let $x[n]$ be N-Length sequence and $X[k]$ be it's N-point DFT. Now, DTFT is defined for infinite length sequences. So, lets derive DTFT of a finite length $x[n]$.

$$X(e^{j\omega}) = \sum^{N-1}_{n=0} x[n]e^{-j\omega n},$$ now write IDFT of $X[k]$ in place of $x[n]$.

$$X(e^{j\omega}) = \sum^{N-1}_{n=0} (1/N.\sum^{N-1}_{k=0}X[k]e^{j 2\pi k/N})e^{-j\omega n},$$ now bring summation w.r.t n inside,

$$X(e^{j\omega}) = \sum^{N-1}_{k=0} (X[k] (\sum^{N-1}_{n=0}1/N e^{j 2\pi k/N}e^{-j\omega n}))$$

$$ = \sum^{N-1}_{k=0} (X[k] (\sum^{N-1}_{n=0}1/N e^{-jn(\omega - 2\pi k/N)}))$$

So, basically, $$X(e^{j\omega}) = \sum^{N-1}_{k=0} X[k]. \Lambda(w - 2\pi k/N),$$ where $$\Lambda(w) = 1/N \sum^{N-1}_{n=0} e^{-jn\omega}.$$

What this means is each sample of $X[k]$ is multiplied a $2\pi k/N$ shifted copy of $\Lambda(\omega)$ and added together. Basically, $X[k]$ is interpolated by a continuous-$\omega$ and $2\pi$-Periodic function $\Lambda(\omega)$. And this function is not a Sinc-function but something else. Sure it looks like Sinc and it will approach to Sinc in the limit.

Further , $$\Lambda(\omega) = \frac{1}{N}*e^{-j\omega\frac{(N-1)}{2}} \frac{\sin(N\omega/2)}{\sin(\omega/2)}.$$ Plotting this function in $[-\pi,\pi]$ is below:

>> w = -pi:0.0001:pi;
>> y = 1/64 * sin(w*64/2)./sin(w/2);
>> plot(w,y)

I repeat, it is not a sinc interpolation. Sinc is not $2\pi$-Periodic function. There is no way we can get a DTFT by interpolating with sinc.

Practically, you can get DTFT by interpolating with the MATLAB snippet I have provided which approximates $\Lambda(\omega)$ function.

What you can check yourself is extending the above plot to $[-4\pi:4\pi]$ and see that it indeed is Periodic function.


The expectation of $x[n]$ is

$$E(x[n])=E(A+w[n])=E(A)+E(w[n])=A+0=A\tag{1}$$

because $A$ is deterministic and the noise $w[n]$ is assumed to have a mean of zero. So in that formula you sum up $N$ terms that are equal to $A$, and then you divide that sum by $N$, leaving you with the value $A$.

As for the expectation of the estimated value $\hat{A}$, note that $\hat{A}$ is a random variable because it necessarily depends on the data $x[n]$, which are random, so computing its expected value (and its variance) makes sense.


It's a bad translation of step response, i.e., the filter's response to a unit step $u[n]$ at the input. Note that the step response is the convolution of the unit step with the filter's impulse response:
$$s[n]=\big(u\star h\big)[n]=\sum_{k=0}^nh[k],\qquad n\ge 0\tag{1}$$
where I've assumed causality, i.e., $h[n]=0$ for $n<0$.
From the given step response and Eq. $(1)$ it should be easy to derive the impulse response $h[n]$. You can also immediately draw a conclusion about its length. Then, by considering the symmetry of $h[n]$, I'm sure you know how to derive the group delay of the filter.


Smoothing "font characters" outlines in bitmap image

I think the open-source command-line program potrace, might do what you want.  It converts bitmaps to bezier curves and has a bunch of options allowing you to trade off smoothness and accuracy.

The open-source Inkscape vector (svg) editor has potrace built in (under the Path->Trace Bitmap menu option.)  The result of applying your example in Inkscape, (I used the mode Inkscape calls "Brightness cutoff" with the default threshold of .45, since your image was already black-and-white rather than grayscale) scaling the size of the resulting bezier curves up a bit, and then having Inkscape output a bitmap is:

One of the cute side-effects of converting to bezier and then back to grayscale is that you get nice anti-aliasing around the curves.

With grayscale input images I find that the "Edge detection" mode works a little better.  There are also options for choosing less or more smoothing, and less or more speckle elimination.


Confusion in Channel Encoding and Convolutional Encoding

Like other comments, I don't really understand your question. What I am trying to do is to list basic things so that other ones can suggest edit because I find it is much easier to write in answer part :). Please let me know if it makes you less confused.

The general communication system is

At transmitter, before modulator is digital, modulator is the block to convert from digital world to analog world because physical signal is transmitted in analog (at least at our normal scale). Modulator converts bit sequences to voltage (NRZ), amplitude and phase of sinusoid (QAM), frequency (FSK), phase of of sinusoid (MSK, GMSK), ...

Signal is transmitted in physical channel always suffers degradation (noise, fading, ...) which is random. To combat this negative effect, channel coder adds redundant information to help us to recover the correct signal. For example adding 9 supplementary bits instead of sending only one, there is the chance that more than 5 bits are correctly decoded and with majority rule, the transmitted data is recovered. The operation is called channel coding. There are a lot of channel codes : turbo code, LDPC, convolutional code, ...

Thus, channel coding and modulation are two different steps that have different roles. For more details, you need to read the recommended resources in the comments.

It seems that you understand "convolutional coding" as modulation. Maybe it comes from the fact that signal transmitted over channel is modeled as convolution operation ? People usually use "convolutional code" to talk about a specific family of code in channel coding step.

About the source coder block in the figure, on the contrary of channel coder that adds redundant information to combat random degradation, source coder is here to remove redundant information from source. Huffman code, Shannon code, ... are famous. Applications such as winrar, winzip, ... are well-known. Source coder and channel coder have different roles, but they are related in the sense of "optimality" because they both deals with "information". Thus it is not always correct to say that source coder and channel coder are independent. You can search for "joint source channel coding" if you are interested. The reason why the myth exists and an counter-proof can find in this article Source Channel seperation Theorem. Note that it could be an advanced topic, depending on your background.


Which channel coding technique is suitable for wireless SPIHT compressed image: Reed-Solomon or convolutional code?

The channel coding used has little to do with your source coding. So, this depends on how you model your channel (and noise, and interferers), and what kind of errors you're willing to tolerate. Saying "SPIHT-compressed image" answers neither!