Why is it a bad idea to filter by zeroing out FFT bins?

Zeroing bins in the frequency domain is the same as multiplying by a rectangular window in the frequency domain.  Multiplying by a window in the frequency domain is the same as circular convolution by the transform of that window in the time domain.  The transform of a rectangular window is the Sinc function ($\sin(\omega t)/\omega t$).  Note that the Sinc function has lots of large ripples and ripples that extend the full width of time domain aperture.  If a time-domain filter that can output all those ripples (ringing) is a "bad idea", then so is zeroing bins.

These ripples will be largest for any spectral content that is "between bins" or non-integer-periodic in the FFT aperture width.  So if your original FFT input data is a window on any data that is somewhat non-periodic in that window (e.g. most non-synchronously sampled "real world" signals), then those particular artifacts will be produced by zero-ing bins.

Another way to look at it is that each FFT result bin represents a certain frequency of sine wave in the time domain.  Thus zeroing a bin will produce the same result as subtracting that sine wave, or, equivalently, adding a sine wave of an exact FFT bin center frequency but with the opposite phase.  Note that if the frequency of some content in the time domain is not purely integer periodic in the FFT width, then trying to cancel it by adding the inverse of an exactly integer periodic sine wave will produce, not silence, but something that looks more like a "beat" note (AM modulated sine wave of a different frequency).  Again, probably not what is wanted.

Conversely, if your original time domain signal is just a few pure unmodulated sinusoids that are all exactly integer periodic in the FFT aperture width, then zero-ing FFT bins will remove the designated ones without artifacts.


How to compute bit level soft decisions from an M-ary symbol?

In each case, the soft-decision output is just the matched filter's output.

For example, say you transmit the bit sequence $1,1$ with codeword $1,0,1,1$. You transmit these bits using, say, QPSK symbols $0.5+0.5j,0.5-0.5j$. After matched filtering in the receiver, you will get two noisy complex numbers, for example $0.55+0.45j,1.2-0.49j$. A hard-decision decoder will turn those numbers into QPSK symbols and then bits before FEC. A soft-decision decoder will use those numbers without any further processing (except maybe some quantization).


Let me propose a majorizing date for this discovery, at least in written form: Rader, C. M. and Gold, B. (1967) Digital filter design techniques in the frequency domain, Proceedings of the IEEE. Such a discussion appears on page 164, under the name of "frequency sampling filters", involving the combination of "elemental" and "comb filter". They say "The phase versus frequency is exactly linear except for discontinuities of $\pi$ radians. These discontinuities occur where the magnitude response is zero".

A lower bound could be Lerner, R. M. (1964) Band-pass filters with linear phase, Proceedings of the IEEE, which seems to keep analog.


The conjugate product of the two filters form a frequency discriminator with an error output that is proportional to a frequency offset at the input.

This is due to the relationship that the cross correlation in one domain is the conjugate product in the other domain (you can skip past the derivation unless interested):

$$c(\tau) = XCORR[x(t)y(t)]$$
$$c(\tau) = \int_{t=-\infty}^{\infty}x^*(t)y(t+\tau)dt$$
$$C(\omega) = \mathscr F\{y(t)\}= \int_{\tau=-\infty}^{\infty}\bigg[\int_{t=-\infty}^{\infty}x^*(t)y(t+\tau)dt\bigg]e^{-j\omega \tau}d\tau$$
$$=\int_{t=-\infty}^{\infty}x^*(t)\bigg[\int_{\tau=-\infty}^{\infty}y(t+\tau)e^{-j\omega \tau}d\tau\bigg]dt$$

From the shift property of the Fourier Transform, the terms inside the brackets are:

$$\int_{\tau=-\infty}^{\infty}y(t+\tau)e^{-j\omega \tau}d\tau = Y(\omega)e^{j\omega t}$$

Therefore:
$$C(\omega) =\int_{t=-\infty}^{\infty}x^*(t)\bigg[Y(\omega)e^{j\omega t}\bigg]dt$$

$$ = Y(\omega)\int_{t=-\infty}^{\infty}x^*(t)e^{j\omega t}dt  = Y(\omega)\bigg[\int_{t=-\infty}^{\infty}x(t)e^{-j\omega t}dt\bigg]^*$$

$$ =Y(\omega)X^*(\omega) $$

The above shows that the cross correlation of two functions x(t) and y(t) in the time domain is equivalent to the conjugate multiplication of their Fourier Transforms (in the frequency domain). Similarly, and this is the important point to explaining the operation of the frequency discriminator, is that the converse is also true: the cross correlation of two functions in the frequency domain is equivalent to the conjugate multiplication of the transform of the same two functions in the time domain.

Consider the following as example spectrums of the "Matched Filter" and "Frequency Derivative Matched Filter" in the block diagram given by the OP. Notably the Frequency Derivative Matched Filter is a Band-Edge Filter.

When there is no frequency error, the input spectrum would be aligned with both of these resulting in the following spectrums at the outputs of these filters. Notably, observe that convolving these two output spectrums will result in a zero result specifically due to the positive and negative balance in the Frequency Derivative Matched Filter Spectrum!

When there is a positive frequency error in the input signal as depicted in the graphic below, the input spectrum would be offset by the frequency error resulting in the following spectrums at the output of the two filters. In this case specifically, convolving the two output spectrums would result in a negative value corresponding to the frequency error.  Similarly a negative frequency error would result in a positive output value after the convolution in frequency of the two spectrums.

The above shows how the convolution of a Matched Filter Spectrum and a Frequency Derivative Matched Filter Spectrum results in an output value that is the error term for a frequency discriminator. Additionally we showed that convolution in frequency is equivalent to the complex conjugate multiplication in time. Thus the complex conjugate multiplication of the two filters in the block diagram form a frequency discriminator as the error signal into the loop.

The frequency discrimintor would operate as described without the scaling by the Signal to Noise Ratio (SNR) and Tanh() operation.

The additional scaling by the SNR, and the Tanh() operation are done to perform maximum likelihood frequency estimation (MLE). This shows the ideal implementation (in a maximum likelihood sense) but to my understanding the SNR scaling is rarely done in practice (but could be if an SNR estimation is performed as part of the receiver operation), and the Tanh() operation is typically approximated with a piecewise linear operation to be the signal for the small signal case when the signal is less than ±1 or limited to the sign of the signal (+1/-1) otherwise.  So basically this path in the block diagram is hard limited for better performance with a limit that is dependent on the SNR of the received signal (a SNR dependent hard limiter is a typical operation within MLE implementations).

I believe that fred harris (fh) points out an advantage to this approach is how the derivative matched filter response is concentrated at the band edge where as we see through the convolution operation most of the information regarding frequency offset resides (when using a complex conjugate multiplier as a discriminator). The filter naturally blanks out the center portion of the spectrum which is not contributing further the our frequency error answer but would be contributing noise otherwise. Other advantages to this approach are that the symbol timing clock (timing recovery) can also be derived. I have not studied it personally in detail nor implemented this approach, but I believe fh makes the case that this implementation results in much lower pattern noise versus other carrier recovery (and timing recovery) approaches.


Minimum phase conversion with Cepstrum method, how to scale the result?

One possibility is to take the FFT magnitude of both, run a linear regression between the two magnitude vectors, and adjust the gain until the slope of the regression fit is 1.


What's the similarities and differences between Wigner transform and wavelet transform?

There is a nice tutorial: Time-Frequency Toolbox For Use with MATLAB tthat compares many transformations, notably continuous wavelet transform and Wigner and its avatars.
Basically, the CWT belongs to atomic linear representations (a decomposition onto wavelet atoms), while Wigner-Ville is an
energy distribution (to split the energy of signals over  time and frequency).
Since the tutorial comes with a toolbox and many examples, you can easily reproduce it on your signals.


If your data isn't white or for certain TED there is a lot of jitter.

If you're using a pre-amble, why not use a least-squares or grid search for timing detection and then heavily damp your NDA detector? Or don't even bother and just lock it down until the next burst.


Deciding on the FFT length for water waves in flume (wave facility)

So, the best news first:

Nyquist is your friend! As long as your sampling rate is sufficiently high, all the information contained in a band-limited analog signal can be represented by the digitized signal.

In other words: you can't sample too fast, and even if the ratio of sampling rate and signal frequency is irrational, you don't lose information.

So, you can follow the sample first, ask questions later philosophy¹.

Then, yes, if you want your periodic signal to occupy singular bins in your DFT, the DFT length needs to be an integer multiple of your period. And yes, FFTs can nowadays be implemented with a lot less restricted lengths². Your $10^N = 2^N\cdot5^N$ certainly is covered by the FFTW.

I'd like to encourage you to not only consider the DFT as frequency estimator.

MUSIC is able to calculate the spectrum at arbitrary frequencies (ie. it's an algorithm that you can ask "hey, what's the amplitude at f=0.123141414288?") instead of giving you an frequency-uniform amplitude vector, and if you know that you're looking for a fixed number of "peaks" in the spectrum, ESPRIT might be your thing. The accuracy of these algorithms is not bound by an orthogonal basis of frequencies, but on your ability to estimate the signals autocorrelation – which you can increase by adding more observation.

So, maybe you'd want to do the following:

Apply the DFT to your signal to get a "rough idea" of where your energy is in frequency
use MUSIC to generate a (pseudo-)spectrum around these frequencies of interest with a very fine resolution
As soon as you know how many tones to expect in the spectrum, throw ESPRIT at the problem to find exactly that number of dominant frequencies in your signal

¹: I'll need to make "SDR: Sample first, ask questions later" t-shirts.

²: The FFTW documentation (as far as I know) doesn't really specify, but in a half-sentence claims that things work best for lengths that can be factorized into prime factors $\le7$; one of the papers quoted on the homepage somewhere mentioned you can also have up to two factors of 11 and one factor of 13 in there to get an optimized FFT. That info is from the back of my head, might be wrong or outdated, and even if totally wrong: an absolutely naive DFT done by a matrix-vector multiplication will not put an unbearable load on your PC.
So simply use the DFT lengths you need, unless the vectors you need to transforms are in the millions of values.
And, directly from the fftw.org home page: Arbitrary-size transforms. (Sizes with small prime factors are best, but FFTW uses O(N log N) algorithms even for prime sizes.)
So, assuming you're doing a $10^6$ transform, and FFTW couldn't do better than $\mathcal O(N\log N)$, that'd be in the order of millions of floating point operations per transform. Assuming that your PC is about 20 years old and can do let's say 10 Million of these per second – do you have to worry at all?


How to normalize predesigned IIR filter coefficients to have unity gain?

In general, divide each of the $b$ coefficients by the maximum gain of your filter. This is equivalent to adding an attenuator before or after the original filter.

In this particular case, the maximum gain seems to occur at frequency 0 (where there is a double zero-pole cancellation).

At frequency 0, the gain of your filter is 36. So divide $b[]$ by 36.

Edit: the math

The frequency response comes from evaluating $H(z)$ at $z = e^{j\theta}$. Substituting and operating, and ignoring the phase, you get:

$|H(e^{j\theta})| = \frac{\sin(\theta/2)^2}{\sin(6\theta/2)^2}$

Taking the limit as $\theta \rightarrow 0$ you get $H(e^{j0}) = 36$.

Note that this is a simpler (and well known) filter applied twice. It corresponds the moving average of 6 samples (with the 1/6 factor missing).

So it is like two filters in cascade, each of which sums 6 consecutive samples. Each is a low-pass filter with gain 6 at 0 frequency.


Extracting different frequency bands (alpha, beta, gamma) from brain activity data

MATLAB's fft returns the FFT of the columns of a matrix X, so by transposing your input matrix before feeding it to the algorithm should work fine: F = fft(X.'); and then you will have the FFT for each channel of yours.

The FFT can be used to obtain the power spectrum. Here is a link with examples and explanations.

Once you have your power spectrum, then there are plenty of ways to extract different frequency bands (mostly based on filtering).


Design FIR filter with half zero coefficions

Oh, these are "Nyquist (M) filters", in your case M=2. You'll very often find them in mutlirate systems, especially as L-th-band-filters (i.e. filters that only let through 1/L of the Nyquist bandwidth) in decimators.
This is very fundamental for multirate systems, not really "research paper stuff", so your best bet is hence a good book on Multirate Systems. Fliege has written a classic "Multirate Digital Signal Processing", harris did, too, with "Multirate Signal Processing For Communication Systems" and Vaidyanathan seems to have written a popular book called "Multirate Systems and Filter Banks", but I haven't ever read that.


Please check this answer, which describes a few approaches to the same problem. Given that bird song is a monophonic signal (only one fundamental frequency at any point in time - as opposed to polyphonic) - and given that the timbre is irrelevant, the most interesting feature to extract for this classification task is a pitch contour.


Detecting presence of music in ambient sound

This is a tough problem. The music identification service Shazam actually has a database of acoustic fingerprints for songs. Music is composed of many variable acoustic components and there are a slew of algorithms that may be used to estimate several components. One simple example would be the zero crossing rate, which is a time domain algorithm. I would refer to this link on the subject for a brief overview of pitch detection algorithms and go from there.


Is the PSD of u1 and u2 computable ? Does it only depends on p1 and p2 ?

Yes. The Fourier Transform is a linear operation so the spectrum of the sum is simply the sum of the spectra. The PSD is just magnitude squared of the Fourrier Transform but since the frequency intervals are disjoint the squaring doesn't create any type of interaction between the two original spectral. So it's simply
$$P_{sum}(f) = \begin{cases}
 & p_1\text { if } f_{11} < f < f_{12} \\
 & p_2 \text{ if }  f_{21} < f < f_{22} \\
 & 0 \text{ otherwise }
\end{cases} $$


Sampling: question regarding sampling period

You sample at times $nT_s$, so the values of the sampled signal are

$$x_d[n]=x(nT_s)=\sin(10\pi nT_s)\big |_{T_s=0.1}=\sin(n\pi)=?$$

Can you take it from here?


Different amplitude for different time resolution with FFT

As you can see in MATLAB documentation for FFT
If you have Y = fft(x) then
$$Y_k = \sum_{j=1}^{N} x_j \cdot exp\left(-\frac{2 \pi i}{n} (k - 1)  (j - 1)\right)$$
$$x_j = \frac{1}{n} \sum_{k=1}^{N} X_k \cdot exp\left(+\frac{2 \pi i}{n} (k - 1)  (j - 1)\right)$$
Notice that the FFT is just a sum and the IFFT is divided by $n$, with this definition you know that  fft(ifft(x)) = x, but as the number of terms increase $Y_k$ will increase.
Your signal is a fixed duration pulse and you want a fixed amplitude in the frequency domain, so it is like an aperiodic Fourier transform
$$ \int x(t) exp\left( \frac{2\pi i}{n} u t \right) dt $$
A zero order approximation for this can be obtained by computing
A=20*log10(abs(fftshift(fft(AA*dt))));

For the pulse of width 1, max(A) will be 0dB this way.


Pole-zeros of a real-valued causal FIR system

Note the difference between the zeros at $0.3 \pi$ and at $0.8 \pi$.

The first one is clearly a zero crossing, much like $abs(x)$ at $x=0$.

At $\theta = 0.8 \pi$, however, the curve is tangent to the horizontal axis, much like $x^2$ at $x=0$. So you have a doulbe zero here.

So your zeros are:

2 zeros at $z = e^{\pm j 0.3 \pi}$
2 double zeros at $z = e^{\pm j 0.8 \pi}$


What is the third coordinate in "Projective Transformation" of an image and how to build the matrix of the inverse transform?

You have to understand that projective transformation is not linear.
It's $(x,y) -> (X(x, y)/Z(x,y) , Y(x, y)/Z(x,y))$, where X, Y, Z - linear functions.

$w$ and $w'$ are essential here - they help represent this non-linear transformation as linear - matrix multiplication - and that way this linear representation is a homomorphism (preserve "chaining" of operations).$w$ for input image is usually taken as $1$;
You get $w'$ from matrix multiplication ($w' = gx + hy + iw$ ) and the end result of transform is $(x'/w', y'/w')$

Generally speaking you can't just assign $i=1$, but due to homogeneity you can divide all the matrix by $i$ and get new matrix, representing the same transform with right-left element $1$

I would advise you to read some textbook on projective geometry for computer vision to better understand this subject. Those questions have some deep implications.


It may not be sufficient by itself, but since one of the problems has to do with lighting variations, a shadow removal pre-processing step may help.  The technique I'm thinking of is described in the following paper:

  "On the Removal of Shadows From Images", G. D. Finlayson, S. D.
  Hordley, C. Lu, and M. S. Drew, IEEE Pattern Analysis and Machine
  Intelligence (PAMI), Vol.28, No.1, Jan, 2006, pp.59-68.
  http://www.cs.sfu.ca/~mark/ftp/Pami06/pami06.pdf

The first part of the process produces an illumination invariant grayscale image, which is probably what you want in this case.  Here's a shot of the example they give in the paper:

(source: datageist.com)

In fact, taking the process one step further to produce a "chromaticity" image may be what you need to cleanly distinguish between the live and dead leaves.  Once again, an example from the paper:

(source: datageist.com)

The catch, however, is that the camera needs to be calibrated first.  Assuming that's possible, some combination of the representations they describe will probably make the other methods you're using more effective.


Carrier frequency offset in Bell 103 FSK system


3 Band Crossover with Linkwitz Riley Filter

A 2nd-order Linkwitz-Riley crossover needs 1st-order allpass phase compensation.

The Lowpass half of the LR crossover is:

$$\dfrac{1}{s^2+2s+1} $$

The Highpass half of the LR crossover is:

$$\dfrac{s^2}{s^2+2s+1} $$

With even-order LR crossovers, in order to sum-to-allpass one of the two sections has to have inverted polarity. Thus:

$$\dfrac{s^2}{s^2 + 2s +1} - \dfrac{1}{s^2 + 2s +1}= \dfrac{(s-1)(s+1)}{(s+1)(s+1)} = \dfrac{(s-1)}{(s+1)}$$


Why am I not getting a flat phase for a Gaussian pulse when doing a Fourier transform in Python?

The group delay will be flat assuming the filter is linear phase (which means the phase is linear vs frequency). Group delay is the derivative of phase with respect to frequency. Either unwrap your phase and detrend it or use the group delay function that is part of scipy.signal to evaluate directly. Any FIR filter (not just the Gaussian) will be linear phase if (and only if) the filter coefficients are symmetric or anti-symmetric.


attenuation of the central difference differentiator

Your filter forgets that in the discrete differences equation, you divide the difference in value by the difference in $x$; so, you'd need to multiply the result with ncycles*2*np.pi/nsamples to get the same.
Other than that, you're not compensating for the group delay of your FIR filter, either, so you want to shift the result by half the filter length, too.


You could implement some convolutions as recursive filters. Consider moving average as an example, to calculate the output signal, instead of summing over the averaging window, you could add the new input and subtract the part of input which exited the window from the output signal to obtain the output of next time step. Also some filters could be implemented as lattice.

But there is not a general method to implement all filters recursively or as lattice.


How to compare and measure the similarity between two between two signals?

The DTW approach is applicable to signals where there is an acceleration or deceleration between the observations of similar signals during the data capture and will show how similar or not the signals are besides the distortion from the "time-warping".

I am not sure that would be the best approach for this case given nothing is indicating there will be a change in the rate of time between the sensors (although to a second order the sampling clock in one will not be exactly matched to the sampling clock in the other so DTW could be useful if the sampling clock offset exceeded the coherence bandwidth of the capture, which is $1/T$ where $T$ is the length of the capture for which we will use correlation, which I will describe next).

This looks like an ideal application for the cross correlation function, which will show the correlation between the two waveforms for every time offset between the two. This is done by first removing the mean from each waveform, and then multiplying the two resulting zero-mean waveforms together element by element and summing the result, repeating for each possible sample shift between the waveforms. These results can be scaled by the product of the standard deviation of the two waveforms, which would normalize the result similar to what is done in the Pearson correlation coefficient.

Another similar approach is to use the Wiener-Hopf equations to solve for the equivalent "channel" between the two waveforms, which I have explained in more detail at this post which evaluated the signal received by two microphones:

Compensating Loudspeaker frequency response in an audio signal


frequency response of a signal in which two signal added in

What you ask for should rather be called a frequency spectrum or Fourier transform of the resulting signal. The term frequency response is reserved for linear time invariant systems.

Anyway, the answer to your question can be found using the z-transform of the signals involved. First let's expand the two source signals into a series of unit delays:

$$S_1(z) = \sum_{k\in\mathbb{Z}} a_k z^k$$
$$S_2(z) = \sum_{k\in\mathbb{Z}} b_k z^k$$

Your combined signal can be written as

$$F(z) = \sum_{k\in\mathbb{Z}} a_k z^{2k} + \sum_{k\in\mathbb{Z}} b_k z^{2k+1}$$

and we can plug in $z=\exp(i\omega)$ to get the Fourier transform:

$$F(\omega) = \sum_{k\in\mathbb{Z}} a_k \exp(ki2\omega)+\sum_{k\in\mathbb{Z}} b_k \exp(ki2\omega)\exp(i\omega)$$

Using the Fourier transforms from the z-transforms above we can see that this actually equals:

$$F(\omega) = S_1(2\omega) + S_2(2\omega)\exp(i\omega)$$

You can interpret this result as the additive combination of the two original signal spectra, each scaled in frequency direction by $1/2$ and one phase shifted proportionally to $\omega$.

Since your original signals are of finite length, your initial spectra are sinc-shaped, centered at the same frequency and with a constant phase shift. I'll leave working out the details to you.


What you did by saying a $0$ bit is mapped to a voltage of zero Volts, and a $1$ bit is mapped to $5$ Volts is exactly what it means to map bits to symbols. Your symbols are 'zero Volts' and '$5$ Volts'. A bit is just a unit of information, and you have to decide which type of analog signal you choose to transmit the source information. You always have to map information bits to some analog signal in order to transmit it over an analog channel.
There are many ways to choose appropriate symbols. Your suggested mapping from bits to symbols requires a baseband channel that can also transmit DC, because your transmitted signal will have a DC value of about $2.5$ Volts, assuming that $0$s and $1$s are equally likely. There are baseband channels that are AC-coupled, i.e., you can use them to transmit very low frequencies, but not DC. In that case you could choose a positive voltage for transmitting a $1$, and a negative voltage of equal magnitude for transmitting a $0$ ("binary antipodal signaling"). There are many standard ways to encode bits for transmitting them over a baseband channel. Here you can start reading up on these so-called line codes.
Note that you can also map more than one bit to a symbol. E.g., you could choose voltages $-3$, $-1$, $1$, and $3$ to encode all possible two-bit sequences. This makes your symbol rate (baud rate) half the bit rate, which requires less channel bandwidth. There's a trade-off between bandwidth and noise immunity. For a fixed transmit power, the symbols are now closer to each other than before, so it is more likely to confuse one symbol for another when there's additive noise on the channel.
There are also channels that have a bandpass characteristic (such as a radio channel). In that case you need to use carrier modulation, which shifts the spectrum of the transmitted signal to the desired frequency band. You can use two orthogonal carriers at the same frequency, such as $\sin(2\pi f_c t)$ and $\cos(2\pi f_ct)$ to transmit two symbols at the same time (quadrature modulation). In that case, symbols are represented by complex numbers (because you have two symbols at a time) and you get a two-dimensional symbol constellation.
It is common to use constellations with many more than just two symbols, e.g., $256$-QAM, which uses one symbol to encode $8$ bits.
All the above gives a very basic and simplistic view of how encoding and modulation works. If you want to dig a bit deeper you should get a good textbook on digital communications.


How to get the MLSEEqualizer work with IFFT and FFT blocks in the communication system?

Here is I modified it, and it's ok now:

% Modified Code: IFFT and FFT blocks are used but the system shows worse BER results.
%   Equalize a QPSK signal transmitted through a dispersive channel using MLSE

qpskMod = comm.QPSKModulator(0,'SymbolMapping','Binary');
qpskDemod = comm.QPSKDemodulator(0,'SymbolMapping','Binary');
% Channel coefficients
chCoeffs = [.986; .845; .237; .12345+.31i];
N=512; % Modulated signal length
mleq = comm.MLSEEqualizer('TracebackDepth',10,...
                  'Channel',chCoeffs, 'Constellation',[1 1i -1 -1i]);
% Create an error rate calculator
ber = comm.ErrorRate;
for n = 1:10
  data= randi([0 3],N,1);

 % Modulate the data and convert it to time domain.
  modSignalx = step(qpskMod, data);
  modSignal=sqrt(N)*ifft(modSignalx);

  % Introduce channel distortion
  chanOutput = filter(chCoeffs,1,modSignal);

  % Equalize the channel output and demodulate
  eqSignalx = step(mleq, chanOutput);
  eqSignal=(1/sqrt(N))*fft(eqSignalx);

  demodData = step(qpskDemod, eqSignal);
  % Compute BER
  a = biterr(data(n), demodData(n));

b(n)=a(1);
end


Non Linear Remapping of the Pixel Values

Yes, using a non-linear mapping function on the intensity. Let us for the moment suppose you have a grayscale image (because of the words "black and white", "shadings"), with luminance $0\leq p_{i,j} <1$.

The piecewise linear mapping: $f( p_{i,j}) = \frac{p_{i,j}}{2}$ if $p_{i,j}<\frac{1}{2}$, and  $f( p_{i,j}) = \frac{3}{2} (p_{i,j}-\frac{1}{2})+\frac{1}{4}$ if $p_{i,j}\geq\frac{1}{2}$ can do the job.

This is a simple example of "Non-linear Contrast Enhancement". If that works for you, you can try with smoother functions.


Does Short-time Fourier transform impact quality of signal?

STFT (short-time Fourier Transform) should not impact the quality of the reconstructed signal if properly implemented.
Furthermore, any signal $x[n]$ should be exactly recovered from its STFT analysis coefficients by a proper implementation of the associated inverse STFT;i.e., The forward and inverse pairs should exactly undo each other.
One should keep in mind about the overlap parameter and window functions however;i.e., if the window includes zeros, the inverse may not be possible at some overlaps. If the window does not contain any zeros, then the inverse is always computed.


If the mean of a random process is constant, does it imply the process is first order stationary?

In the usual sense of the term, first-order stationarity means that the first-order distribution of all the random variables is the same: each $X_t$ has the same CDF, and so the same pdf (or pmf) too if the random variables are continuous (or discrete). If the random variables have a mean, then they all have the same mean. But, a sequence of independent
identically distributed Cauchy random variables is a first-order stationary process (as well as a strictly stationary process for that matter) but it cannot be said to have constant mean because none of the variables have a mean at all.

Since your response to my query on the problem statement is that your definition of first-order stationarity is the  standard one, then the answer to the question in the title of your question, viz.

  If the mean of a random process is constant, does it imply that the
  process is first-order stationary?

is NO, and an example of a process that has constant mean but is not
first-order stationary is a sequence of constant-mean random variables whose variances depend on time.

The answer to the question in the text, viz.

  Is a WSS process necessarily first-order stationary?

ia also NO.  Consider the random process
$\{X(t)\colon X(t)= \cos (t + \Theta), -\infty < t < \infty\}$
where $\Theta$ takes on four equally likely values $0, \pi/2, \pi$ and $3\pi/2$ and so $X(t)$ takes on four equally likely values
$\cos(t)$, $\cos(t+\frac{\pi}{2}) = -\sin(t)$,
$\cos(t+\pi) = -\cos(t)$, and
$\cos(t+\frac{3\pi}{2}) = \sin(t)$.
  It is easy to see that in general $X(t)$ and $X(s)$ have different distributions and so the process is not first-order stationary. On the other hand, $E[X(t)] = 0$ for all $t$ while
\begin{align}
E[X(t)X(s)]&= \left.\left.\frac 14\right[2\cos(t)\cos(s) + 2\cos(t+\pi/2)\cos(s+\pi/2)\right]\\
&= \left.\left.\frac 12\right[\cos(t)\cos(s) + \sin(t)\sin(s)\right]\\
&= \frac 12 \cos(t-s)
\end{align}
and so the process is wide sense stationary, but its first-order distribution does depend the choice of $t$ and so the process is not
first-order stationary.

Note: the above proof uses a very commonly accepted definition of wide sense stationarity as found in numerous books, online resources, Wikipedia etc, viz. a random process $\{X(t)\}$ is called a wide sense stationary (WSS) process if

$E[X(t)]$ is a constant, and
$R_{X}(t,s) = E[X(t)X(s)]$ depends only on the difference $t-s$ of the two arguments and not on the individual arguments $t$ and $s$.


The easiest way to understand the concept of filter order is to look at the Laplace (continuous time) or Z-transform (discrete time) of the filter.  The filter order is the order of the transform's denominator.

For example, the Laplace transform of an integrator is $\frac{1}{s}$.  Extending the time of the integration does not change that it is an integrator, it just changes the bounds of the integration.  If you were to make it a double integration, though (two integrators in series), then the Laplace transform would be $\frac{1}{s^2}$ and you would have a 2nd order filter.

Averaging can be considered a partial integration with a gain factor, so your integrator with averaging would be a variant of a double integrator, making it a second order filter.


Exponential decaying Swept / chirp signal

Frequency is the derivative of phase, so you need to start with the expression for frequency, which in your case is
$f=200\exp(-\alpha*t)$

At time $t=0$ this will have the frequency of 200, but we need to solve for $\alpha$ so that the frequency at time $t=2$ is 20.  Solving for $\alpha$ gives $\alpha=1.1513$.
To get the correct expression for the phase you need to integrate the above expression for the phase i.e. solve

$f =\frac{1}{2\pi}\frac{d\phi}{dt}$

The following matlab code generates this signal.

fs=500; %sampling frequency
t=(0:2*fs)/fs; %Time vector
alpha=1.1513;
sig = exp(-j*2*pi*200/alpha*exp(-alpha*t));

Note that just because the instantaneous frequency behaves this way doesn't mean the spectrogram will follow it exactly. It will get better if you have a higher time-bandwidth product.


Shape of structuring elements for morphological gradients

I looked up in Jaehne, Gonzalez, Soille (the one you've posted as well as Mathematical Morphology and Its Applications to Image and Signal Processing) and some other special morphological papers and haven't found neither any design criteria for the structuring element nor any special hints why it has to be symmetrical.

Personally I think that a symmetrical SE is good for the same symmetrical effects on the object you want to modify.
With my existing experience I would not use a non-symmetrical SE, because I can't take it for any object or scenario and I don't know how it will react for other cases.

Nevertheless it's an interesting question and I am trying to get an answer.


How to segment femoral and tibial knee cartilage individualy?

Looks like you have to do some motion tracking of points to separate the two. Is it possible to take images before and after some motion? That way you could model how the cartilage moved as a function of tibia or femoral motion, and you might be able to separate what portions of the cartilage belong to the tibia and femur.


You just triggered one  old itch, so this is a partial (both meanings) answer only. For discrete-time, uniformly sampled STFT, one often uses the usual fixed-window with $3/4$, $1/2$, $1/4$ overlaps, and somehow uses  direct closed form inverses. However, as long as there is redundancy (and completeness), there exists a infinite option for inverses.

In a   collaborative work (with Jean-Christophe Pesquet and Jérôme Gauthier), following other works on vector-frames, we used the formalism of complex oversampled multirate filter banks

to model analysis STFT, with the aims of

checking conditions under which the STFT were invertible, given a hop between time-frames,
building "better" inverses (better localized in time or frequency) with optimization.

Our aim at that time was double:

for data in higher dimensions $d$, reduce the STFT redundancy $r^d$ as far as we could,
-for 1D data, reduce the number of channels and increase the hop as we could as similar quality for real-time event detection.

All was based on polyphase matices, as we did consider windowed transform only as a special case;  so the main outcomes were the following:

given any redundant analysis filter bank (complex or not, windowed-DFT like or not), check whether it is invertible
in most cases, anything redundant is almost surely inversible,
windows with negative values can be used, liked Kaiser windows, see On the Optimization of Oversampled DFT Filter Banks, 2007

you can reduce the degrees of freedom in the oversampled synthesis to perform better concentrated inverses after selection/denoising for (close-to) real-time computations,
ensure the Hermitian character in the inverse for scalar threshold.

There were remaining issues, among which: the difficulty to go back to the orthogonal setting when the redundancy tends to one. Weirds behaviors were observed, possibly related to Farey sequences, but I have no definitive clue. And
References can be found in Optimization of Synthesis Oversampled Complex Filter Banks, 2009. Some related codes is available at SURE-LET Optimal oversampled Complex Filter Banks synthesis toolbox. We do plan a better release, but it has been remaining upcoming for months.

However, we did not investigate the non-COLA effect, and the necessary condition is not stated. My belief is that it is not necessary (due to redundancy), but this should be proven.


Weird results when using MATLAB's $\tt butter$ function?

The answer to your question can mainly be found in the Matlab documentation. The butter function can design filters either in the analog domain as well as in the discrete domain. In order to make an analog filter, you can do so by giving the string 's' as a third parameters. The transfer function can be defined by the zeros/poles/gain or by numerator/denominator coeefficients. The two available syntaxes are

[z,p,k] = butter(n,Wn,'s')
[b,a] = butter(n,Wn,'s')

When designing a filter in the discrete domain, the normalized frequency is a ratio of the cutoff frequency over the biggest frequency that you can represent according to Nyquist theorem. For instance : if you have a sample time of 1ms, hence a sampling frequency of 1Khz, the biggest frequency you can represent is 500Hz.  Designing a filter with a normalized frequency of 0.5 will put the cutoff frequency at 250Hz, since 250/500 = 0.5

In your question, you specify that the cutoff frequency didn't seems to match the -3dB point. Can you give us more detail about your issue. Maybe show us the work you have done until now ?


The question has an answer on Math SE. The argument goes as follows:

\begin{eqnarray}
\int_{-\infty}^\infty \left| \frac{\sin t}{\pi t} \right|dt &=& \frac{1}{\pi}\sum_{n=-\infty}^\infty \int_{n\pi}^{(n+1)\pi} \left| \frac{\sin t}{t} \right| dt \\
&>& \frac{1}{\pi} \sum_{n=-\infty}^\infty \frac{1}{|n+1|\pi}\int_{n\pi}^{(n+1)\pi} |\sin t|dt \\
&=& \frac{2}{\pi^2} \sum_{n=-\infty}^\infty \frac{1}{|n|}\\
&\rightarrow& \infty
\end{eqnarray}
where in each interval $[n\pi, (n+1)\pi]$ we lower bound $1/t$ with $1/(n+1) \pi$.


A system having impulse response $ h(t)=u(t) $ stable or not?

Your derivation is correct. The impulse response is not absolutely integrable, hence the system is not BIBO stable. From the corresponding transfer function $H(s)=1/s$, you can see that there is a single pole at the origin. Systems with single poles on the imaginary axis, like the integrator in your example, are also called marginally stable. These systems are not BIBO stable, transients do not decay, and certain inputs can cause unbounded outputs.


harmonic algorithm for extending compressed bandwidth of audio file

there is also this answer in this comp.dsp thread.

high-pass filter your input signal (keep the unmolested input around for later).  the cutoff frequency should be adjustable.  i would think a 2nd-order HPF with a Q approx= 0.707 (Butterworth) would be good enough.
square that HPF output,
add (maybe with an adjustable gain coefficient) that squared HPF back to the original input and output that audio sample.

see if that brightens up your day.


Is an FIR filters 'transition width' frequency dependent?

You're right: the required filter order is approximately inversely proportional to the desired transition bandwidth $\Delta\omega$, regardless of the cut-off frequency. This is reflected in the empirical formulas for estimating the required filter orders for the Kaiser window design method as well as for the Parks McClellan equiripple design of low pass filters:

$$\begin{align}M&=\frac{-20\log_{10}\delta-8}{2.285\Delta\omega}&\qquad\text{(Kaiser Window)}\tag{7.63}\\{}\\
M&=\frac{-10\log_{10}(\delta_1\delta_2)-13}{2.324\Delta\omega}&\qquad\text{(Parks McClellan)}\tag{7.104}\end{align}$$

where $\delta$ is the maximum allowable approximation error (Kaiser), and $\delta_1$ and $\delta_2$ are the maximum pass band and stop band errors, respectively (Parks McClellan).

The equation numbers above are taken from the book Discrete-Time Signal Processing by Oppenheim, Schafer, Buck, 2nd edition.


They are both valid:

$10 \cdot \log_{10}$ is for power ratios
$20 \cdot \log_{10}$ is for voltage/pressure/etc ratios (these are called "field quantities" or "root-power quantities")

See Field, power, and root-power quantities and Decibel definition

Digital values are just unitless numbers, so you need to think about what they ultimately represent.  Typically they're root-power quantities (like audio pressure waveforms), but could conceivably represent power waveforms.


Radar signal processing - Difference between range gate/range bin

If we wanted to add a nuance, the "gate" is the actual switching on of the receiver for the duration of observing the reflected signal (which can be as narrow as a single sample, which even a single sample for analog to digital conversions is an integrated observation over a sample time), while the range bin is which "gate" you are in. (Imagine a receiver that can have several gates covering the full possible range of the target). This originated from pulsed radar systems which is what my description above would apply to but can be shown to be equally applicable to other radar methods such as FMCW where an ambiguity in position can occur based on the repetition rate of the waveform.

For tracking, I believe this is where the term "Early-Late Gate" originated: In this use the gate is subdivided into an early half and late half (or can be described as separate gates, one slightly earlier and the other slightly later in vicinity of known target location), and the central portion of the gate position is adjusted earlier or later based on the difference between the early and late gate.  Here too we are "gating" the received signal to limit our observation to a narrow time interval (or even one sample).


What's the formula for choosing FFT sample rate and window length to make sure the frequency I want is in the center of a bin?

The concept you need is to have a whole number of cycles within your sample frame.  The cycle count will be the bin index.

For the math, consider carrying units:

$$ \frac{Cycles}{Frame} = \frac{Cycles}{Second} \cdot \frac{Seconds}{Frame} $$

Where:

$$ \frac{Seconds}{Frame} = \frac{  \frac{Samples}{Frame} }{ \frac{Samples}{Second} } $$

Solving for your sampling Rate:

$$ \frac{Samples}{Second} = \frac{  \frac{Samples}{Frame} }{ \frac{Seconds}{Frame} } $$

$$ \frac{Samples}{Second} =  \frac{ \frac{Samples}{Frame} \cdot \frac{Cycles}{Second}  }{ \frac{Cycles}{Frame} }  $$

If you have a pure tone, there isn't a need to center it on a bin to find your parameters, although it is much easier to do it that way.

To find the frequency, you can use the formulas in these blog articles:

Exact Frequency Formula for a Pure Real Tone in a DFT
Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT
Improved Three Bin Exact Frequency Formula for a Pure Real Tone in a DFT

Once you have the frequency, you can calculate the amplitude and phase with the formulas in this article:

Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1

You should also realize that you don't need to call a full FFT.  You simply need to calculate the two or three bin values you are interested in which is generally a lot quicker.


How does minimum-latency partitioned convolution reverb work when you receive input samples in chunks, rather than one at a time?

This is typically done using a segmented overlap add method or sometimes also refered to as a block convolver.

Let's assume your block size if 512 (makes the numbers a little easier).

Chop up your impulse response into 32 blocks of 512 samples each. Zero pad each block to 1024 samples and FFT. You know have 32 filters $H_0(z) ... H_{31}(z)$
On each new input block: zero pad to 1024 samples and FFT. Keep the most recent 32 signal spectra around, so you have $X_N(z), X_{N-1}, ... X_{N-31}(z)$
Multiply the spectra and sum them, $Y(z) = \sum_{k=0}^{31} H_k(z) \cdot X_{N-k}(z)$
Inverse FFT the result: you get 1024 time samples. The last 512 are the "overlap". Store for the next block.
Take the first 512 samples of your inverse FFT and add the overlap from the previous block. That's your output.

Repeat steps 2-5 for each new block.


Contradiction between complex baseband and real-valued baseband

As a sidenote:

When you demodulate your signal $x_{\mathrm{rf}}(t)$  you say that there is some phase difference $\phi$. This translates to crosstalk (i.e., the quadrature and in-phase components influence each other at the receiver). It's easy to see that if $\phi = 0$ (namely, we know exactly when to sample the received signal), there is no crosstalk at all:

$$y_i = \sum_{k=0}^K \left( g_k x_i[n-k] + g_k \tan(0) x_q [n-k] \right) =  \sum_{k=0}^K  g_k x_i[n-k]$$

In the baseband domain, this means that the filter $c_k$ is purely real (if it was purely imaginary, then there wouldn't be crosstalk either, you just have to invert the channels to retrieve the information correctly).

The answer:

Whether there is crosstalk or not, both of your equations get to the same result, but with different names. Using your notation, $a_k$ in the first equation would be equivalent to $g_k$ in the second one, and $b_k$ to $-g_k\tan(\phi)$. I believe that your confusion arose from thinking of $a_k$ and $b_k$ as two different filters (which is correct), instead of thinking of them as the real and imaginary part of a single filter ($c_k$).


How to find a maximal possible value of a cross correlation

Yes, but only if you know the mean and variance (or standard deviation) of the signal. The maximum correlation will be the (variance of the signal + the mean-squared) times the total number of samples.   If you do not already know the variance, calculating it would be similar to calculating the auto-correlation in terms of processing so nothing is really saved in that case. But in some cases you may have prior knowledge or ability to more easily estimate the variance (from a smaller sample set for example if the process is white and stationary).

To check my logic, consider what is involved to calculate correlation for a sampled signal when it is aligned with itself (for the maximum autocorrelation) compared to calculating the variance for the signal:

To correlate two signals: You multiply the two signals sample by sample, and then sum the result.

The equation to calculate the variance is similar, but you must subtract the mean from each sample before doing the product, and then you must divide by the total number of samples. (Important to note, this is not the unbiased estimate of the variance which scales by N-1, but the calculation for the 2nd Moment about the mean, which scales by N). The variance therefore is just a normalized autocorrelation at lag=0 with the mean removed.

To put very simply, CorrMax is just 1/N * the second moment. (The second moment is the average of the squares or mean-squared value). Variance is the 2nd moment with the mean removed (about the mean), so we have to add it back in as we've done above if we were starting with that particular factor known. However consider if we already know the mean-squared value, (or the square of the rms value since rms is commonly used), then the relationship is the mean-squared value times the number of samples.

Therefore in summary:
To determine the maximum correlation for a sequence X that is N samples long,  assuming you already have either the mean and variance (2nd moment about the mean), use:

$$ CorrMax(X) = N(\sigma_X+\mu_X^2)$$

Where:

X is a vector of length N

$ \sigma_X $ is the second moment of X about the mean

$ \mu_X $ is the mean of X

Or alternatively calculated from the mean-squared value for the sequence,

$$ CorrMax(X) = N(X_{rms}^2)$$

Where:

X is a vector of length N

$ X_{rms}^2 $ is the second moment of X (mean-square of X)


Understanding the Difference Between MAP Estimation and ML Estimation

You have a set of message set $m_i$, $0 \le i \le N-1$. (For example, QPSK will be $N=4$). For the transmitted message $m_i$, the corresponding symbol vector is $\textbf{x}_i$, and the received symbol vector is $\textbf{y} = \textbf{x} + \textbf{w}$, where $\textbf{w}$ is the AWGN at the receiver. The above is a simplified baseband model assuming a simple Line-Of-Sight(LOS) channel without any delay.

At the receiver, after observing $\textbf{y}$ you arrive at a decision for the transmitted symbols, $\tilde{\textbf{x}}_i$. The decision you make is such that the probability of error $P_e = P(\tilde{m_i} \ne m_i)$ is minimum. In other words, the Probability of being correct $P_c = 1 - P_e$ has to be maximized. Using a hypothetical rule you have made a decision $\tilde{m}_i$ based on $\textbf{y}$, so
$$
P_c(\tilde{m} = m_i;\textbf{y}) = P(\tilde{\textbf{x}} = \textbf{x}_i;\textbf{y}) =
P(\textbf{x}_i|\textbf{y})P(\textbf{y})
$$
Here I have made use of formula $P(AB) = P(A|B)P(B)$. Event A is the event that $x_i$ was transmitted. Our decision is correct if $\tilde{x_i}$ is interpreted as $x_i$. We are trying to maximize the probability of getting this decision correct. In order to maximize the above term among all message set $0 \le i \le N-1$, the term $P(\textbf{y})$ can be ignored since it does not depend on $i$. Hence

$$
\tilde{m_i} = argmax_i\,\, P(\textbf{x}_i|\textbf{y})
$$
You are deciding the $i$ which maximizes the above probability after observing $\textbf{y}$. Hence the whole method above is called Maximum A-posteriori Probability decision because you are maximizing a-posteriori probability $P(\textbf{x}_i|\textbf{y})P(\textbf{y})$

The term above can be re-written as $P(\textbf{x}_i|\textbf{y}) = P(\textbf{y}|\textbf{x}_i)P(\textbf{x}_i)/P(\textbf{y})$.

If we assume equal probability for all messages in the set $\textbf{m}$ then we can ignore the term $P(\textbf{x}_i)$ since they are equal for all values of $i$. Hence
$$
\tilde{m_i} = argmax_i\,\, P(\textbf{y}|\textbf{x}_i)
$$
This is the Maximum Likelihood detection rule. We are maximizing the likelihood probability $P(\textbf{y}|\textbf{x}_i)$.


What interesting audio can be synthesized with ifft

You can take a look at Photosounder. The idea behind is simple: Take an image as if it were an spectrogram and compute its IFFT, the result is  always interesting at least! It tends to work better with abstract images and geometric shapes, a great example is the video on the upper-right corner, "The sound of fractals".

I would never try to replicate a real instrument using this method, but you can encourage your fellows to try it as a challenge...


Relation between time domain, DTFT domain and frequency domain

If this is a homework problem, shield your eyes.

Otherwise, it is a simple units issue.

S has units of 1k samples per second.

Your frequency is expressed in radians per sample.

There are $ 2 \pi $ radians per cycle.

So the puzzle now is to find the relationship between $\frac{samples}{second}$, $\frac{radians}{sample}$, $\frac{radians}{cycle}$, and your goal of $\frac{cycles}{second}$,


Corresponding frequency domain filter of spatial domain filter in image processing?

So, the $M\times N$-DFT is very well-defined to be (aside from the usual ambiguity about the factor up front)

$$
X[k,l] = \frac{1}{\sqrt{MN}}\sum_{n=0}^{N-1} \left(\sum_{m=0}^{M-1}x[m,n]\cdot e^{-jmk\frac{2\pi}{M}}\right)e^{-jnl\frac{2\pi}{N}}
$$

In your $M=N=3$ case, that becomes

$$
X[k,l] = \frac{1}{3}\sum_{n=0}^{2} \left(\sum_{m=0}^{2}x[m,n]\cdot e^{-jmk\frac{2\pi}{3}}\right)e^{-jnl\frac{2\pi}{3}} \tag1
$$

Furthermore,

$$ x = \begin{bmatrix}
0 & 1& 0 \\
1& 0 & 1\\
0 & 1&0 \\
\end{bmatrix} $$

so that all there's a lot of zeros.

Simply expand the sum from $(1)$, filling in the values of $x$.

$$\begin{array}
\,3X[k,l]&=& &\left(e^{-jk\frac{2\pi}3}\right)e^{j0} & \text{first row}\\
&& +& \left(e^{-j0}+e^{-jk\frac{4\pi}3}\right)e^{-jl\frac{2\pi}3}& \text{second row}\\
&&+ &\left(e^{-jk\frac{2\pi}3}\right)e^{-jl\frac{4\pi}3}&\text{third row}\\
&=&&e^{-jk\frac{2\pi}3}\\
&&+&e^{-jl\frac{2\pi}3} + e^{-j\frac{2\pi}3\left(l-2k\right)}\\
&&+&e^{-j\frac{2\pi}3\left(k-2l\right)}
\end{array}$$

So,

$$X[k,l]=\frac13\left(e^{-jk\frac{2\pi}3}+e^{-jl\frac{2\pi}3} + e^{-j\frac{2\pi}3\left(l-2k\right)}+e^{-j\frac{2\pi}3\left(k-2l\right)}\right)$$

Set in all values for $k$ and $l$, and get your $3\times3$ DFT output .


What are F0 counters in Speech signals? What is represented by F0?

Check this source

  To begin with, let’s remember what the fundamental frequency is and in
  what tasks it may be needed. The fundamental frequency, which is also
  referred to as F0, is the vibration frequency of the ligaments when
  pronouncing voiced sounds. When pronouncing unvoiced sounds, for
  example, by whispering or uttering hissing and whistling sounds, the
  ligaments do not vibrate, which means that this characteristic is not
  relevant.
  *Note that the split on voiced/unvoiced sounds is not equivalent to the split on vowels and consonants.

  The variability of the fundamental frequency is quite large, and it
  can vary greatly not only between people (for lower average male
  voices the frequency is 70–200 Hz, and for women it can reach 400 Hz),
  but also for one person, especially in emotional speech.

  F0 is used in a wide range of solutions:

Emotion Recognition
Sex determination (male/female voices)
Speaker deterioration, or splitting the speech into phrases
In healthcare, detection of the pathological characteristics of the voice (for example, using the acoustic parameters of Jitter and

  Shimmer): F0 could be used for the detection of signs of Parkinson’s
  disease 1; Jitter and Shimmer can also be used for Emotion
  Recognition


Confusion about subtle difference between discrete-time and continuous-time

The problem is that $$e^{j 2 \pi t} \neq 1,$$ except in the case when $t$ is an integer. However, in that case, you have the discrete-time formula: $n$ is assumed to be an integer. In most DSP literature, wherever you see $n$ used as in this case, you can safely assume it's an integer.


What you want is basically a delay
$$x(t-\tau) \Rightarrow e^{-j\omega \tau}$$
If you can live with $k$ being quantized you can simply use a single tap delay FIR filter. If you need more granularity you need to implement a fractional dealy. I would cascade a shot fractional FIR filter with a long single tab bulk delay filter.
If $k$ is positive the filter becomes a non-causal. The best way to deal with this depends a bit on your application. If you do off-line processing, you just need to delay all other signals paths by the same amount.


A quick scribble on paper shows that FMCW radar has a resolution of

$$
\Delta d = \frac{c}{2b}
$$ with $c$ being the speed of wave propagation (aka speed of light, essentially), $b$ being the bandwidth, and $2$ coming into play because the wave has to trave both ways. For your $\Delta d\overset != 1\mathrm m$,

$$
1\mathrm{m} = \frac{3 \cdot 10^8 \frac{\mathrm{m}}{\mathrm{s}}}{2b}
$$

which means you'd only need $b\ge 150\, \mathrm{MHz}$ in my calculation.

That's something that software defined radio devices can do nowadays.

Hence, to answer your questions:
My questions:

  How, if at all, high-resolution, high-bandwidth digital radars are realised in practice?

Simply be mixing the signal down to baseband (or a low IF), and then sampling it with an ADC.

  what would make sense in case of FMCW radar is to mix the transmitted signal and the target return signal and digitize only the beat frequency for Inverse FFT and further analysis.

Yes, that's the way it's done in most "small" FMCW radars, for example automotive radars.

  In that case, could you still call it a digital or software defined radar?

Puh. Digital: yes.
Software defined: well, if everyone asked me, no, but paper doesn't blush, so there's probably someone out there who calls this kind of radar "software defined". Please don't do that.


About output signal digital band pass filter

What you are seeing is not "noise", but the step response of your filter.

Note that the filter is started with a null state (unless you explicitly call filter with state variables).

Your signal starts at value 10, so you will get something very similar to 10 times the filter's step response.

If you want to eliminate this transient, then you can either:

Subtract 10 from your input signal, so that it starts at 0.
Prepend your signal with enough samples so that the transient does not overlap with your signal of interest. These samples can be all 10 in value, for example.
Multiply your whole signal by a window function; something like yw =
y .* hann(length(y))'. This way the whole signal is smoothed at the edges and you minimize the transient.


What are the advantages of using fft2 over fft

They perform two different mathematical operations

FFT executes a 1-dimensional Discrete Fourier Transform one each column of the input matrix (or the first non-singleton dimension)
FFT2 executes a full 2-dimensional Discrete Fourier Transfrom on the entire matrix.

Which one you want to use depends on your specific application. One is not inherently "better" than the other, they are just two different things. In particular

$$\mathcal{F}_{2D}(X) = \mathcal{F}(\mathcal{F}(X)')'$$


How to determine if a filter is bandpass/stopband from its pole-zero diagram in z-domain

For strongly selective filters that have well defined pass and stop bands and narrow transition bands and which, therefore, can approximate ideal brickwall filters to some desired degree, poles and zeros are generally distributed around and closer to the unit circle. Normally you would like to have all of them inside the unit circle but due to other reasons, such as linear phase FIR impulse response, some zeros can be put outside of unit circle too.

Effective frequency of a single pole/zero is the angle of its location from real axis counter clockwise. If that's a zero, it will stop frequencies nearby. If it's a pole it will amplify freqencies nearby.

By distributing those poles and zeros in groups you will be defining pass bands (that include band of poles) and stop bands (that include band of zeros), or bandpass filters as mixture.

When the magnitude of a single zero is away from unity, its selectivity decreases;i.e., pass-stop band distinction dissapears. When the magnitude is larger than one, it will also have an overall gain despite being a zero. If the magnitude is close to zero (zero at origin) then the gain approaches one.

Similarly, when the magnitude of a single pole is away from unity, its selectivity decreases;i.e., pass-stop band distinction dissapears. When the magnitude is larger than one, it will also have an overall attenuation despite being a pole. If the magnitude is close to zero (zero at origin) then the attenuation approaches one.


How should I pre-process a real valued signal in order to use Kay's estimator?

The tool for converting real signals into their analytic representation is the Hilbert transform.
Suppose your signal were a projection of some helical rotation with variable amplitude onto the real-time plane, like in the image below.

Source
The Hilbert transform produces such complex signal given its real part. It is a linear transform and it's very easy to do in frequency domain. Without going into its mathematics and derivations too much in depth, the Fourier Transform imaginary part of your signal is the same as the one of your real signal multiplied by $-j$ (rotated by 90 degrees). By symmetrical properties of real signal, you get the following relationship:

All your negative ferquency components become 0.
Your DC component stays the same.
All your positive frequency components double

In Matlab, for example, you would do the following:
a = rand(1,201);

hilbert_a = ifft( [ 1, 2*ones(1,100), zeros(1,100)] .* fft(a) );

or simply use built-in the hilbert function.


Will decorrelating an audio signal help reduce interference patterns at low frequencies?

You will simply change the interference pattern. By modifying the phase at different frequencies you will simply change the phase relationship between the two speakers and thus move the locations of any pressure nulls/nodes at that frequency. It would be pure luck as to whether this would 'reduce' the interference pattern; also depending on your definition of the reduce in this context.

With a single loudspeaker you can generally control the sound pressure at one location. Two speakers would give you control at two discrete points. N speakers, at N points. In a reverberant environment you will also need to identify the room response as well to achieve this.

The above rule can be broken, by using some kind of optimisation (least-squares as probably the simplest) to optimise the freqeuency response at N control locations toward a certain target (flat frequency response?), using M speakers, where M < N. You wouldn't get very far with 2 speakers though.

Read up on sound zones for some interesting background material on this and some rather extreme speaker configurations :)


What is the name of this color transformation?

I do not know the name of your transformation $T$, but it can't be the Karhunen-Loeve transform: The Karhunen-Loeve transform is image specific and its transformation matrix is a rotation matrix (see for example the KLT in Mathematica).

Given the values in your matrix, it seems to attempt to separate between intensity (first row) and color (the other two rows, spanning the remaining hyperplane) while being  optimized for speed.


Why do we always use square kernels for filters?

If the kernel is odd then you can center it on top of a pixel, which is nice. Lots of kernels are symmetric, so that's probably why they are often square. You can certainly use non-square or even kernels though.


Correct signal analysis flow (Wind Turbine vibration)?

I think the question as is, is unanswerable, but here I am making an attempt of correct it with these comments as answer.

  Should I upsample before or after rectifying?

For the algorithm, if you are finally assessing against a threshold in magnitude, you should not care about down|upsampling. Do not add that.

  Does this process look correct?

Strictly speaking, you are taking the algorithm from your document. Since your document has 200 pages, this is someway hardly to say if you are taking correctly the approach taken by the author.

  Should I do another LP filter to "clean" up the signal?

This depend on the specific fault you want to capture for the bearing mechanics. From there, you can decide which frequencies to focus, and how to assess them.

  Should I do another LP filter at a lower cutoff frequency to "clean" up the signal?

Again, please check specifically how your algorithm is built. Once you have it implemented, we could study modifications or alternatives.


From a programming point of view I don't see any possibility of avoiding the copy operation for 1024 values because the CUFFT library expects a pointer to an array containing NX (FFT size) values. If you don't initialize all of them, the result of the FFT operation will be random.

From a DSP point of view you could easily use a 32- or even 24-FFT without loosing information. Zeropadding in time domain corresponds to an interpolation in frequency domain, where the interpolated frequency bins can be different from the actual spectrum. You should check if you really need the frequency resolution of a 1024-FFT.


What is definition of independent random variable

No, NO and $\mathbf{NO}$.

Knowing that $E[X]=E[Y]=0$ and $E[XY]=0$ (or more generally that $E[XY]=E[X]E[Y]$ or equivalently, $E[XY]-E[X]E[Y]=0$) does not help in the least in proving or deducing that $X$ and $Y$ are independent random variables.  The conditions stated in the previous sentence are necessary for $X$ and $Y$ to be deemed independent random variables, but they are a very far cry from being sufficient for independence to hold. That is, if $E[XY]=E[X]E[Y]$, then we might suspect that $X$ and $Y$ might be independent random variables, but we cannot jump to the conclusion that they are independent random variables based on just the evidence that $E[XY]=E[X]E[Y]$.

  Definition: $X$ and $Y$ are said to be independent random variables (sometimes more prolixly mutually independent random variables) if their joint distribution function (joint CDF) $F_{X,Y}(u,v) = P\{X \leq u, Y\leq v\}$ satisfies
  $$F_{X,Y}(u,v)=F_{X}(u)F_{Y}(v) ~\text{for all real numbers}~ u ~\text{and}~v.\tag{1}$$
  Here, $F_{X}(u)=P\{X \leq u\}$ and $F_{Y}(v)=P\{Y\leq v\}$ are the marginal CDFs of $X$ and $Y$ respectively. It is important to remember the qualifier to the equality $F_{X,Y}(u,v)=F_{X}(u)F_{Y}(v)$: it must hold at all points in the $u$-$v$ plane.. Inanities such as "$X$ and $Y$ are independent when $X=5$ and $Y=3$ but not otherwise" are not permitted.

If $X$ and $Y$ are discrete random variables, then $(1)$ is equivalent to their joint _probability mass function (joint pmf) $p_{X,Y}(u_i, v_j)$ factoring into the product $p_X(u_i)p_Y(v_j)$ of the marginal pdfs
for all $i$ and $j$.

If $X$ and $Y$ are jointly continuous random variables, then $(1)$ is equivalent to their joint _probability density function (joint pdf) $f_{X,Y}(u, v)$ factoring into the product $f_X(u)f_Y(v)$ of the marginal pdfs
_for all real numbers $u$ and $v$.

Now, as the OP correctly points out, if $X$ and $Y$ are independent, then it is true that $E[XY]=E[X]E[Y]$, that is,
$$X, Y ~\text{independent random variables} \implies E[XY]=E[X]E[Y]$$ but

  the reverse implication
  $$E[XY]=E[X]E[Y]\implies X, Y ~\text{independent random variables}$$
  $$\mathbf{DOES~NOT~HOLD}$$

Still don't believe a word of all this malarkey? Consider discrete random variables $X$ and $Y$ such that
$$p_{X,Y}(1,0) = p_{X,Y}(-1,0) = p_{X,Y}(0,1) = p_{X,Y}(0,-1) = \frac 14.$$
Verify that

$X$ and $Y$ take on values $\pm 1$ with equal probability $\frac 14$ and value $0$ with probability $\frac 12$.
$XY=0$ always and so $E[XY]=0$ also, and that $E[X]=E[Y]=0$ also.
The joint pmf does not factor into the product of the marginal pdfs for at any of the four points where the joint pmf is nonzero.


How to compare the syllable extracted from the audio

A quick google for "bird voice recognition" yielded

https://peerj.com/articles/488/

which has an abstract describing that the authors refined the results they got using the old-school Mel-frequency cepstral coefficient comparison, as it has been common in (human) speech recognition for decades. I could start explaining what cepstrals and Mel-frequencies are, but there's so much literature out there that you literally just have to click read the first thing that hits you when you google for that. It's a relatively well-researched topic.


Circular convolution of a non causal signal

For an $N$-point circular convolution you can think of each signal as being periodically extended with period $N$. For your example with $N=4$ that would mean that the two sequences are

2 1 1 -1 and 2 -1 0 0

where both now start at index $n=0$. The result of the cyclic convolution is

5 0 1 -3

which is just a cyclic shift (by $2$) of the (correct) result that you obtained.


Difference between frequency sampling and windowing method

The Windowing method for filter design uses samples of the impulse response for the filter coefficients, truncated to the filter length (and then possibly, but not required, tapered by a higher performing window).

The Frequency Sampling method for filter design uses the IDFT (Inverse Discrete Fourier Transform) of the desired frequency response for the filter coefficients, using the number of samples across the frequency response to be equal to the number of coefficients in the filter.

The point of confusion is that in practice the impulse response is often approximated by an IDFT of the target frequency response, while an IDTFT (Inverse Discrete Time Fourier Transform) would actually be required. The IDTFT can be approximated by interpolating the frequency response first prior to taking the IDFT such that the time span of the impulse response increases. (Allowing the impulse response to sufficiently decay).

Therefore, if you take the IDFT of the sampled frequency response with the number of samples equal to the number of filter coefficients, you are using the Frequency Sampling method. If you take the IDFT of a higher number of samples in the frequency response, and then truncate to your filter length, you are approximating the Windowing method to the extent that the impulse response matches the actual impulse response for your target filter.  (For rectangular response filters there is no need to do an IDFT for the Window method; just sample a sinc function since that it easy to calculate and represents the actual impulse response).

(For further description of how the DFT and the DTFT are related see my response to:
Sine Wave Aliasing during IFFT)

This is a useful reference as well:

http://blogs.mathworks.com/steve/2010/03/15/the-dft-and-the-dtft/

Frequency Sampling

Very simple. Not very effective (except where selective precision is needed).

Coefficients of the FIR filter are the inverse DFT (discrete Fourier Transform) of the sampled frequency response desired.  This approach is very simple, but is not very effective for typical filter designs; much longer filter lengths are required over optimized approaches such as least squared and equi-ripple (Remez-exchange) filter design approaches. To be noted, the solution will be exact at the frequency sample locations, with ripple in between. When an application requires an exact solution at discrete frequency locations (and you do not care what happens in between) this can be a practical approach.

Windowing

Simple. Non-optimal.

Coefficients of the FIR filter are determined by an approach that approximates the inverse DTFT (Discrete-Time Fourier Transform); specifically the frequency response is interpolated (as a true IDTFT would be of a continuous frequency function), an IDFT is actually still calculated and the result of the IDFT (samples of the impulse response) is windowed by multiplying by a Window function of choice (I always liked the Kaiser window with its two "knobs" N and Beta). Multiplying the impulse response by a rectangular window is the estimation of the IDTFT process; as zero padding a time domain function will result in additional samples of the DTFT in the frequency domain.  The frequency domain is interpolated such that the impulse response determined from the DFT decays to insignificant values. Note that this results in the time domain samples that are samples of the actual ideal impulse response, simply truncated in time (this highlights the second view of what the difference is between the two approaches; here we see that prior to truncation and windowing, the impulse response used is samples of the continuous time impulse response since we over-sampled in the frequency domain to eliminate the effects of time domain aliasing.

Demonstration Example

Here is a specific example demonstrating the difference between Frequency Sampling and Windowing approach. Consider a target filter design with 99 taps and the following desired response for the normalized frequency range of $[0: N-1] * 2\pi/N$:

desired_resp= [ones(1,10) zeros(1,80) ones(1,9)];

For the frequency sampling method, the coefficients of the filter are simply the IFFT of the desired resp (rotated such that the dominant tap is centered in the filter, which would have naturally occurred in the IFFT if we had included the expected filter delay).

As shown in the frequency response for these coefficients, this will result in an exact solution at the frequency sampling locations (which can be useful in certain applications), but due to the fact that the impulse response has not fully decayed at the edges of the filter, there will be significant ripple in between. This also means that the impulse response shown is not truly samples of the impulse response for the target frequency response, as we are seeing the effects of folding in the time domain (equivalently to the effects of aliasing in the frequency domain when sampling in the time domain- it is the same process). This is specifically what differentiates the Windowing method of filter design!

In contrast, the Windowing method uses 99 samples of the actual impulse response. We can approximate this by interpolating our desired frequency response (at a higher sampling rate in frequency, which corresponds to a longer time span in time), or certainly in an exact approach, we can mathematically calculate the impulse response and sample that directly, skipping the IDFT operation all together [That in fact would be the exact definition of the Window approach, but in practice it can be sufficiently approximated by a IDFT if the target response is sufficiently interpolated. To demonstrate the difference, I have interpolated our target response by a factor of 10, resulting in the following estimate of the impulse response (989 samples) by taking the IDFT of the desired frequency response after it has been interpolated (note the result still suffers from time domain aliasing compared to using the formula for a sinc function as the impulse response which would be precise, but this is to show that even though both solutions use the IDFT, the result of taking an IDFT of 99 samples directly versus taking 99 samples of a 989 sample IDFT are very different):

Below is the comparison of the coefficients in both approaches showing how the difference in values occur at the tails; the Window method better approximates the actual impulse response, which is still decaying past the length of our filter. The Frequency Sampling method causes the values beyond the length of the filter to fold over, increasing the values of the coefficients, leading to increased ripple in the frequency domain.  Thus truncating the result (rectangular Window) using the Window method results in a higher performance filter, and further reducing the tails with a more advanced window function can lead to further performance at the expense of increasing the transition band.

The center 99 samples of this impulse response are chosen as the filter coefficients (this is the rectangular window). Below shows the result for both filter designs. Further, the performance of the windowing approach can be further improved by selecting the 99 taps used with a higher performance window such as Hann, Kaiser, Blackman-harris, etc.

Zooming in on the passband further shows how the frequency sampling method results in an exact answer at the frequency sampling locations but has overall more error than the window approach with the same number of taps.

Note on FIR Filter Design related to the Windowing and Frequency Sampling Approach

I believe the above two approaches are insightful, but I would rarely use either for an actual filter design. It would add interest if readers commented on applications where either of these two approaches are actually the "best choice". My design approaches of choice are the Least Squared and the Remez-Exchange (equi-ripple) algorithms. Most of the time I go with Least Squared as that results in the best statistical performance (best SNR solution), while the Remez will meet a target specification consistently throughout the passband and stopband.

See below plots of comparative results for a 31 tap FIR filter solution with a 9.6KHz passband and 12KHz stopband showing how a Window design (using the Hamming Window) compares to Least Squared, and how Equi-ripple (Remez) compares to Least Squared; giving insight into why I typically go with Least Squared.

Also see Matt L's response to Difference between Equiripple & Least Squares design for FIR digital filters.


Understanding the relation between instantaneous frequency and phase of a sinusoid

An analogy: a car is at $x=0$ at $t=0$ and moves with velocity $v$ (constant). The famous formula for this is:

\begin{equation}
    x = v \cdot t
\end{equation}

By integrating velocity $v$ one gets position $x$, by derivating position $x$ one gets velocity $v$.

Now, take position to be a point in a circle (for this is convenient to use radians instead of raw meters along the circumference, but this does not really matter). Take velocity to be how fast that point moves around the circle (for this is better to choose an unit like radians per second to make it independent of the circle radius, but this does not matter neither). This would be the phasor with constant rotational speed, a.k.a sine wave. So, we can reinterpret the formula like:

\begin{equation}
    \theta = \omega \cdot t
\end{equation}

And the same logic as above applies: By integrating frequency $\omega$ one gets phase $\theta$, by derivating phase $\theta$ one gets frequency $\omega$.

You can verify by dimensional analysis that they both are the same formula actually (well, actually, radian is said to be a dimensionless unit, something like a 'dozen' or a 'mole'... that is another story):

\begin{equation}
[Length] = [Length/Time] \cdot [Time]
\end{equation}

When you integrate with respect to time, one less $[Time]$ appear in the denominator. The derivative does the opposite.


In opecv function Hough circles how does parameter 1 and 2 affect circle detection and how can I adjust them to increase accuracy?

Intuition for parameters of HoughCircles:

image: 8-bit, single channel image. If working with a color image, convert to grayscale first.

method: Defines the method to detect circles in images.  Currently, the only implemented method is cv2.HOUGH_GRADIENT, which corresponds to the Yuen et al. paper.

dp: Resolution of the accumulator array.  Votes cast are binned into squares set by dp size.  Set too small and only perfect circles are found, set too high and noise collaborates to vote for non-circles.

minDist: Minimum distance between the center (x, y) coordinates of detected circles. If the minDist is too small, multiple circles in the same neighborhood as the original may be (falsely) detected. If the minDist is too large, then some circles may not be detected at all.

param1: Is a number forwarded to the Canny edge detector (applied to a grayscale image) that represents the threshold1 passed to Canny(...). Canny uses that number to guide edge detection: http://docs.opencv.org/2.4/modules/imgproc/doc/feature_detection.html?highlight=canny

param2: Accumulator threshold value for the cv2.HOUGH_GRADIENT method.  The smaller the threshold is, the more circles will be detected (including false circles). The larger the threshold is, the more circles will potentially be returned.

minRadius: Minimum size of the radius in pixels. Don't set minRadius and maxRadius far apart unless you want all possible circles that might be found in that range.

maxRadius: Maximum size of the radius (in pixels).  Don't set minRadius and maxRadius far apart unless you want all possible circles found in that range.

If you don't know radius, you must auto-tune hyper parameters:

If you don't feed the right radius (+- 15 pixels) to HoughCircles, the output is next to useless because it returns  circles everywhere and nowhere without rhyme or reason.

However there is a trick if you don't know the radius and don't know the optimal dp resolution and accumulator array threshold.  HoughCircles is quick (linear time?), put HoughCircles in a loop, start out with a big minDist, fine dp resolution, and high accumulator voting threshold, and start guessing the largest radius of circle you're willing to find, HoughCircles returns nothing.  Each iteration, slowly decrease radius, decrease threshold, increase dp resolution.  Eventually HoughCircles will find a pixel-perfect circle.

Do this process a second time and approach from the other side, loop HoughCircles starting out guessing a minimum radius, low resolution, high voting threshold, eventually the smallest and highest quality circle will be found.  If the two circles agree on a common circle, you can be very confidant the circle found is high quality, otherwise increment/decrement the hyperparameters at different ranges and speeds.

HoughCircles is mathematically precise and exactly correct, it has been proven to find pixel perfect circles even when you don't know radius.  With hyper-parameters auto-tuned correctly for your images, the circles found are pixel perfect.

Not only can you use HoughCircles to find the middle of the circle, you can use it to define the outer lip and inner lip of the circle's edge gradient (sub pixel accuracy). For example, the outer lip or inner lip of a coin's circumference.

Many passes of many iterations must be used to auto-tune hyper parameters uniquely for every image to slowly discover which circles are the one you want.

Prove pixel perfect accuracy:

If performance must be proved as sub-pixel perfect for an industrial setting, you can't get around intuitively understanding what the algorithm is doing so you can fiddle with its implementation to suit your needs.  For in depth guide, See: "Computer Vision: A Modern Approach" by David Forsyth and Jean Ponce, 2nd Edition, Chapter 10 on the Hough Transform.


Confusion about shannon's entropy of the grayscale image after lossy compression

Shannon's entropy only makes sense when you define it over a probability distribution. You're using a probability distribution where each pixel is sampled independently from the histogram of pixel intensities. Clearly, real pictures don't come from such a distribution. If they did, they would all look like white noise.

Also note that individual elements don't have entropy. Only distributions have entropy. So if you have two images it doesn't make sense to compare their entropy. It only makes sense to compare the entropy of the random distributions that you sampled the images from.


Should I ever pick the continuous wavelet transform over the discrete one? DWT vs CWT vs STFT

Standard 2-band discrete wavelet transforms have some subsampling, thus they tend to be much less precise in localization. Plus, their shapes are limited, and the ones with finite support are slightly dissymmetric (all but Haar), which is problematic for symmetric pulse/edge detection.

If you don't care about computations, and are interested in visualization, matching, feature detection, CWT might be easier.
If you are interested in  denoising, compression, restoration, DWT are often more appropriate.

Yet, there are many useful schemes bridging the gap between DWT and CWT that can be used in dedicated applications.

A reference paper: Paul S. Addison, 2005, Physiol. Meas. Wavelet transforms and the ECG: a review
References for fast CWT implementations
Matlab code: continuous wavelet transform and inverse
C code: CWT Library: Continuous wavelet transform library
C++ code: Wavelet - A library for online estimation of the Continuous Wavelet Transform


What tools do you use to create signal processing plots?

Extending a bit on the answer of @Rob. If you are writing a paper in LaTeX (and friends), it is (imho) highly desirable to draw the diagrams also in a LaTeX related way. This comes with a number of advantages:

Direct and easy integration, using either \input or including as PDF using \includegraphics
The same fonts, font sizes, colors, etc. as in the rest of the document (main advantage in my opinion)
One workflow. You can put everything under the same version control (e.g. Git, Mercurial, ...) and do not require any other software: only a LaTeX installation - which is needed anyways.
No need to manually move stuff around using the arrow keys, no automatically changing grids, no vector-based programs which can't export eps, ... (yes, I'm looking at you, Visio!)

Now there are a number of LaTeX related tools for this. Rob mentioned TikZ in his answer. There are also Metapost and PSTricks which also have the capabilities to get the job done. If you want to stick closely to the LaTeX syntax you seem to already know, you might want to chose TikZ.

To do especially plots using TikZ, you might want to have a look at the PGFPlots package, which offers extensive plotting capabilities. All options are well described in the manual, and there is a huge example gallery. For example to get a plot as the first one in your diagram, all you need is

\documentclass{article}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}
\begin{axis}[
        width = 8cm, height = 4cm,
        xlabel={$\omega$},ylabel={$X(\omega)$},
        xtick = {-1,1},
        xticklabels={-B,B},
        ytick = {0,1},
        axis lines=middle,
        xmin = -4, xmax = 4,
        ymin = 0, ymax = 1.2
    ]
        \addplot[] coordinates {(-1,0)(0,1)(1,0)};
\end{axis}
\end{tikzpicture}
\end{document}


Issues in generating AR model with a constraint

I was able to produce the AR(p) model with much of zero-coefficients. I use Scilab rather than Matlab, therefore I'll just describe the technique:

Generate some appropriate data with the desired spectral properties - the AR model will mimic these spectral properties.
Fit the AR model of the desired order to your data.
Find the desired number of the coefficients with the lowest values and reset them to zero  (something like a(find(abs(a) < stdev(a) / 2)) = 0 - valid in Scilab).

The resulting model should (may?) still has a similar spectral properties & satisfy your requirements considering zero-coefficients.


You have to clearly define what you mean by an "analog FIR filter". "No poles" is not correct because (discrete-time) FIR filters do have poles; they are just all at the origin of the $z$-plane (for causal FIR filters). Note that filters without poles do not exist. Take as an example the discrete-time transfer function

$$H(z)=1-az^{-1}\tag{1}$$

with a zero at $z=a$. Obviously, $H(z)$ has a pole at $z=0$. Similarly for

$$H(z)=z-a\tag{2}$$

which also has a zero at $z=a$, but a pole at $z=\infty$.

If you mean transversal filters implemented as a tapped delay line, then you can find analog versions of it, see e.g. this paper.

If you literally mean a filter with a finite impulse response (FIR) then, technically speaking, such a thing does not exist in the continuous-time world. But of course, any stable system has an impulse response that decays sufficiently fast such that for practical purposes it may be considered as zero (or unmeasurable) after a certain time. The discussion if such a system should be considered as FIR or not is moot.


In your case you probably want to calculate the SNR as mean over standard deviation.

signal=rand([256,192,330]); %demo data
SNR = mean(signal,3)./std(signal,[],3);
SNRdb = 10*log10(SNR);

this way you obtain different SNR values per pixel. 256x192 pixels in 330 frames.

To get the values for each of the 330 frames instead you must first reshape your matrix.

signal=reshape(signal,[],330);  %pixel values merged
SNR = mean(signal)./std(signal);
SNRdb = 10*log10(SNR);


Using cross-correlation and finding the peak is one way of finding the time delay between two signals.
This assumes that the two signals are periodic and have the same fundamental frequency.

When the two signals have close, but not equal fundamental frequency, the lag position of the maximum cross correlation value will change in time.  If the lag of the maximum does not change (and some thresholds are satisfied so you're not just cross-correlating two very weak signals with the same 60 Hz hum in the background), the two signals are synchronous.  If it changes slowly, that indicates the degree that the two signals have different frequencies.  If the lag requires T seconds for the lag to slide one period (which would be the distance between lags of maximum correlation) in one second, the two signals differ in frequency by 1 Hz.


Fast Fourier Transform for varying period (Order analysis)

There are a few things going on here:

  Does this thought process look correct?

Strictly, No.

The main reason for this is that you are not performing a resampling step that is the differentiating point between the Discrete Fourier Transform (DFT) and order analysis.

  Is it correct to say that the final fft plot is in the frequency domain? From this link, ...

Yes, it is correct to say that the final FFT plot is in the frequency domain. But, the order plot would be in a quantity that has no units. It is basically a multiplication factor.

A little bit about Order Analysis and how it is different to the DFT:

OK, so you have a rotating shaft which you sample at some $Fs$ using an accelerometer. Let's use the figure provided, at 25kHz.

If the shaft is rotating at 1 Revolution Per Minute (RPM) then 1 revolution fits in $60 \cdot 25000$ samples. If you continued sampling for a full 2 minutes, you would have $2 \cdot 60 \cdot 25000$ samples that would also happen to describe 2 revolutions. The data for the first revolution are in indices 0-1499999, the data for the second revolution are in indices 1500000-2999999 and so on.

Now, we throttle the motor up and get the shaft at 2 RPM. Now, $60 \cdot 25000$ samples describe 2 revolutions! The data for the first revolution are in indices 0-749000, the data for the second revolution are in indices 750000-1499999 and so on.

More importantly, what has happened from the point of view of the data that we are interested in (the data about 1 revolution), the sampling frequency has been halved.

At 1 RPM we observed $60 \cdot 25k$ samples for one revolution. At 2 RPM we observed $\frac{60 \cdot 25k}{2}$ samples for one revolution. At 3 RPM...and so on.

At 1200 RPM, one shaft rotation takes $\frac{1}{1200} \cdot 60 = 0.05$ seconds and fits in 1250 samples (at 25kHz).

So what?

Imagine that there is a "bump" (for whatever reason, a bearing, uncentered load, structural resonance, whatever). At 1200 RPM, you can hear that thing banging across the floor. But at 1 RPM, the discontinuity is so "shallow" that it is spread across millions of samples, if it is even there. At low RPM, the elasticity of the materials might allow the power of this discontinuity to scale linearly. But at high RPM, the non-linearities kick in. Now, the shaft cannot bang off of its housing (it can, but by then it is too late) and you can see this because the structure will give rise to additional harmonics. So, it is not only that at high RPM things become louder, but also, new things show up, that were not observable at low RPM.

So, we have a problem here, because with increasing RPM we are trying to fit more information at increasingly less and less and less bandwidth.

Part of the answer here is to increase the sampling frequency.

But, we would still be left with no way to normalise the data for the increase in speed. Therefore, if we were to run DFT on a block of data, we would not be sure how much of a rotation or where does interesting behaviour occurs.

One approach to this would be to physically modulate the sampling frequency by the RPM.

This is what Order Analysis does, but it does it in software. What you are doing prior to DFT is resampling the signal (as a function of RPM), so that, no matter what the RPM does (going up or down), 1 shaft revolution sits in exactly the same amount of samples.

By doing this, you can observe consistent parts of the spectrum as the RPM increases (e.g. harmonics that are predictable as per the linked article) and whether interesting components start to appear at specific RPM.

A fundamental assumption of the DFT is that the signal at its input is periodic. So, although you may be passing it a voice recording of someone saying "one, two, three", what the DFT "sees" is "one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, one, two, three, ..." and analyses this as a whole.

This is the reason why the spectrum looks like two peaks connected with a curving ribbon on the first spectrum from the figure of the book you are posting.

What has happened in the second picture is to adjust the timing of the samples based on the RPM, so that each revolution now takes up exactly the same number of samples.

There is nothing additional hapenning to that sinusoid, except modulating its frequency in a way proportional to the RPM "signal". If you adjust for the RPM, what you see is basically just the same sinusoid.

A better way to visualise this is via the spectrogram. Here is a sweeping sinusoid:

And here is its spectrogram:

(Both images are from the Wikipedia article).

What the spectrogram is, is repeated applications of the FFT in short overlapping blocks throughout the signal. Therefore, its x-axis is time, its y-axis is frequency and the cell colour is proportional to the strength of that harmonic.

If you apply a spectrogram to your data, you will notice that the "run-up" (or run-down) shows up clearly as increasing (or decreasing) "lines". The reason why there might be more than one lines is because there are additional components to the rotation that also happen faster (or slower), depending on the RPM.

What we really want is to get rid of the "run-up" (or run-down) component of the DFT.

What we are interested in is the spectrum associated with 1 shaft revolution at some fixed RPM. This is equivalent to one vertical slice of the spectrogram at some time instant (which however in our case is proportional to RPM).

With order analysis, what you are doing is basically "rotating" the spectrogram, at an angle that matches this increasing line. If you did this, in this particular example, you would see a straight line. (Technically, it is not a rotation, but a "skewing", but just for illustration purposes, think of it as rotating this line so that it is parallel to the x-axis).

This is what this re-sampling (mentioned in the figure you attach) refers to.

After you do that, the x-axis of the spectrogram would be RPM, the y-axis of the spectrogram would be Order and the cell colour would be proportional to the strength of that harmonic.

For more information please see here and here

Hope this helps.


Image 2D Real Cepstrum with DFT, Is `ifftshift` Needed?

If you have image in the spatial domain, in order to calculate its DFT transform you should use fft.
Once you have its DFT in order to get back to the spatial domain use the function ifft.

Either way, the DFT transform gives you the data in the [0, 2pi] axis.
Use fftshift to move it into the [-pi, pi] domain as seen in the paper.

Good Luck.

MATLAB Code

This code worked for me:

mInputImage = imread('3.png');
mInputImage = double(mInputImage) / 255;
mInputImage = mean(mInputImage, 3);

mImageDft = (fft2(mInputImage));

mImageDft = log(1 + abs(mImageDft));

mImageCepstrum = abs(fftshift(ifft2(mImageDft)));

mImageCepstrum = mImageCepstrum / max(mImageCepstrum(:));

figure();
imshow(mImageCepstrum);

Enjoy...


Show That the Power Spectrum Density Matrix Is Positive Semi Definite (PSD) Matrix

Pay attention that for Scalar Random Process the Power Spectrum Density is non negative.

Namely, let $ y \left[ n \right] \in \mathbb{R} $ a WSS Random process with its Auto Correlation function given by:

$$ {R}_{y, y} \left[ m \right] = \mathbb{E} \left[ y \left[ n \right] y \left[ n - m \right] \right] $$

Then the Power Spectrum Denisty is:

$$ {S}_{y, y} \left( f \right) = \sum_{m = -\infty}^{\infty} {R}_{y, y} \left[ m \right] {e}^{-j 2 \pi f k} \geq 0 $$

Then, by defining $ z \left[ n \right] = \boldsymbol{v}^{T} \boldsymbol{x} \left[ n \right] $ one would get:

$$ \begin{align*}
0 \leq {S}_{z, z} \left( f \right) & =  \sum_{m = -\infty}^{\infty} \mathbb{E} \left[ z \left[ n \right] {z}^{T} \left[ n - m \right] \right] {e}^{-j 2 \pi f m} \\
& = \sum_{m = -\infty}^{\infty} \mathbb{E} \left[ \boldsymbol{v}^{T} \boldsymbol{x} \left[ n \right] \boldsymbol{x}^{T} \left[ m - m \right] \boldsymbol{v} \right] {e}^{-j 2 \pi f m} \\
& = \boldsymbol{v}^{T} \left( \sum_{m = -\infty}^{\infty} \mathbb{E} \left[ \boldsymbol{x} \left[ n \right] \boldsymbol{x}^{T} \left[ n - m \right] \right] {e}^{-j 2 \pi f m} \right) \boldsymbol{v} \\
& = \boldsymbol{v}^{T} {S}_{x, x} \left( f \right) \boldsymbol{v} \\
& \Rightarrow {S}_{x, x} \left( f \right) \succeq 0
\end{align*} $$

Some remarks regarding simulating it in MATLAB:

If one use MATLAB's xcorr() to calculate the Auto Correlation one should use iffstshift() to shift the function to be "Symmetric" to MATLAB in order to have the DFT Real and Non Negative. This is due to the fast MATLAB's fft() expects the first sample to be of index $ 0 $ (See the lags output of xcorr()).


Continuous Time Signal and Discrete Time Signal - Connection Between Periodicity and Discretness

This is a 2 parts answer:

Given a discrete representation of the Spectral Data does necessarily represent a periodic function in the time domain?
Well, it is by definition. Look, the DFT is a fancy name to a scaled Discrete Fourier Series, namely combination of harmonic signals. By the nature they are all an integer multiplication of a base frequency. Now think about a linear combination of finite harmonic signals which all are multiplication of known frequency -> The are all periodic in the lowest frequency. And so is the signal.
Can all periodic signals be represented by a Fourier Series?
Now you are asking about convergence and there are few requirements about signals to be represented by a Fourier Series which deals with the ability to integrate them over the period. Another issue, when you say Fourier series, remember an infinite number of signals are considered.


Compression of Multivariate Phase Signals

Your signal is easily described by two parameters (the initial value and the frequency), the sawtooth may be linearized, by stitching each linear segment to the end of the other, and then you can use a linear regression.
The function that stitch the linear segments of a sawtooth of amplitude $2 \pi$ is the unwrap, there is a function with the same name in octave/matlab if you want to use them.
Here is a code snipped that you could use to get the frequency of your columns in python
import numpy as np
omega = np.arange(2, 8).reshape(1, -1) # a row vector
t = np.linspace(0, 10, 1001).reshape(-1, 1) # a column vector
# this is a matrix as you described
matrix = np.angle(np.exp(1j * omega * t))
linearized = np.unwrap(matrix, axis=0);
# this function will get the frequency of your underlying signal
omega_ = np.mean(np.diff(linearized, axis=0), axis=0)

Then for the phase you could use
phase = np.mean(linearized - omega_ * t)


Problem solved.

Looks like calculation of Os ($\omega_s$) was the culprit for this problem. Proper equation is:

Os=Om*sqrt((gm^2-g1^2)*(1-gm^2))/(1-gm^2);

instead of:

Os=Om*(sqrt((gm^2-g1^2)*(1-gm^2))/(1-gm^2));