number of subcarriers for an OFDM system

You've got that right – when designing an OFDM system, you pick the number of subcarriers usually just large enough to make the subchannels look flat.

You avoid overly large subcarrier numbers, because more subcarriers means longer symbols, and you need to restrict the symbol length to stay below coherence time of the channels.

Now, 3 is indeed a choice that is very small. I wouldn't build an OFDM system with 3 subcarrieres, simply because the smallest cyclic prefix that you can build is then $\frac 13$ of the symbol length – and chances are, you need more than one $\frac1{B_{subcarrier}}$ as guard interval to "capture" the whole channel impulse response.

Then again, your channel seems to be very benevolent, since it has such a high coherence bandwidth, so maybe that's really the case.

However, in a situation where you'd do 3- or 4-subcarrier OFDM, you'd simply not do OFDM alltogether and would probably either simply use a single-carrier system with a proper equalizer (might be easier, even), which will not waste as much spectrum as OFDM, or you'd go for filterbank multicarrier systems, because you will have a very bad time¹ trying to achieve synchronicity in such an OFDM system. Schmidl&Cox will be very ambiguous prior to first lock, and the SNR of any timing symbols will be less-than-good with just three values to work with per OFDM symbol.

¹ yes, that's a pun.


Transfer function of a nonhomogeneous difference equation

The last difference equation is not a linear system due to the addition of the constant $\gamma$, therefore it does not have a transfer function.


How to remove the boundary effects arising due to zero padding in scipy/numpy fft?

Can there be any changes made in the above method to solve this
  boundary problem?

Yes there can be, but it all depends what you want to occur at the boundary. If you don't want it to be smoothed then simply replace the edge with the original signal for a sensible number of samples, depending on the width of your gaussian. You gaussian might be > 0 for say 10 samples, so replace the 5 samples at the start and end of the signal with the original samples.

  Can you suggest a way of removing this zero padding in the scipy fft
  source code?

Yes. Pad it yourself. If you want to consider the wrap around, as if your signal is truly periodic, then fabricate the wraparound effect for a power of 2 in size by carefully padding your own signal. Just as a simple idea, multiply your window by 2 and find the nearest power of 2, then fill the middle of the copy with zeros to pad it out

[1234567890] [12345] - zero pad here - [67890]

Now when you apply the convolution, you should get the same effect for the first window's worth of the result as if the convolution kernel had wrapped around the original signal. To illustrate, take your gaussian and zero pad it at the end, and slide through the convolution to see for yourself that it yields the same thing

 original   left    padding    right
[1234567890 12345 000000000000 67890]
[xgaussianx 00000 000000000000 00000]

so the value for the right most sample comes from this:

 original   left    padding    right
[1234567890 12345 000000000000 67890]
[sianx00000 00000 000000000000 xgaus]

Which is equivalent to what we actually want

[1234567890]
[xsianxgaus] (gaussian shifted by 5 samples and wrapped around)

and the value for the left most sample comes from

 original   left    padding    right
[1234567890 12345 000000000000 67890]
[00000xgaus sianx 000000000000 00000]

Which is also equivalent to what we actually want

[1234567890]
[sianxxgaus] (gaussian shifted by -5 samples and wrapped around)

All this relies on deciding what it is you actually want to happen. If you understand the mathematical process of convolution, and you understand exactly what it is you want to calculate, you can fabricate the conditions under power of 2 scenarios to ensure you get fast computation. Generally speaking most people just opt for the replace because their signal windows are not periodic and have sharp discontinuities at the edges.


Measure 3-D length using smart-phone's accelerometer

Your model may be misleading, as bias and noise (including non-linearity and quantization noise) get double integrated, thus error increases quadratically and without bound over time.


The higher the sidelobes of a signal frequency spectrum, the greater the "cross-talk" between transmission channels

Convolution theorem states convolution in time domain is the same as multiplication in Fourier domain, also the multiplication in time domain is a convolution in Fourier domain.

So when you truncate a signal (multiply it by a window), your signal's spectrum will be convoluted with the window's spectrum which widen your signal's spectrum depending on the main lobe width and also replicate signal spectrum depending on the side lobe height which is called frequency leak.
if you consider different channels over different bandwidth, when you truncate the signal the frequency leak will cause channel cross-talk.


I don't think you can completely get rid of histogram or threshold-based binarization, since the former is to achieve line segmentation, while the latter is to extract the letters. The Radon horizontal projection is used for line segmentation, and the center line can be used to approximate the baseline of each segment. Yet this is somehow equivalent to the horizontal histogram.

There are some improvement based on the basic histogram. In Zheng et al's paper  A model based line detection in documents, they applied a directional single connected chain strategy that considers the global textline orientation and merges the non-overlapped small segments into one line. The algorithm provides more accurate base line detection than the pure histogram-based line segmentation.

Besides, in Arabic script where short characters and large diacritics exist, the Tarik et al's paper A Novel Baseline Detection Method of Handwritten Arabic-Script Documents Based on Sub-Words that measures the vertical distance between the semi-overlapped subwords baselines.

Yet in those studies in the attempt to increase the accuracy, histogram and/or threshold were still a necessity.


What is meant by "arbitrary" in the context of digital filter design?

The term "arbitrary" in this context refers to the desired frequency response, and what it means is that the respective design method accepts any desired response. So unlike required by many standard design routines, the magnitude response does not need to be piecewise constant or linear, and, more importantly,  the phase response need not be linear or minimum-phase. The corresponding design problem is a complex approximation problem because the desired frequency response which must be approximated is a complex-valued function.


For your second question, valid means the range of output where there is full overlap between the convolving signals.

Note that for two causal signals, $x[n]$ of length $L_x$ and $h[n]$ of length $L_h$ the convolution output $y[n]$ will be of length $L_y = L_x + L_h -1$.

However among those output $y[n]$ samples, the first $L_d = \min\{L_x,L_h\} -1$ and the last $L_d= \min\{L_x,L_h\} -1$ samples will be computed with partial overlaps between $x[n]$ and $y[n]$ (whicih are the edge samples) but the remaining portion (that's in the center) will be be computed by full overlap with the shorter sequence being fully inside the longer sequence.

In your case; $L_x = 4$ and $L_h=3$ then the first and last $2$ samples will be discarded and the center $2$ samples will be returned.


If we focus on finite-duration impulse response filters and apply the
MSE criterion to optimize the coefficients of the feedback and feedforward filters.

The equalizer output can be expressed as:
\begin{equation}\hat{I}_{k}=\sum_{j=-K_{1}}^{0} c_{j} v_{k-j}+\sum_{j=1}^{K_{2}} c_{j} \tilde{I}_{k-j}\end{equation}

where $\hat{I}_{k}$ is an estimate of the $k$ th information symbol, $\left\{c_{j}\right\}$ are the tap coefficients of the filter, and $\left\{\tilde{I}_{k-1}, \ldots, \tilde{I}_{k-K_{2}}\right\}$ are previously detected symbols. The equalizer is assumed to have ( $K_{1}+1$ ) taps in its feedforward section and $K_{2}$ in its feedback section.

The MSE criterion result in a mathematically tractable optimization of the equalizer coefficients; based on the assumption that previously detected symbols in the feedback filter are correct, the minimization of MSE:
\begin{equation}J\left(K_{1}, K_{2}\right)=E\left|I_{k}-\hat{I}_{k}\right|^{2}\end{equation}
This leads to the following set of linear equations for the coefficients of the feedforward filter:

\begin{equation}\sum_{j=-K_{1}}^{0} \psi_{l j} c_{j}=f_{-l}^{*}, \quad l=-K_{1}, \ldots,-1,0\end{equation}
where
\begin{equation}\psi_{l j}=\sum_{m=0}^{-l} f_{m}^{*} f_{m+l-j}+N_{0} \delta_{l j}, \quad l, j=-K_{1}, \ldots,-1,0\end{equation}

The coefficients of the feedback filter of the equalizer are given in terms of the coefficients of the feedforward section by the following expression:

\begin{equation}c_{k}=-\sum_{j=-K_{1}}^{0} c_{j} f_{k-j}, \quad k=1,2, \ldots, K_{2}\end{equation}

The values of the feedback coefficients result in complete elimination of intersymbol
interference from previously detected symbols, provided that previous decisions are
correct and that $K_2 > L$


What does the one-tap equalizer used with OFDM accomplish?

The OFDM signal as a whole is affected by frequency selective filtering. It is usually designed such that the subcarrier bandwidth is smaller than the channel coherence bandwidth. This yields you frequency flat fading for each subcarrier which can be described by a single complex multiplication and equalization can equally done with just one tap.

Equalization in OFDM usually happens in the frequency domain. Let $X$ be the vector of subcarriers in the frequency domain and $x=IDFT(X)$ the corresponding time domain signal. Then, if the impulse response of your channel is a vector $h$, the received signal $y$ will be $y=x*h$ where $*$ referes to convolution. Then $Y=DFT(Y)=X \cdot H$ with $H=DFT(h)$ and $\cdot$ refering to element-wise multiplication. However, this is only true if $y$ is the circular convolution of $x$ and $h$ while in reality the convolution is linear. Therefore, the cyclic prefix is prepended (it needs to be at least as long as $h$) so that the linear convolution contains the circular convolution.

Then, $Y=DFT(Y)=X \cdot H$ with $H=DFT(h)$ and $\cdot$ as stated before. Since each element of $Y$ and $X$ refers to one subcarrier it can be seen that each subcarrier is multiplied with just one complex element from $H$ and hence equalization just becomes $Z=Y/H$ where again $/$ referes to element-wise division. This describes zero-forcing equalization which might not be the best choice if additive noise and interferers are present but it transports the idea.

And this is nothing else but one-tap equalization for each individual subcarrier.


Background music and human voice occupy the same frequency ranges, so they cannot be separated by filters. See this Q&A for a few suggestions of techniques to separate human voice from background music.


It's easiest to look at this in the frequency domain. We have
$$\mathcal{F} (x*y) = X \cdot Y$$
$$\mathcal{F} (x \star y) = X' \cdot Y$$
where ' is the complex conjugate
So we get
$$\mathcal{F} (f \star (g * h)) = F' \cdot G \cdot H $$
$$\mathcal{F} ((f \star g) * h) = (F \cdot G)' \cdot H = F' \cdot G' \cdot H  $$
These are only identical if G is real.
In your case you let's simplify the notation to $f = S_0, g = S_0, h = H$ $S_0$ has a unit impulse as autocorrelation, which means that the spectrum of the autocorrelation is just unity, i.e. $1$
$$\mathcal{F} (f \star f) = \mathcal{F} (\delta) = 1 = F' \cdot F$$
We use that property and can show that
$$\mathcal{F} (f \star (f * h)) = F' \cdot F \cdot H = H$$
since $F' \cdot F = 1$


Intuition behind the Gaussian Filter in Image Processing

As far as your comment suggests, you are using a LP Gaussian Filter.
LP = Low pass.
This means it passes the low frequency component of the image, and prevents the high frequency component.

In image processing, you can imagine the high frequency component as the places in the image where you have a sudden jump in pixel values. Like moving from a black background and suddenly to white one.

These components are associated with edge detection normally, but at the same type they are associated with noise. It's because noise has a high frequency.
How come?
There is for example a type of noise known as "Salt and Pepper" where there appears spikes in the image. You see the it as a weird dot in the image.
These normally have way different value than the surrounding.

With this sort of a filter (Low Pass Gaussian Filter) you can suppress the noise, but you will lose some sharpness in the image due to killing some parts of the edges.
Therefore; this is called smoothing. It's because you tend to make the image smooth in a way that pixel values are not changing rapidly but rather in an elegant way.

You can imagine this as averaging. It's because you are substituting the values in a certain region of the image by it's average. That will help in smoothing the image.


As pointed out by niaren, you have to multiply by two, but this also means that your whole magnitude response is shifted up by 6dB, also in the stopbands. So you will have 6dB less stopband attenuation. Of course, the attenuation relative to the passband gain remains constant. You have to take this into account when designing your filter. Say you need 60dB stopband attenuation in absolute terms, not relative to the passband gain, then you need to design a 0dB gain filter with 66dB stopband gain. After multiplying the coefficients by 2 you'll get 6dB passband gain and 60dB stopband attenuation.


Why is the size of results from FFT half the size of the input, while that is not the case in image processing?

No.

Size of FFT in 1D is not half of the signal length. It's still the signal length (assuming FFT length was of the signal length). However for real data, due to conjugate symmetry, half of this is redundant and can be dicarded to reduce storage. Note however that, when it's required to process the FFT data, you have to construct the full length from the half and proceed.

This is also the reason why spectral analysis software display only positive frequencies for real data; i.e, negative frequencies will be a mirror copy (of the magnitude) of positive frequencies.

But you don't have to discard the half. You can just retain it.

For image processing, the symmetry of FFT for real input data again exist and if you wish you can also dicard half of the image FFT data. Whether this will be employed or not depends on the intentions of the package.


Basically you have to do a shading correction. This process removes known artifacts from an image (or a series of images) based on a reference image (or "compensation image" as it is sometimes called) containing only those artifacts. Such a reference images can be obtained by photographing a uniform object (i.e. a perfectly white wall or perfectly black darkness) under good conditions so that only the signal from the defect camera chip is recorded.

I have done this example with Mathematica, the basic principle can be implemented anywhere.

Given a set of images (I used grayscale images since I have no infrared images) without the defective signal

img1=

img2=

img3=

and an artificial defect camera chip signal aka the reference image (I pulled this out of my hat).

noise = GaussianFilter[Binarize[Image[RandomReal[1, {299, 499}]], 0.95], 5]

Now we generate the noisy images.

img1n = ImageAdd[img1, noise]
img2n = ImageAdd[img2, noise]
img3n = ImageAdd[img3, noise]

Then we simply subtract the reference image from the noisy images and get the actual images without the noise.

ImageSubtract[img1n, noise]
ImageSubtract[img2n, noise]
ImageSubtract[img3n, noise]

If you happen to have a noisy image and a perfectly fine image with exactly the same scene (and of course image sizes)

img1n =

img1=

then you can simply generate the reference image yourself by subtracting the latter from the former.

ImageSubtract[img1n, img1]


Teager-Kaiser Operator vs. Hilbert Transform

This paper may be of interest:

David Vakman, "On the Analytic Signal, the Teager-Kaiser Energy Algorithm, and Other Methods for Defining Amplitude and Frequency." IEEE Trans. Signal Processing. (1996)

Summarising from the paper:

$$\Psi(u) = a^2w^2 = [u'(t)]^2 - u(t)u''(t)$$
$$\Psi'(u') = a^2w^4 = [u''(t)]^2 - u'(t)u'''(t)$$

where $u(t)$ is the signal, and $a$ and $w$ iare the amplitude and frequency estimates to solve for.

$$a(t) = \frac{\Psi(u)}{\sqrt{\Psi(u')}}$$
$$w(t) = \sqrt{\frac{\Psi(u')}{\Psi(u)}}$$

But this means if at some point in the signal $u'(t)=1$ and $u''(t)=u''''(t)=0$, for example, then $\Psi(u) = 1$ and $\Psi'(u) = 0$ and thus $a(t) = \infty$ and $w(t) = 0$.

The paper thus explains this causes spikes in amplitude and zero frequencies with the Teager-Kaiser algorithm.

For your signal, the reverse must be true. If $\Psi(u) = 0$ and $\Psi'(u) = 1$ then you would get spikes in frequency when amplitude is zero. This would occur when $[u'(t)]^2 = u(t)u''(t)$

The Hilbert transform uses the entire signal to compute. It is implemented by doubling the positive frequencies and setting the negative frequencies to zero. Therefore the energy of the signal does not change (Parseval's theorem). Thus, intuitively, there can be no infinite spikes, as this would require infinite energy.


Remove background noise from audio file (python or matlab)

Since you're adding white noise, the highpass and lowpass filtering will almost  not remove the noise in the frequency band where you want to keep your signal, so you will always have some background noise with this highpass and lowpass filtering strategy.

Not sure if this helps, it depends on the signal-to-noise ratio:
If you can clearly distinguish the noise from the signal in the spectrum (something similar as in the second figure of the Noisy Signal example in Matlab's documentation of the fft), you could set a threshold and make the spectrum with an amplitude below that threshold equal to zero before taking an inverse Fourier transform to get back to the (denoised) time-domain signal.

By the way, GuassianNoise = np.random.rand(len(FourierTransformation)) produces uniformly distributed noise. To get Gaussian noise, you need randn instead of rand.


How to produce a high-pass filter from a low-pass one?

The result will indeed be a high pass filter. From your difference equation, the transfer function of the low pass filter is

$$H_l(z)=\frac{\beta}{1-(1-\beta)z^{-1}}\tag{1}$$

with $\beta=1/\alpha$.
Note that this is actually a leaky integrator, not a classic low pass filter, because its frequency response does not have a zero at Nyquist.

The high pass filter has the transfer function

$$H_h(z)=1-H_l(z)=(1-\beta)\frac{1-z^{-1}}{1-(1-\beta)z^{-1}}\tag{2}$$

which has a zero at DC (i.e. at $z=1$), as it should. However, it is scaled, i.e. its value at Nyquist ($z=-1$) is not $1$, but $2(1-\beta)/(2-\beta)$. So in order to have a 0 dB gain at Nyquist, you need to scale your high pass filter with the inverse of that factor:

$$H'_h(z)=\frac{2-\beta}{2(1-\beta)}H_h(z)=\left(1-\frac{\beta}{2}\right)\frac{1-z^{-1}}{1-(1-\beta)z^{-1}}\tag{3}$$

EDIT:

You can compute the 3dB cut-off frequencies of the low pass filter (1) and of the high pass filter (2) by solving

$$|H_l(e^{j\omega_c})|^2=\frac12\quad\text{and}\quad |H'_h(e^{j\omega_c})|^2=\frac12$$

After some algebra this gives for the low pass filter

$$\cos\omega_c=1-\frac{\beta^2}{2(1-\beta)}$$

and for the high pass filter

$$\cos\omega_c=\frac{1}{1+\frac{\beta^2}{2(1-\beta)}}$$

where $\omega$ is the normalized frequency in radians:

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

with the sampling frequency $f_s$.


Determining sampling rate for signal reconstruction

Sampling a signal requires that it should be bandlimited. When it's not, you should apply a preprocessing analog bandlimiting (aka anti-aliasing) filter to enable its proper sampling. There are various concerns on how to choose the most practical anti-aliasing filter. The most important concern is the filter cutoff frequency which determines the minimum sampling rate to record the signal.

Your application is about capturing capacitor charge/dischage signals which are of spiky nature and hence posses very large bandwidths; quite large sampling rates might be necessary. In very rough terms, the dischaging signal's bandwidth is inversely proportional to the time-constant of the discharging network, assuming a first order simplified model.

If you know this time constant, then you should select the lowpass filter whose time constant is smaller than that. This will make sure that the anti-aliasing filter's bandwidth will be larger than the signal of interest's bandwidth.

So you should firt determine the time cosntat of your observed network, to predict a first order estimate of the minimum sampling rate required.


Samples per second in a row vector in Matlab?

Duration is just the length of the signal. half makes the duration half, and double makes it double!
Pitch is the auditory sensation of the frequency. If you double the length of a signal (given the rate of playback remains the rate of sampling) then it will sound like a lower frequency. Similarly, when you make it half, it will sound like a higher frequency.

The code you can use to hear it:

fs = 8000;  % sampling frequency
duration = 1; % seconds
f0=100; % frequency
t=0:1/fs:duration; %time axis
wave = 10*sin(2*pi*f0*t);
sound(wave,fs) % original tone
pause
half_wave = wave(1:2:end);
sound(half_wave,fs) % half duraton, double pitch
pause
double_wave = interp1(linspace(0,1,numel(wave)),wave ,linspace(0,1,numel(wave)*2-1));
sound(double_wave,fs) % double duration, half pitch


Amplitude spectrum (transfer function) of signal?

HINT:
If $X(j\omega)$ is the Fourier transform of a signal $x(t)$, try to remember (or look up) how $X(j\omega)$ is related to the Fourier transform of the derivative $x'(t)$.
Some people are more familiar with this in the Laplace transform domain: the time derivative corresponds to which operation in the Laplace transform domain? Replace $s$ with $j\omega$, and you'll understand the given figure.


Zero-phase vs Linear-phase vs Nonlinear-phase System

The phase seems to be either 0 or $\pi$, which correspond to real positive and real negative frequency responses respectively.

So this can be considered a real (zero phase) system.


It means that $x(n)$ and $y(n)$ can be complex numbers rather than just real numbers. As a result there is no conjugate symmetry in the Frequency domain for these signals. Note that part (c) the asterisk in $x^*(n)$ usually means to take the complex conjugate. If $x(n)$ and $y(n)$ were real, i.e. not complex, then (b) and (c) would be equivalent.


Why remove energy at Nyquist before ifft?

you don't have to set $X\left(\frac{N}{2} \right)=0$ if you don't want to.  it will correspond to this component:

$$ X(k)\frac{1}{N}e^{j 2 \pi \frac{nk}{N}}\bigg|_{k=\frac{N}{2}} = X\left(\frac{N}{2} \right)\frac{1}{N}(-1)^n $$

but when you sample some $x[n]$, FFT it and find that $X\left(\frac{N}{2} \right) \ne 0$, you do not know the phase of that Nyquist component before it was sampled.  the Nyquist component actually aliases with itself in such a way that the quadrature part of the negative frequency cancels the quadrature part of the positive Nyquist frequency component.  such would look like:

$$ = \frac{X\left(\frac{N}{2} \right)}{N}(-1)^n = \frac{X\left(\frac{N}{2} \right)}{N}\left( \cos(\pi n) + A \sin(\pi n) \right) $$

$A$ could be anything and you wouldn't know it because $\sin(\pi n)$ is always zero for any integer $n$.

so if you're doing audio, i would likely set both the Nyquist and the DC components to zero and also insure the complex-conjugate symmetry with the other bins before iFFT.  just to make sure that, within rounding error, the imaginary part of your result is zero (the result is purely real) and that you have an unambiguous magnitude and phase for every non-zero component.


What's the essential bandwidth of the unit step function?

there can be no essential bandwidth of the unit step function.  defining such would lead you to self-contradiction.

suppose you have a function $x(t)$ which has Fourier Transform

$$ X(f) = \int\limits_{-\infty}^{+\infty} x(t) e^{-j 2 \pi f t} \ dt $$

and from $X(f)$, you defined some consistent measure of bandwidth $B_x$ in such a way that is independent of the amplitude of $X$ and is dimensionally consistent with $f$.

now suppose you scale the time argument of $x(t)$ to speed it up or slow it down:

$$ y(t) = x(a t) \quad \text{for } a > 0 $$

we know that the the Fourier Transform of $y(t)$ would be

$$ Y(f) = \frac{1}{a} X\left(\frac{f}{a} \right) $$

and the bandwidth of $Y(f)$ defined in the same manner as $B_x$ would be

$$ B_y = a B_x $$

so, suppose $a>1$, then you speed up $x(t)$ by the factor $a$ and the bandwidth is increased by the same factor.

but if you "speed up" (or slow down) the unit step function, you still have the very same unit step function:

$$ u(t) = u(a t) \quad \text{for } a > 0 $$

where

$$ u(t) \triangleq \begin{cases}
1, & \text{if } t > 0 \\
0, & \text{if } t < 0 \\
\end{cases} $$

so scaling time for the unit step does not change it, yet if it had a bandwidth, scaling the time must change that bandwidth.  a contradiction results.


How do I retrieve texture using GLCM and classify using SVM Classifier?

Are you using Matlab?  If so then you would either need the Bioinformatics Toolbox, that includes an SVM classifier, or you can download libsvm, which has Matlab wrappers for training and testing.

Then you are going to need some labeled data. Are you classifying liver tummors as opposed to healthy liver?  Then you would need images of liver tumors and healthy livers, each labeled as such.

Then you need to compute some features. What those are, depends on the nature of the problem.  Texture features seem like a good start.  Consider using co-occurrence matrices or local binary patterns.

Edit:
As of the R2014a release there is a fitcsvm function in the Statistics and Machine Learning Toolbox for training a binary SVM classifier. There is also fitcecoc for training a multi-class SVM.


Constructing a tone through additive synthesis

You say:

  I'm not after any particular timbre, certainly I'm not try to emulate a
       piano, I'm simply after a harmonic tone that is pleasing to the ear.

but I think something that is pleasing to the ear is probably harder to quantify than emulate a piano!

Bill Sethares has done lots of good work (both theoretically in signal processing and practically as a musician).  Bill wrote a book specifically about these ideas.

You can hear an example of his work, called C to shining C here. There is a discussion about it on Wikipedia.

The central idea is that consonance and dissonance are related to psychoacoustics: how the brain perceives sound.  As Bill says:

  For instance, the perception of "timbre" is closely related to (but also quite distinct from) the physical notion of the spectrum of a sound. Similarly, the perception of "in-tuneness" parallels the measurable idea of sensory consonance. The key idea is that consonance and dissonance are not inherent qualities of intervals, but are dependent on the spectrum, timbre, or tonal quality of the sound. To demonstrate this, the first audio track on the accompanying CD plays an example where the octave has been made dissonant by devious choice of timbre, even though other, non-octave intervals remain consonant. In fact, almost any interval can be made dissonant or consonant by proper sculpting of the timbre.

So, now that I've said it's too hard to answer your question, I'll try to give some pointers.

The main parameters to look at are the amplitudes ($A_k$), decay parameter ($\alpha_k$), (relative) phases ($\phi_k$) and frequency ($\omega_k$) relationships between your tones:

$$
s(t) = \sum_{k=0}^{K-1} A_k e^{\alpha_k t} \sin( \omega_k t + \phi_k)
$$

For a harmonic sound, you probably want to select $\omega_k = (k+1) \omega$ where $\omega$ is your fundamental frequency.

For "brighter" sounds, you want $A_k$ for $k>1$ to be larger (i.e. the higher harmonics have greater amplitude).

You can do flanging or phasing by modulating the phase:
$$
\phi_k(t) = \Phi_k \sin(\omega_\phi t)
$$
where $\omega_\phi$ is much lower frequency than $\omega_0$.


Bilinear interpolation implemented by convolution

Honestly I haven't read that article you linked to, but as long as you want a convolution kernel for 2D bilinear interpolation, then the following should help.

Bilinear interpolation gives a crude result which can be sufficient in case the application does not require a perfect output otherwise. Then the following simple triangular 2D convolution kernel implements a bilinear interpolation of the ratio (after compression by $M$) given by $r = \frac{L}{M} > 1$.

L = 9;          % upsampling ratio
M = 5;          % downsampling ratio (M < L assumed)

h   = (1/L)*conv(ones(1,L),ones(1,L));  % 1D linear interpolator.
hBL = h'*h;                             % 2D bi-linear interpolator.

I = im2double(imread('Cameraman.tif')); % Read some image into double
S = size(I);                            % size of image

Ie = zeros( L*S(1), L*S(2) ) ;          % expand the image by stuffing zeros
Ie( 1:L:end , 1:L:end ) = I;            % assign the original pixels

Iup = conv2(Ie,hBL);                    % linear upsample (interpolate) by L
Iint = Iup(1:M:end,1:M:end);            % downsample by M

figure,imshow(I);title('original image');
figure,imshow(Iint);title('upsampled by L/M times');


Advantage of using nonideal filter in purely digital application

Your assumption #2 is false.  Zeroing bins does not produce an ideal filter with unity gain in the passband.  Far from it if you look closely (between FFT bins).


What could be causing these humps every 10 Hz on my PSD estimate data?

These are indeed artefacts from the stimulation

x = zeros(1,10*4000)
x(1) = 1;
x(401) = 1;
[power,f]=periodogram(x,hamming(length(x)),[],4000);
scatter(f,log(power),1);
xlim([0 200]);


Everything you see are numerical artefacts. A "low pass" filter with Nyquist as cut-off frequency is simply a connection. However, you're asking for a recursive filter of order $2$ or $3$, so what you get is a filter with identical numerator and denominator polynomials (i.e., pole/zero cancellation). Check your filter coefficients, a and b should be identical.

So what you have is an unnecessarily complicated implementation of the identity operation, which results in numerical errors. The numerical errors are caused by the "filtering" operation (made worse by the double filtering implemented in filtfilt), and by the combination of FFT and IFFT.

Also note that you should plot the absolute value of spec_w, because due to numerical inaccuracies spec_w might have a small imaginary part and I don't know what your plot function does with a complex-valued vector.


Why is Error Vector Magnitude (EVM) measured as an RMS value?

You've got it – zero mean noise still has power. That power happens to be its variance. RMS of such signals thus happens to be the noise amplitude's standard deviation – and give a sensible number to assess link quality and error probabilities.


Online boosting (training on the fly) for face detection

If you are asking about online learning using OpenCV's Cascade Classifier, I'm pretty sure it is not possible. The reason is that Cascade Classifier consists of number of simpler classifiers which are decision trees in case of OpenCV's face detection system. I don't know exactly, what decision tree algorithm is used, but as far as I know there are no algorithms supporting online learning (such an algorithm would incorporate tree restructuring which is quite inefficient). So, if the base "bricks" do not support online learning, the whole "building" will not support it either.

On other hand, if you ask whether it is possible in general, then the answer is "yes". All you need is to find algorithm that supports online learning (such as artificial neural networks) and adapt the whole face recognition system to support it. You may still use cascade structure and haar-like (or even newer local binary pattern) features, but base classifiers should be replaced.

However, in practise I would say "don't do it". Or at least search for simpler options. It is quite a lot of work, both - research and development.


Invariances of FFT-based Image-Registration vs. SIFT-Features

I don't understand whether these processes are also invariant to object-alterations!

They are not.

How would the extracted fft-features look, if I alter the object (scratches, marks, dents etc.)? Would I still be able to match them properly?

Yes. But within certain limits.
The Fourier-Mellin transform is based on the Discrete Fourier Transform (DFT) whose coefficients depend directly on its input image's values. Any difference in some $I(x,y)$ will be propagated to $\mathcal{F}(I(x,y))$ and subsequent stages.
In the Fourier Mellin transform, this propagation is even longer because the transform's output is the normalised Discrete Fourier Transform (DFT) of the cartesian-to-polar followed by log transformed amplitude DFT. That is, two DFTs, one on the image and one on circular paths around the amplitude spectrum of the image. But the end result is the same. Changes in $I(x,y)$ propagate to the transform's output. Both transforms are inversible by the way.
In the context of your application, this would introduce some "fuzziness" around the centroid of your classes (For example "Part is good", "Part is defective"), which would be expected. If the texture differences (dents, scratches, etc) are in the very high frequencies it would not even matter for techniques like the one described in the Ishiyama paper where the match is based on part of the spectrum (low frequencies, or, the overall shape of the part).

The extracted features should be invariant against translation, rotation, illumination and partial object-alterations, because one use-case could be failure analysis.
...
Another idea would be to extend the existing feature-matching with Bag of Features to speed up the matching process. Would this be more suitable? Is there another, "better" approach that I haven't considered yet?

It is impossible to achieve all of these goals just at the level of the features.
First of all, what you are dealing with here is probably a One class classifier. You know what "good" looks like and you want to allow parts that are within very close proximity to this estimation of good.
The second point to clarify is the nature of the images. How do defects present themselves on the image? What does a deffective part look like? What does a good part look like? What are the processes that cause these defects? Are you interested in a "yes/no" decision from an image of a part or are you interested in localising where the defect is on the image?
To give you an example of the kind of reasoning you might want to use here: Suppose that we want to quantify how flat a surface that is produced by a planing machine is. We can illuminate it with a source at a very shallow angle that would cause any dents or pits or scratches, to appear as dark streaks or spots on the surface image. Then, simply by thresholding the image and measuring how many "dark pixels" appear in it gives us an estimate of how smooth that surface is. Notice here, there was no reason for fancy features, etc. Only knowledge of the process and the output we are after.
Finally, illumination differences in an industrial environment are easily controlled. And if you have a challenging subject like metalic parts that can produce reflections, you definitely want to control the lighting. You can apply some sort of equalisation or overall normalisation of the brightness, but it is impossible to control for ALL of the artifacts that reflections might produce (epsecially from a single image).
Hope this helps.


How can I calculate the frequency using different color bands?

From a practical perspective, there's good news and there's bad news.  The good news is that, in a sense, consumer grade digital cameras [at least] already do something like what you're describing, so there's a large amount of information available about the process, which is called demosaicing.  The bad news is that, in "standard" designs, the various RAW color channels actually have lower resolution than the physical sensor array, and the color channels in non-RAW images have already been reconstructed.

The reason for this is that the physical sensors only detect light intensity (not color), and therefore each available physical "pixel" is paired with a color filter, according to some pattern. For example, the Bayer Filter is a common pattern which allocates 50% of the available spatial resolution to the green channel, but only 25% each to the red and blue channels.  The various spatial patterns are called Color Filter Arrays (CFAs).  In effect they model the idea of separate RGB sensors, just in a more cost-effective manner.  If there is some way of exploiting "spectral crosstalk" from the materials used in a given CFA, it would be handled in the context of demosaicing.

That means, if you're dealing with a camera which uses that type of sensing scheme, the most reasonable thing you can hope for is access to the camera RAW data (i.e. un-demosaiced), from which you can apply a higher-quality demosaicing algorithm to reconstruct the image.  From the wiki:

  When one has access to the raw image data from a digital camera, one
  can use computer software with a variety of different demosaicing
  algorithms instead of being limited to the one built into the
  camera...The differences in rendering the finest detail (and grain
  texture) that come from the choice of demosaicing algorithm are among
  the main differences between various raw developers; often
  photographers will prefer a particular program for aesthetic reasons
  related to this effect.

The following link seems to be a good overview of various reconstruction algorithms:

A Study of Spatial Color Interpolation Algorithms for Single-Detector Digital Cameras

In higher-end cameras, however, the RAW color channels might not even be spatially interleaved.  Two of the different higher-end schemes (beam splitting, temporal interleaving) are reviewed here.

One type of sensor that I think works more like what you had in mind is the Foveon X3 Sensor. The following paragraph from the link gives a good overview of the resolution issues:

  For example, the dimensions of the photosite array in the sensor in
  the Sigma SD10 camera are 2268 × 1512, and the camera produces a
  native file size of those dimensions (times three color layers). This
  amounts to approximately 3.4 million three-color pixels. However, it
  has been advertised as a 10.2 MP camera by taking account of the fact
  that each photosite contains stacked red, green, and blue color
  sensing photodiodes, or pixel sensors (2268 × 1512 × 3). By
  comparison, the dimensions of the photosite array in the 10.2 MP Bayer
  sensor in the Nikon D200 camera are 3872 × 2592, but there is only one
  photodiode, or one pixel sensor, at each site. The cameras have equal
  numbers of photodiodes, and produce similar RAW data file sizes, but
  the Bayer filter camera produces a larger native file size via
  demosaicing.

I think this is what you're looking for, assuming I understood the question.


Convolution of an Image with a Kernel That Is a Product of Two Functions

It is equivalent to a bilateral filter. You can move the kernel to each local part of image to do the dot product, then sum them up to get the result in a loop. There is also a fast implementation with O(1) in matlab programming for your reference. And the relevant publications are:

K.N. Chaudhury, D. Sage, and M. Unser, "Fast O(1) bilateral  filtering using
trigonometric range kernels," IEEE Transactions on Image Processing, vol. 20,
no. 11, 2011.
K.N. Chaudhury, "Acceleration of the shiftable O(1) algorithm for bilateral filtering
and non-local means," arXiv:1203.5128v1.

I assume the Gaussian kernel is position shift-invariant. So normalization means the sum of elements inside the kernel equals to 1. For a 3*3 Gaussian kernel, it is

1/16  2/16  1/16
2/16  4/16  2/16
1/16  2/16  1/16

The bilateral filter can also be normalized with the K factor of integral of the product of  Gaussian kernel and intensity kernel. You can try both to observe the effect of normalization on your image data.


Exam Question - end-to-end delay from video capture to video playout

You can't start playing any B frame until the next following I or P frame has been scanned, compressed, transmitted, received, and decompressed.


Why are FIR filters still stable even though they contain poles?

FIR filters contain as many poles as they have zeros.  but all of the poles are located at the origin, $z=0$.

because all of the poles are located inside the unit circle, the FIR filter is ostensibly stable.

this is probably not the FIR filter the OP is thinking about, but there is a class of FIR filters called Truncated IIR filters (TIIR) which may have a pole on or outside the unit circle that is canceled by a zero at the same location.  the simplest example of this is the moving sum or moving average filter.  but, from an I/O perspective, these TIIR filters are FIR.

but i wouldn't naively guarantee "stability".  using control-system language, the TIIR filter are not "completely observable" and may appear stable because its impulse response appears finite in length, but inside the filter states could be going to hell, and with finite numerical precision, that internal instability will eventually show up at the output.

we have to disabuse ourselves from the notion that "FIR filters have no poles".  ain't true.


The problem appears to be that your windowing function doesn't preserve complex conjugate symmetry. So ifft(y1) has a significant imaginary part. By discarding this through the real() operation you induce a significant error which results in the discrepancy. To verify try

z = ifft(y1); plot(imag(z));


Probably the simplest OCR algorithm would be the case where the characters are from a known font, in perfect rows with no distortion or rotation, and then you can find occurrences of each character in the text by doing cross-correlation with the known characters from the known font.

Or even simpler if you can create the font before it's used, with cheaper ways to recognize characters, like the MICR fonts or OCR fonts.


FFT of white noise with different sampling rate

The variance of sampled white noise that is filtered by a brick-wall lowpass filter is reduced by a factor equal to the ratio of the filter cutoff frequency to 1/2 the sample-rate.  You start off with the same variance for the 100mhz and 10mhz cases, but the reduction factors are different due to the different sample-rates. Or to put it another way, the 100mhz-sampled signal has less power/Hz than the 10mhz sampled signal, so when you compare them over the same bandwidth, the filtered 100mhz signal has less power.


Interference Distribution Assumption in NOMA IoT Scenario

interference intensity is different at each sub-channel but varies with same distribution

You're already using the term intensity, so this sounds like you're most likely implicitly using a Poisson Process, which is a stochastic typical model for "arrival processes", e.g. packets appearing in a network. It's based on the idea that the probability of an event occurring within a time period is simply proportional to the length of that period, and does not depend on the start time of the period.

Really check your literature thrice for them indications of where the term "intensity" came from. It's totally possible that something in there confirms the above (i.e. they say that the intensity $\lambda$ defines the factor between expected number $N$ of occurences and the duration of observation $T$, $\mathbb E[N(T)] = \lambda T$); it's also totally possible that other assumptions make this model impossible ("we assume interference to happen in bursts/clusters…").

Notice how careful I am about proposing this: it is a mathematically very convenient model, and it's pretty accurate to model a lot of naturally occurring behaviours, but it can badly jump in your face if you forget it's likely a simplifying behaviour.

Take this example (which I made up right now, but I think it's realistic enough to make my point):

  Our packets are 0.1 s long. Interference emissions are much shorter, 0.001 s; they either appear within a packet (leading to a broken packet), or they don't. We model the occurrence of the interference to be a Poisson process.

This works well if the interferers are, for example, independent sensors that send a short packet when someone enters a room. In essence, your model represents the reality of when interference bursts happen accurately.

Even when the interferers are actually doing ALOHA for medium access, due to the packets being so short, the chance that you'll see them clustering is so low that you can still assume disjoint time periods (i.e. your packet transmissions) have independent interference arrival probabilities.

Now, if you happen to make your packets shorter, e.g. 0.008 s, you'll start "seeing" that ALOHA of the interferers mean that if there is indeed a correlation between interference occurrence during a packet and during a following packet.

If you made assumptions for your NOMA system based on the Poisson-ness of the process (like, you're doing CSMA/CD with deterministic backoff, you don't use an overly long spreader to distribute your info,…), you will run into a real congestion problem although your channel isn't actually heavily occupied.

So, treat your model assumptions with care. If in doubt, back them up with a quick simulation of what happens if they aren't fully met; if you can show that your overall system is robust about small model deviations, that greatly increases the credibility of whatever you're writing/designing.


Designing an IIR to generate a specific data sequence

Yes, there are time-domain design methods for IIR filters. One of the best-known ones is Prony's method. It is well described in the book Digital Filter Design by T.W. Parks and C.S.Burrus (ch. 7.5).

If the desired sequence is (right-sided) periodic, then there's a very simple and exact solution possible (which is usually an exercise for DSP students).


Transform the signal into frequency domain

$$F(f) = \frac{1}{\left|a\right|} \left[~\frac{\Pi(f-\nu_{c})}{2} + \frac{\Pi(f+\nu_{c})|}{2} ~\right]$$

As you know the fourier domain representation of the hilbert operator is

$$H(f) = -j\text{sgn}(f)$$

Applying this transform means you just multiply $H(f)$ to $F(f)$. That gets you

$$\hat{F}(f) = -j\text{sgn}(f) \cdot \frac{1}{\left|a\right|} \left[~\frac{\Pi(f-\nu_{c}}{2} + \frac{\Pi(f+\nu_{c})|}{2} ~\right] $$

$$ \hat{F}(f) = \frac{1}{\left|a\right|} \left[~\frac{\Pi(f-\nu_{c})}{2j} + \frac{-\Pi(f+\nu_{c})|}{2j} ~\right] ~~ \text{Remember} -j = \frac{1}{j}$$

which is same as

$$ \hat{F}(f) = \frac{\Pi(f)}{\left|a\right|} \ast  \left[~\frac{1}{2j} \delta(f-\nu) + \frac{-1}{2j} \delta(f+\nu) ~\right] $$

Inverse Fourier Transformation of this gives you

$$\hat{f}(t) = \text{sinc}(at)\sin{\left(2\pi\nu_{c}t\right)}$$

Hope this helps!


You should photograph a known calibration target with your system, (like this one or the one that I put in the thumbnail).

This will allow you discovering the OECF of your system. You can either calculate it by yourself, given the technical data sheet or use third-party software like Imatest, ImageAnalyzer or others.


Writing a Discrete Fourier Transform program

When applied correctly, the DFT yields exactly the same result as convolution, it is no approximation (neglecting rounding errors due to limited precision of floating point calculations, of course).

To calculate the circular convolution of vectors $x$ and $y$ (both of length $N$):
$$
x \circledast y = \operatorname{IDFT}_N\left[\operatorname{DFT}_N[x]\cdot\operatorname{DFT}_N[y]\right]
$$
To calculate the linear convolution first append $N-1$ zeros to both vectors yielding the zero-padded vectors $\tilde x$ and $\tilde y$. Then:
$$
x * y = \operatorname{IDFT}_{2N-1}\left[\operatorname{DFT}_{2N-1}[\tilde x]\cdot\operatorname{DFT}_{2N-1}[\tilde y]\right]
$$

To implement the [I]DFT I highly recommend the FFTW library. I promise that it will take you less time to learn its usage than to write your own implementation. For the special case of $x$ and $y$ containing only $1$ and $0$ their might be a more efficient algorithm, though. You would have to look into the Cooley-Tukey algorithm that is used for the FFT and see if you can simplify it for your needs. Whether it gets simpler will also depend on how you define computational complexity - if you define it in terms of additions and multiplications a custom FFT implementation probably won't buy you much. But if you go down to the hardware implementation level you can save some multipliers in the first butterfly stage by replacing them with multiplexers. Anyway, it is well known that the FFT algorithm complexity is of order $N\log N$ as you say. I would like to think that you don't have to show this again.

If you'd like to play around with IFFT and FFT I recommend Octave (an open source Matlab clone).


How to denoise a conversation between two people?

I think what you're attempting to do is a case of blind source separation (check Wikipedia). There are many techniques to solve this kind of problem, maybe you could check those on the Wikipedia and use the one that fits your needs best.

Particularly, last year one of my teachers showed us a demo of the ICA algorithm, and we got to play round with it in MATLAB. It seemed really effective, and there's a research group that even created this MATLAB toolbox for the ICA algorithm. Here you can check a simple demo (called cocktail party problem) of how ICA works and how effective it is.

Since you say you're new to DSP, I'm not sure if ICA or the other techniques would be too advanced for you, but maybe this will put you in the right direction. Good luck.


The "best" detector is a highly subjective subject, and there is (in my opinion) not a definitive answer here. I work in radar processing and I've used everything that you've mentioned in one form or another.  All of these CFAR methods are tools, and all of them have good/bad applications.  For instance, asking what the "best" wrench is in the toolbox highly depends on what you're going to use that wrench for.  This is no different.  OS CFAR can help with non-stationary noise floors which are varying quite a bit from CUT (cell under test) to CUT, but sometimes CA-CFAR can work just fine and is much simpler to implement.

In summary, they're all tools, and therefore no single CFAR detector is inherently "better" than any other one; it all depends on your system, the noise floor, the interference environment, the clutter environment, etc.


The filtering distribution $p(x_k|y_1,...,y_k)$ is computed via the Kalman filter and given a linear-time invariant system is parametrized by the mean $E[x_k|y_1,...,y_k]$ and covariance matrix.

When doing smoothing we want to compute $p(x_k|y_1,...,y_N)$ with $N \ge k$, in this case the mean we compute is $E[x_k|y_1,...,y_N]$.


How one apply correctly FFT in image denoising

The trick inside the paper is the following:

What you want to compute is $\sum_{i \in W} |I(x+i)-I(y+i)|^2$, where $I$ is an image, $x$ and $y$ two noisy pixels and $i$ is a 2D offset used to define a patch.
Expanding the expression yields: $\sum_i I^2(x+i) + \sum_i I^2(y+i) - 2 \sum_i I(x+i)I(y+i) = A + B - 2C$.
$A$ and $B$ are computed using a squared integral image, i.e., an integral image from the squared original image.
$C$ is the convolution between the two patches centered on $x$ and $y$. Thus, it can be computed in the Fourier domain, where it becomes a multiplication. You get the value of $C$ by computing the Fourier transform of the patch around $x$, the patch around $y$, pointwise-multiplying these results and taking the inverse Fourier transform of the multplication result.

The Fourier transform is obviously a 2D transform since you are working with 2D data. What you obtain for a given patch is a 2D array of complex values.

Additional notes

In my opinion this article is not the best NL-means speedup strategy.
Experiments I did way back in 2007/2008 show that pre-selection of patches are better (both in terms of speed and quality of the results).
I have started blogging about these here, but unfortunately I am looking for time to finish the posts.

The original NL-means papers mention blockwise implementations that can be interesting. There are fundamentally 2 ways in implementing NL-means:

writing a denoising loop for every pixel in the image
writing a denoising loop for each patch, then back-project the patches to form an image.

The first impolementation is the original approach, because in 2005 memory and multicore CPUs were expensive. I chose on the other hand number 2 on recent hardware in the past 2 years. It depends on your typical image size and if you want to be able to compute domain transforms like DFT/DCT (as in the proposed paper and in BM3D).


Wiener Filter Additive Noise Uncorrelated

Let me change the notation slightly to make things clear. Let $d[n]$ be the desired signal, $x[n]$ is the noisy input signal, and $v[n]$ is the noise. Furthermore, we assume that

$$x[n]=d[n]+v[n]\tag{1}$$

and that the desired signal $d[n]$ and the noise $v[n]$ are uncorrelated and at least one of them has zero mean, which implies that their cross-correlation equals zero:

$$R_{dv}[n]=E\{d^*[k]v[k+n]\}=0\tag{2}$$

From $(2)$ it follows that the cross-correlation between the desired signal $d[n]$ and the input signal $x[n]$ is given by

$$\begin{align}R_{dx}[n]&=E\{d^*[k]x[k+n]\}\\&=E\{d^*[k]\left(d[k+n]+v[k+n]\right)\}\\&=E\{d^*[k]d[k+n]\}+E\{d^*[k]v[k+n]\}\\&=E\{d^*[k]d[k+n]\}=R_{d}[n]\tag{3}\end{align}$$

From $(3)$ it follows that the cross-power spectral density of $d[n]$ and $x[n]$ equals the power spectral density of $d[n]$:

$$S_{dx}(\omega)=S_{d}(\omega)\tag{4}$$

With the given assumptions, the frequency response of the non-causal Wiener filter can be written as

$$H(\omega)=\frac{S_{dx}(\omega)}{S_x(\omega)}=\frac{S_d(\omega)}{S_d(\omega)+S_v(\omega)}\tag{5}$$


Modelling zero padded OFDM as circular convolution between channel and FFT frame

The answer lies in relationship between Linear Convolution and Circular Convolution. Even Zero padding is not required to model circular convolution using linear convolution. Zero-padding is only required when we are sending OFDM symbols consecutively to cope up with ISI. My point will get clearer in a minute. I will try to show mathematically, by taking an example, where $\vec{d} = d[0], d[1], d[2], d[3], d[4]$ are time-domain samples, meaning IDFT of 5 data symbols in context of OFDM. And the channel is 3 Taps so, $\vec{h} = h[0], h[1], h[2]$. Now, the channel is an LTI system so it takes linear convolution of input $\vec{d}$ with channel impulse response $\vec{h}$, which can be given by the following equation:$$y_{lin}[n] = \sum^{m=2}_{m=0}h[m]d[n-m],$$ignore noise for the sake of this question. Linear convolution output $\vec{y}$ will be of length $5+3-1=7$. Expand the convolution in matrix form to get the following: $$\vec{y_{lin}} = \begin{pmatrix}y_{lin}[0]\\y_{lin}[1]\\y_{lin}[2]\\y_{lin}[3]\\y_{lin}[4]\\y_{lin}[5]\\y_{lin}[6]\end{pmatrix} = \begin{pmatrix}h[0]&0&0&0&0\\h[1]&h[0]&0&0&0\\h[2]&h[1]&h[0]&0&0\\0&h[2]&h[1]&h[0]&0\\0&0&h[2]&h[1]&h[0]\\0&0&0&h[2]&h[1]\\0&0&0&0&h[2]\end{pmatrix}.\begin{pmatrix}d[0]\\d[1]\\d[2]\\d[3]\\d[4]\end{pmatrix}$$This is not circular convolution because the $H$ matrix is NOT Circulant yet. Now remove last 2 rows and add them to top 2 rows to get the following: $$\vec{y_{circ}} = \begin{pmatrix}y_{circ}[0]\\y_{circ}[1]\\y_{lin}[2]\\y_{lin}[3]\\y_{lin}[4]\end{pmatrix} = \begin{pmatrix}h[0]&0&0&h[2]&h[1]\\h[1]&h[0]&0&0&h[2]\\h[2]&h[1]&h[0]&0&0\\0&h[2]&h[1]&h[0]&0\\0&0&h[2]&h[1]&h[0]\end{pmatrix}.\begin{pmatrix}d[0]\\d[1]\\d[2]\\d[3]\\d[4]\end{pmatrix},$$Notice that the last 3 values of $\vec{y_{circ}}$ are unchanged from $\vec{y_{lin}}$, only the first 2 values are changing. And, now the new $H_{circ}$ matrix becomes Circulant. So, the trick lies in removing last $(L-1)$ values of received samples, $y[N], y[N+1], y[N+2], ...,y[N+L-2]$, and adding them to top $(L-1)$ received values, $y[0], y[1], y[2], ..., y[L-2]$, before taking $DFT$ of this new vector $\vec{y}_{circ}$ obtained from original received vector $\vec{y}_{lin}$. Now in context of OFDM, Zero-Padding is required otherwise the last $(L-1)$ samples of this OFDM symbol will get added to the first $(L-1)$ samples of next symbol, and then we will not be able to do the trick which allows us to convert linear convolution into circular convolution. Now, like conventional OFDM RX Chain, DFT can be done on $\vec{y_{circ}}$ to get the following:$$DFT\{ \vec{y_{circ}} \} = H[k].D[k], $$ where $H[k] = DFT\{ \vec{h} \}$. So, basically, the above described $(L-1)$ Zero-Padded system is equivalent to adding Cyclic Prefix in absence of Noise, if we do that extra processing on received samples $y[n]$, in the sense it makes the channel single tap in frequency domain. Disadvantages: First, because we are going to add last $(L-1)$ recieved samples to first $(L-1)$ samples of received vector, hence, in presence of AWGN the noise variance for the first $(L-1)$ data samples will be doubled (Because 2 samples of uncorrelated noise will get added too). However, in absence of noise, both Zero-padding and Cyclic-Prefix will perform same. Second, Zero-padding is anyway required so we are going to waste channel utilisation, but practically speaking, sudden shutdown of TX modulator is not possible. There will be transient and that will ruin the time-domain samples at transmission. And, since TX modulators cannot support complete shutdown of power, zero-padding based OFDM systems are very rarely implemented.


There is some confusion here. Frequency Response and Group Delay are properties of a System not the properties of a signal. You can only calculate the spectrum of a signal (magnitude and phase). A "signal" doesn't have group delay unless it's the impulse response of an LTI system.


There are a few simple things that you could try:

sliding-window median filtering with an appropriate window size
linear-regression in a sliding window
robust linear regression to the total signal, if you know the signal duration etc well enough

EDIT: Further idea after a comment from OP on another answer

Since you want the slope at 160000 be the same as before (i.e. no flat section in your signal), you could use a robust linear regression technique such as RANSAC with an appropriate number of sample points to analyze/filter your data.


Absolute integrability is not a necessary condition for the Hilbert transform to exist. Simple counterexamples are the Hilbert transforms of $\cos(\omega t)$ and $\sin(\omega t)$, given by $\sin(\omega t)$ and $-\cos(\omega t)$, respectively, if $\omega>0$.

Since the functions $u(t)$ and $\textrm{sgn}(t)$ are related by scaling and an additive constant, and since the Hilbert transform of a constant is zero, the Hilbert transforms of the two functions - if they exist - must be scaled versions of each other:

$$\textrm{sgn}(t)=2u(t)-1\\\Rightarrow \mathcal{H}\{\textrm{sgn}(t)\}=2\mathcal{H}\{u(t)\}\tag{1}$$

I believe that the Hilbert transform of these two functions can be computed in the following way. [Edit: this turned out to be true. See below and the comments to RBJ's answer.] We know that the Hilbert transform can be written as a convolution:

$$\mathcal{H}\{x(t)\}=\frac{1}{\pi t}\star x(t)\tag{2}$$

This can be rewritten as

$$\mathcal{H}\{x(t)\}=\frac{1}{\pi}\left(\log(|t|)\right)'\star x(t)=\frac{1}{\pi}\log(|t|)\star x'(t)\tag{3}$$

where $'$ denotes differentiation w.r.t. $t$. The rightmost equality in $(3)$ is true because we generally have

$$(f\star g)'=f'\star g=f\star g'\tag{4}$$

With $x(t)=u(t)$ we get $x'(t)=\delta(t)$, and from $(3)$ we obtain

$$\mathcal{H}\{u(t)\}=\frac{1}{\pi}\log(|t|)\star \delta(t)=\frac{1}{\pi}\log(|t|)\tag{5}$$

From $(1)$ and $(5)$ we get

$$\mathcal{H}\{\textrm{sgn}(t)\}=\frac{2}{\pi}\log(|t|)\tag{6}$$

The results given by Eqs $(5)$ and $(6)$ can also be derived in the frequency domain. The Fourier transform of the convolution

$$\mathcal{H}\{u(t)\}=\frac{1}{\pi t}\star u(t)$$

is given by

$$\mathcal{F}\left\{\frac{1}{\pi t}\right\}\mathcal{F}\{u(t)\}=-j\;\textrm{sgn}\{\omega\}\left(\pi\delta(\omega)+\frac{1}{j\omega}\right)=-\frac{\textrm{sgn}(\omega)}{\omega}=-\frac{1}{|\omega|}$$

So we have

$$\mathcal{H}\{u(t)\}=\mathcal{F}^{-1}\left\{-\frac{1}{|\omega|}\right\}$$

From this document (p. 27, Ex. 19.4), we know that

$$\mathcal{F}^{-1}\left\{-\frac{1}{|\omega|}\right\}=\frac{1}{\pi}\log(|t|)$$

confirming Eq. $(5)$.

(Note the different scaling in the cited document due to the use of the unitary Fourier transform.)

EDIT:

The problem can also be solved using the definition of the Hilbert transform of bounded functions:

$$\mathcal{H}\{f(t)\}=\frac{1}{\pi}\left[\mathrm{p.v.}\int_{|\tau|<1}\frac{f(\tau)}{t-\tau}d\tau+\mathrm{p.v.}\int_{|\tau|>1}f(\tau)\left(\frac{1}{t-\tau}+\frac{1}{\tau}\right)d\tau\right]$$

where $\mathrm{p.v.}$ denotes the Cauchy principal value.

For $f(t)=u(t)$ we have

$$\mathcal{H}\{u(t)\}=\frac{1}{\pi}\left[\underbrace{\mathrm{p.v.}\int_{0}^{1}\frac{d\tau}{t-\tau}}_{\displaystyle I_1}+\underbrace{\mathrm{p.v.}\int_{1}^{\infty}\left(\frac{1}{t-\tau}+\frac{1}{\tau}\right)d\tau}_{\displaystyle I_2}\right]$$

If $t$ is outside the interval $(0,1)$, the first integral is simply

$$I_1=\int_{0}^{1}\frac{d\tau}{t-\tau}=\log\left(\frac{t}{t-1}\right)$$

If $t\in (0,1)$ we get

$$I_1=\lim_{\epsilon\to 0^+}\int_0^{t-\epsilon}\frac{d\tau}{t-\tau}+\int_{t+\epsilon}^1\frac{d\tau}{t-\tau}=\log\left(\frac{t}{\epsilon}\right)+\log\left(\frac{-\epsilon}{t-1}\right)=\log\left(\frac{t}{1-t}\right)$$

Since in both cases the argument of the logarithm is positive, both results can be summarized as

$$I_1=\log\left|\frac{t}{1-t}\right|$$

If $t<1$ the second integral is

$$I_2=\lim_{\epsilon\to 0^+}\left[\int_1^{1/\epsilon}\frac{d\tau}{t-\tau}+\int_1^{1/\epsilon}\frac{d\tau}{\tau}\right]=\lim_{\epsilon\to 0^+}\left[\log\left(\frac{t-1}{t-1/\epsilon}\right)+\log\left(\frac{1}{\epsilon}\right)\right]=\\=\lim_{\epsilon\to 0^+}\log\left(\frac{t-1}{\epsilon t-1}\right)=\log\left(1-t\right)$$

If $t>1$ we have

$$I_2=\lim_{\epsilon\to 0^+}\left[\int_1^{t-\epsilon}\frac{d\tau}{t-\tau}+\int_{t+\epsilon}^{1/\epsilon}\frac{d\tau}{t-\tau}+\int_1^{1/\epsilon}\frac{d\tau}{\tau}\right]=\\=\lim_{\epsilon\to 0^+}\left[\log\left(\frac{t-1}{\epsilon}\right)+\log\left(\frac{-\epsilon}{t-1/\epsilon}\right)+\log\left(\frac{1}{\epsilon}\right)\right]=\\=\lim_{\epsilon\to 0^+}\log\left(\frac{1-t}{\epsilon t-1}\right)=\log(t-1)$$

Again, both results can be summarized as

$$I_2=\log|1-t|$$

Combining these results, we obtain the Hilbert transform of $u(t)$:

$$\mathcal{H}\{u(t)\}=\frac{1}{\pi}(I_1+I_2)=\frac{1}{\pi}\left(\log\left|\frac{t}{1-t}\right|+\log|1-t|\right)=\frac{1}{\pi}\log|t|$$

which agrees with our previous result given by Eq. $(5)$.


Only the first criterion is correct when referring to input-output stability. The second criterion is just a way to compute the poles, as mentioned in a comment. All realizable continuous-time systems using lumped elements have poles, which doesn't contradict the stability requirement.
Poles in the right half-plane (RHP) make a causal system unstable. Note that if we let go of the causality requirement, RHP poles don't cause instability but they make the system non-causal.
Note that poles on the imaginary axis cause transients that don't die out, which is undesirable in almost all practical applications. So in practice you want all poles to lie on the open left half-plane, i.e., $\textrm{Re}\{s_{\infty,k}\}<0$, where $s_{\infty,k}$ is a pole of the system.


limit of the integral in the anti-transform

I think you are making this more complicated than it needs to be.

The Fourier and autocorrelation integrals are defined from $-\infty$ to $+\infty$. That's always correct and the saftest way to write them.

If you are integrating over a function with limited support on $[0,T]$ than you can adjust the integration interval to  $[0,T]$ too. The areas outside the support interval is all zero so it makes no difference whether you include them in the integration or not. In this case both versions are identical and both versions are correct. You can use whichever one is more practical for your specific problem or application.

In your specific example: if $u(t)$ has support on $[0,T]$  then it's autocorrelation $r_{uu}(t)$ also has support only on $[0,T]$. That means you can calculate the PSD  by executing the Fourier integral over either $[-\infty,+\infty]$ or $[0,T]$. Both are correct and both will give you the same result.

$S(f)$ has infinite support so you need to use $[-\infty,+\infty]$ for the inverse Fourier integral.

$$S_{uu}(\omega)=\int_{-\infty}^{+\infty}r_{uu}(t) \cdot e^{-j \omega t} d \omega = \int_{0}^{T}r_{uu}(t) \cdot e^{-j \omega t} d \omega$$

$$r_{uu}(t)= \frac{1}{2 \pi} \int_{-\infty}^{+\infty} S_{uu}(\omega) \cdot e^{+j \omega t} d \omega $$


Proof of time-invariance of continuous-time system

Let $y_1(t)$ be the response to the signal $x_1(t)$:

$$y_1(t)=x_1(-t)\tag{1}$$

Now let $x_2(t)$ be a shifted version of $x_1(t)$:

$$x_2(t)=x_1(t-T)\tag{2}$$

The response to $x_2(t)$ is

$$y_2(t)=x_2(-t)=x_1(-t-T)\tag{3}$$

If the system were time-invariant, its response to $x_2(t)$ should be a shifted version of its response to $x_1(t)$:

$$y_2(t)=y_1(t-T)\tag{4}$$

However, from (1) we have $y_1(t-T)=x_1(-(t-T))=x_1(-t+T)$. Comparing this to (3) we see that (4) is not satisfied, and, consequently, the system is not time-invariant. Now that you've seen the proof, try drawing all the signals in order to gain a better understanding.


What is QEEG technically, and how is it computed?

qEEG or quantitative EEG is an umbrella term that encompasses many different signal processing methods that compute a number/statistic/quantitative biomarker from EEG signal traces. Algorithms usually boil down to some form of spectral analysis either using Fourier transforms, or wavelets.

The images in shown in the question are most likely generated through spectral analysis of individual EEG channels. In your figure, it seems like the author chose Delta to mean 1-3.5Hz. Some authors use Delta as all frequencies <4Hz. (https://en.wikipedia.org/wiki/Electroencephalography)

To compute the contribution of each band, one resorts to standard power spectrum analysis techniques. For example, for a given EEG channel $i$, the Delta qEEG biomarker value can be computed as:

$$
\mbox{Delta}_i = \frac{\int_{1}^{3.5} P_{i}(f) df}{\int_{0}^{\infty}P_i(f)}
$$

where $P_i(f)$ is the power spectral density of the $i^{th}$ EEG channel signal. In practice, you can use, say, Welch method to estimate $P_i(f)$ from sampled EEG channel traces, the integrals will get replaced by sum over frequency bins, and you'd only go up to Nyquist.

Repeat this process for each frequency band, for each EEG channel to get sequences $\{\mbox{Delta}_i\}_{i=1}^N$,$\{\mbox{Theta}_i\}_{i=1}^N$, $\{\mbox{Alpha}_i\}_{i=1}^N, \cdots$, where $N$ is the number of EEG channels.

The final step is to map these numbers on the image of the scalp. This is a visualization step which involves interpolation. You can use a tool like this: https://martinos.org/mne/stable/generated/mne.viz.plot_topomap.html#mne.viz.plot_topomap or https://www.nbtwiki.net/doku.php?id=tutorial:plot_a_time-frequency_representation_and_power_spectrum_of_a_single_channel#.WJtLMxC2pJI.


Respiration (breathing) signal movement artifact removal

I figured out how to tackle the problem.

I used the EMD algorithm and then smoothed the signal. I removed the first two IMF[i] and the last IMF[i] and reconstructed the signal again. This removes posture movement.
I then  identified 3sd away outliers and chopped off the signal to pieces. For each piece I calculated the FFT and found the dominant frequency.


Spatial domain vs. Frequency domain filtering of an image. Which one is better?

OK.  In theory, convolution in the time domain is the same as multiplication in the frequency domain.  And filtering a signal with a linear filter is mathematically the same as convolving by its impulse response.
So in the time domain, if you have an image $x_{n, m}$, and a filter with an impulse response $h_{n, m}$, then after filtering $x$ with $h$, you get a signal $y = x * h$.  (Note that especially in image processing, a filter's impulse response is often called a kernel.  Until you're smarter about this than I usually am, just treat "impulse response" and "kernel" as synonyms).
If you let $X = \mathcal F \left \lbrace x \right \rbrace$ denote the Fourier transform, then if $y = x * h$, $Y = H X$.  This turns out to be a happy result from 1-dimensional Fourier analysis, and an even happier result is that if you're using Cartesian coordinates, it applies to 2-dimensional (on 3- or 4-, or 42-dimensional) Fourier analysis, too.
That's the theory.  Basically, it's what you've found in your research*.
So if you have $h$ then you can compute $y$ as $y = x * h$ (or $y = h * x$ -- convolution commutes).  But you can also compute $y = \mathcal F^{-1} \left \lbrace \mathcal F \left \lbrace \mathcal h \right \rbrace F \left \lbrace \mathcal x \right \rbrace \right \rbrace $.  This seems to be dreadfully counter-productive, but for two things.
First, you can pre-compute $H = F \left \lbrace \mathcal h \right \rbrace$.  If it's settled enough, you can even compute it during design time and just have a pre-defined array in code.
Second, there's this thing called the fast Fourier transform (FFT).  From a mathematical perspective, a fast Fourier transform would be just an overly complicated way of doing a plain old discrete Fourier transform except that for large images it uses far fewer operations (hence -- Fast Fourier Transform).
Convolution of an $N \times M$ image by a $N_h \times M_h$ kernel takes -- if I'm getting my math right -- $N M N_h M_h$ multiply-add operations.  Basically, it's the area of the kernel times the area of the image.  Taking the FFT and inverse FFT of that same image area takes $k N M \log(N M)$ multiply-add operations**.  So for a large enough filter kernel, using the FFT can be faster -- particularly because you can use a technique called "overlap-add" to speed things up.
So, theoretically there is no difference between doing your filtering using convolution in the spatial-domain or by multiplying in the frequency domain.
Practically, doing it in the frequency domain is always more complicated, but it can be faster.  Also practically, when you use the FFT there can be some slight loss of numerical precision -- basically, the speed comes at the cost of a bit of extra numerical scrud due to rounding errors and whatnot.  This numerical effect is almost certainly of no concern to you.
Finally -- practically --

it's rare that you have useful filter kernels that are large enough that doing the computation in the frequency domain will speed things up.
Moreover, if all you're interested in is Gaussian filtering, there are some fast almost-Gaussian filtering technique that apply a small (3x3, I think) kernel to an image repeatedly.  Because of the Central Limit Theorem, if you choose a small kernel correctly and filter by it repeatedly, the resulting filtering will tend to approximate filtering by a Gaussian.

* Things may have been a bit clouded by the fact that there's a result that's specific to Gaussian filters, in that if $h_{n, m}$ is Gaussian, then so is $H(\omega_1, \omega_2)$ -- it's just that in general it has a different scaling.  This is not always the case, but if you're looking at some page that doesn't fully explain it, the notion of multiplying by a Gaussian in the frequency domain being the same as convolving by a Gaussian in the spatial domain would look extra-weird.
** I don't remember the exact number, but I think for a naïve FFT implementation, $k=6$, and there's more clever implementations that get it smaller yet.


Obtaining real numbers from FFT with same length as original signal

Since you are reproducing the paper's results, read it carefully.
At the bottom of page 4 it says

This representation is formed by complex numbers, eliminating the imaginary part of each number in the frequency-domain signal. For this transformation, it is needed to calculate the power spectral density (PSD), as shown in Equation (3),
$$
P = \lim_{T\to\infty} \frac{1}{T} \int_0^T |x(k)|^2 dt \tag{3}
$$
where $P$ represents the energy from the signal, $T$ represents the length of the signal lapse and $x(k)$ represents the frequency-domain signal. The spectrum is normalized by the length of the signal.

The PSD is a real function with respect to frequency $k$, so it has the same length of FFT result $x(k)$ if you keepace with Eq. (3).
There are many methods to estimate PSD, but if you want to recreate the results in the paper, you should use exactly the same equation.


Granulometry by morphological opening and closing

When you use the opening, the overall operation is called a granulometry. However, when you use a closing, it's called an anti-granulometry. It's the dual operation. Combined together, they form the Pattern Spectrum.

In your case, if you apply a closing, you will certainly measure the structures between the marbles.

However, it is definitely not recommended to apply a smoothing first, because it will smooth the edges and small particles. Consequently, the particles detection will be less accurate.


Phase differences between DFT and FFT (in MatLab)

you should replace the line

n = (1:length(data))/length(data); % Calculate the n/N factor

with the line

n = (0:length(data)-1)/length(data); % Calculate the n/N factor

Note that in the standard DFT (and FFT), the indices for both $x[n]$, $X[k]$ range in $0 \leq n,k \leq N-1$.


Data representation of largest DCT coefficients

Your question is very accurate. Storing the largest  (1 % for instance) coefficients only from a sparsifying transform (DCT, wavelet, else) is fool's-gold, since you (more importantly the decoder) don't you where those are located.

And storing their binary indices counterbalances the compression gain for the sparsification. So trade-off should be found between theoretical sparsity and storage. For images, in JPEG, a low-frequency to high-frequency parsing scheme (zig-zag) is used. After quantization, lots of zeroed-out coefficients are compacted with End-Of-block symbols.  In JPEG 2000, a multi-scale and bit-wise  path is often  used, and zeroed-out coefficients are compacted with EZM-SPIHT like coding, that perform sort of vector coding of zero bit-plane.

So, some pre-determined way of parsing coefficients, based on prior knowledge on the signal or image structure, is often useful for actual compression, and applicable to specific class of data.


For a centre frequency of $F_C$ for a modulated signal $v(t)$.
$$
r(t) = v_I(t) \sqrt{2} \cos 2\pi F_C t - v_Q(t) \sqrt{2}\sin 2\pi F_C t
$$

Assuming a direct RF sampling at a rate of $F_S=1/T_S$ before downconversion, we plug $t=nT_S$ above.
$$
r(nT_S) = v_I(nT_S) \sqrt{2} \cos 2\pi F_C nT_S - v_Q(nT_S) \sqrt{2}\sin 2\pi F_C nT_S
$$

The expression $F_C nT_S$ above becomes $\frac{F_C}{F_S}n$ to yield
$$
r(nT_S) = v_I(nT_S) \sqrt{2} \cos 2\pi \frac{F_C
}{F_S} n - v_Q(nT_S) \sqrt{2}\sin 2\pi \frac{F_C
}{F_S} n
$$

Perhaps this is the ratio of the centre frequency to the sample rate you are asking about.

EDIT: Downconversion can be embedded in a downsampling operation as follows. In a conventional operation, we will have an LO downconverting the signal to baseband, filter it and then downsample it according to the symbol rate, as shown below.

Here, the filter is operating at a higher rate which is unnecessary when we have to throw $M-1$ out of every $M$ samples anyway. This is a hint for what we need to do.

First, we interchange the operations of downconvesion and filtering so our filter needs to be implemented in passband now. If the lowpass filter impulse response is $h[n]$, the passband filter response is
$$ g[n] = h[n] e^{j2\pi \frac{F_C}{F_S} n}$$

The interchange of downconversion and filtering is shown below.

Next, we slide the downsampling operation past the LO whose frequency gets multiplied by $M$, i.e., $F_C/F_S \cdot M$. We assume that
$$F_C = \frac{F_S}{M}$$
then our LO frequency after sliding it past the downsampler becomes
$$\frac{F_C}{F_S}M = \frac{F_S}{M}\frac{M}{F_S}=1$$
which implies that no downconversion is actually required!

However, remember that filter is still operating it at a higher rate. So we slide the downsampler past the filter now and implement it as a polyphase filter. The original spectrum is shown below.

This operation changes the spectrum such that all spectral replicas from $-F_S/2$ to $F_S/2$ end up at baseband overlapped on each other. The interesting part is that at the input of the polyphase arms, they possess a unique phase profile across the whole $360^{\circ}$ around the frequency IQ spectrum. A rough sketch is below.

This is only true in an exact manner for the sinusoidal signals, so the rest of the task for a signal with bandwidth is performed by the filter. When implemented as a polyphase filter, this filter operates with the delays matching those of the spectral replicas shown above. A final summation at the end cancels everything out leaving the desired spectrum.

This also implies that if the desired spectrum is at $F_S/M$ and we want to employ the initial lowpass filter $h[n]$ instead of the bandpass filter $g[n]$, a corresponding rotation is required after each polyphase arm to line everything up before the final cancelation.


Ok, its a trick question...its $Y(z) = X(z)$ because $y[n]$ and $x[n]_L$ are identical sequences.

$$
y[n] = \begin{cases}
x[n]_L & 0 \le n \le L-1\\
0&otherwise
\end{cases}
\longrightarrow  Y(z) = X(z)
$$


Why are Fourier analysis and transform only applicable for LTI systems?

The family of Fourier Transforms are specificaly developed for analysing frequency contents of the signals for which there is no definition of linearity or time invariance. Hence we can define the Fourier transform of any signal, as long as it's integrable (i.e. stable).

On the other hand we can also define a Frequency Response $H(\omega)$ for a particular class of systems (those being LTI) when they are characterised by a unique signal denoted by $h(t)$ which is called as the impulse response of the LTI system.

The Fourier transform $H(\omega)$ of the impulse response $h(t)$ is important because the complex exponentials $e^{j\omega t}$ are eigenfunctions of such LTI systems and when the input is decomposed into a sum of such complex exponentials, via a Fourier transform, so will be the output, multiplied with a complex gain given by the Frequency Response $H(\omega)$ applied for each particular frequency. This gives a considerable ease for the development of the mathematical theory of signals and systems.

On the other hand consider the following nonlinear system defined as:
$$y(t) = x(t)^2$$
For which you can apply CTFT to both sides to derive the result in frequency:
$$Y(\omega) = \frac{1}{2\pi} X(\omega) \star X(\omega)$$

Also consider the following linear but Time-Varying system defined by:
$$y(t) = g(t)x(t)$$
For which you can also express the relation in Frequency domain via CTFT as:
$$Y(\omega) = \frac{1}{2\pi} G(\omega) \star X(\omega)$$

Which I hope convinces you about the scope of Fourier Transforms and their applications to systems other than LTI.


Is it possible to do deconvolution with two data sets that have different sampling rates?

Resample your data to be close to the same rate. Note that the ratio between the time samples is 1.4990403....   if you resampled the 0.0781 data by a frequency ratio of 3/2, you would result with a time error of 0.000033 ps per cycle. From your sampling rate and block length you can determine how much of an impact this would have to see if that scale is sufficient: 0.0521/.000033 =  157.9 samples to slip $2 \pi$ radians. Resampling with higher ratios to support larger data block sizes: For instance 9369/6250 will get you to 1.49904 with a resulting time error of 3.07E-7 ps per cycle, with a block length of .0521/3.07E-7 = 169,650.6 to slip $2 \pi$ radians.

I have not tried what I am about to suggest next so my thought process could be flawed, but it seems that you could feasibly do the simpler 3/2 resampling ratio and compensate for the remaining frequency error with a quadratric phase versus frequency all pass function. This is because the residual frequency error is the error at the sampling rate, and the error at all frequencies will be in proportion of their frequency to the sampling rate (so a linear function over frequency). Phase is the integral of frequency versus time, but I am thinking of the possibility of a compensating quadratic phase (vs frequency) function to reduce the residual error. Note in this case the derivative of phase versus frequency is delay, so I believe that this would be similar to emulating a phase shift with a time delay; in that it will have the correct phase delay as anticipated at one frequency only, and then will have linear error (in phase) as you deviate from that frequency.  By getting our frequencies close with the 3/2 resampling, and assuming sufficient points in the fft such that the frequency spacing is small I could see the possibility that a smaller overall error is achieved with such an approach. This needs to be pursued further mathematically or experimentally but is a thought. If you do anything along this lines sucessfully please share your results! (Or if anyone sees the flaw with this approach that makes it not feasible please illuminate and I will delete this 2nd paragraph). I did post this technique as an additional question: Compensate for residual resampling error with a quadratic phase vs frequency- possible?


Designing high-pass filters using impulse invariance

You're right, you can't use the impulse invariance method for transforming continuous-time filters with a transfer function that is not strictly proper, i.e., for which the degree of the numerator is not strictly smaller than the degree of the denominator. High-pass and band-stop filters have transfer functions with numerator and denominator polynomials of the same degree, which means that the corresponding partial fraction expansion has a constant term. The latter corresponds to a (scaled) Dirac delta impulse in the impulse response, which obviously cannot be sampled.

Take a first-order high-pass filter as an example:

$$H(s)=\frac{s}{s+a}=1-\frac{a}{s+a}\Longleftrightarrow h(t)=\delta(t)-ae^{-at}u(t)\tag{1}$$

Now you can't determine $h(nT)$, $n\in\mathbb{Z}$, at $n=0$.

In sum, for the impulse invariance method to be applicable we require

$$\lim_{s\to\infty}H(s)=0\tag{2}$$


Understanding the resulting image matrix when differentiating image

Your interpretation is correct: directional derivation operators highlight variation in a given direction. Here, you use the $2$-point discrete derivative in the $x$-direction (along image rows). It may emphasize vertical features.

First, such operators indeed extend the initial image range. However, one often uses the absolute value of the derivative to compute a gradient. With the modulus, the derivative values fall back in the initial image range.

Second, as $[1 \;-1] = -[-1\; 1]$, and since convolution is a linear operator, you will get the opposite result as in the previous case. And it won't matter if you take the absolute value.


We are analyzing a real signal with the DTFT. Since we are using a limited number of samples it's like we are transforming a finite signal.

No. The DTFT takes an infinite discrete time signal as an input but the spectrum is continuous. I think you are confusing the DTFT and the DFT. They are different things.

  As I remember, the FT of a finite signal has an infinite band spectrum.

Yes

  So I suspect that the DTFT of a sampled signal (with a finite number of samples) always produce aliasing (since in the DTFT spectrum each replica is an infinite band). Is this correct?

Yes and no. Sampling a finite signal does with a finite sampling frequency will always induce some amount of aliasing, although in practice you can make it almost always "small enough" so it's acceptable for your application. That has nothing to do with DTFT (or DFT). Aliasing is caused by sampling, not during the analysis of the sampled data. IN other words, the aliasing happens regardless of whether you apply a DTFT or not.


The MATLAB / OCTAVE function freqz(b,a) is typically used for displaying the magnitude and the phase of a Frequency Response function $H(e^{j\omega})$ associated with a discrete-time, linear time-invariant (LTI), system with rational transfer function $H(z)$ (having a corresponding LCCDE representation) given by
$$H(z) = \frac{ \sum\limits_{k=0}^{M} b_k z^{-k} }{\sum\limits_{k=0}^{N} a_k z^{-k} } \tag{1}$$
where the numerator and denominator coefficients $b_k$ and $a_k$ represent such systems.
It's not meant to display spectrum of signals. Of course, an impulse response $h[n]$ associated with an LTI system can also be considered as a signal, and therefore its Frequency Response function $H(e^{j\omega})$ can also be considered as a Fourier transform of a signal ($h[n]$ in this case) too.
However, for such an interpretation to make sense, the signal $h[n]$ should be the impulse response of an LTI system. Stated in other words, to display the spectrum of $h[n] = \cos(\omega_0 n ) $, you should find the coefficients $a_k$ and $b_k$ that constructs a LCCDE whose impulse response solution is $h[n]$.
Consider a first order LCCDE (with initial rest):
$$ y[n] - A y[n-1] = x[n] \tag{2}$$
taking its Z-transform
$$ Y(z)(1 - A z^{-1}) = X(z) \tag{3}$$
and the System Function is
$$ H(z) = \frac{Y(z)}{X(z)} = \frac{1 }{1 - A z^{-1}}  \tag{4}$$
and the impulse response is (assuming causal)
$$ h[n] = A^n u[n] \tag{5} $$
If we let $A = e^{j\omega_0} $ , then the impulse response becomes
$$ h[n] = e^{j \omega_0 n} u[n] \tag{6} $$
From Eq-6, we have a single complex-exponential, by the choice of the coefficient $A =e^{j\omega_0}$. When combined with Euler identity of a causal cosine wave, we can see that
$$ h[n] = \cos(\omega_0 n) u[n]= 0.5( e^{j \omega_0 n} + e^{-j \omega_0 n} ) u[n] \tag{7} $$
Which can also be written with $A$ as
$$ h[n] = \cos(\omega_0 n) = 0.5( A^n + (A^*)^n ) u[n] \tag{8} $$
where $A^*$ is the complex-conjugate of $A$. Using the relation between Eqs.4&5, we can find the Z-transform of $h[n]$ as
$$ H(z) = \frac{0.5}{1- A z^{-1}} + \frac{0.5}{1 - A^* z^{-1} } \\ \tag{9} $$
equivalently
$$ H(z) = \frac{1 - 0.5(A+A^*) z^{-1}}{ 1- (A + A^*) z^{-1} + |A|^2 z^{-2} } \\ \tag{10} $$
we can construct back the differential equation that corresponds to this $H(z)$ as
$$ y[n] - 2\mathcal{Re}\{A\} y[n-1] + |A|^2 y[n-2] = x[n] -\mathcal{Re}\{A\} x[n-1] \\ \tag{11}$$
Plugging $A = e^{j\omega_0 n}$ into Eq-11 yields
$$ y[n] - 2 \cos(\omega_0) y[n-1] + y[n-2] = x[n]-\cos(\omega_0) x[n-1] \tag{12}$$
from Eq-12, you find the corresponding coefficients $a_k$ and $b_k$ as:
$$ b_k = [1, -\cos(\omega_0)] \tag{13}$$
and
$$ a_k = [1, -2 \cos(\omega_0), 1] \tag{14} $$
Now, you can display the Fourier transform of $h[n] = \cos(\omega_0 n)$ by using the freqz(b,a) function, with the coefficeints in Eqs-13,14.
Note that, $h[n]$ is only a causal sine wave valid for $n>0$. It's not the ideal two sided sine wave. Furthermore, I hope you are aware that the ideal sine wave's Fourier transform is a dirac impulse function. And what you can view freqz(b,a) or any other function is actually the spectrum of teh windowed cosine wave of a finite length.
The following MATLAB/ OCTAVE code shows you the implementation of above:
% spectrum of h[n] = cos(w0*n) to be displayed by freqz(b,a) function

N = 128;               % number of samples to evaluate the impulse response

w0 = 3*(2*pi/N);       % frequency of the cosine h[n] = cos(w0*n)
b = [1,-cos(w0)];      % coefficients b[k] of the LCCDE
a = [1,-2*cos(w0),1];  % coefficients a[k] of the LCCDE

h = impz(b,a,N);       % N samples of h[n] in n=[0:N-1] evaluated

figure,stem([0:N-1],h);
title('imp-resp of y[n] - 2 cos(\omega_0) y[n-1] + y[n-2] = (1-cos(\omega_0)) x[n]');

K = 1024;              % number of frequency points to be evaluated.
figure,freqz(b,a,K);
title('Freq-Response of the system: y[n] - 2 cos(\omega_0) y[n-1] + y[n-2] = (1-cos(\omega_0)) x[n]');

The outputs look like :

and

Note that I have chosen a perfect frequency for the cosine wave, so that it displayed almost like an impulse. It will not be the case in general.


Is there any trade-off between the sparseness and residual of the constructed data?

You can reformulate your problems as a minimization problem:

$$P1:  \min_x  \|Ax -b\|_2^2 + \lambda \|x\|_o $$

Where $\|x\|_o$ is the $\ell_o$ norm of x where it is the number of nonzero elements. Looking for a solution to p1 is essentially looking for a solution that best explains the model $\|Ax -b\|$ while having sparse coefficients. The problem is NP non convex NP hard and is extremely difficult to optimize. Another approach is to use the best convex envelope to the $\ell_o$ norm which is the $\ell_1$ norm.

$$P2:  \min_x  \|Ax -b\|_2^2 + \lambda \|x\|_1 $$

Under some conditions on A, the two problems are in fact equivalent in the sense they have the same global optimal solution. There are many efficient solvers for P2 including ADMMS, randomized methods, ISTA and FISTA.

To solve P2 in MATLAB, SPAMS is your friend here.


Anyone know what algorithm the Spice AC Noise Analysis uses?

Noise analysis in Spice (Berkeley Spice) is done by summing up the power spectral density from every noise source in the circuit.
There are a couple caveats.

The circuit is assumed to be linear.
In other words the circuit is first solved for a specific DC operating point then each of the components' equivalent linear noise sources are modeled as thermal noise sources.  Each of the components internal resistance values produce thermal noise whose frequency spectrum is shaped by the LINEAR response of the circuit.
Each source is assumed to be uncorrelated with any other.
That means that the power spectral density of the resulting noise is the sum of the power spectral density of each source individually. This is not necessarily the case when you have matched pairs of transistors that are closely coupled and share phonic coupling.
The 1/f noise is a complex parameter model.
1/f noise is difficult to model but there have been several standard models developed. If you stay away from low frequency in the analysis you can avoid dealing with this issue.

See google for more info on noise models such as Noise Sources in Bulk CMOS.


The RTL-SDR, as most software-defined radios, has an analog front-end that includes a tunable mixer. This mixer can down-convert the frequency band you're interested in to baseband, where it is sampled.

For example, let's say that you want to receive a signal centered around 900 MHz and with bandwidth $2B$ -- that is, the band covers the range from $900-B$ to $900+B$ MHz. You tune the mixer to 900 MHz and (after filtering) you end up with a signal centered around DC, covering the band from $-B$ to $B$ MHz. In general, this is a quadrature signal, that can be sampled at rate $2B$ complex samples per second.

As you see, with this setup, the sampling rate of the device only limits the bandwidth it can capture, without regard for the value of the carrier frequency (within the limits of the filters and mixers, of course).


How to extract IF (instantaneous frequency) from Discrete Hilbert Transform?

I have a nice solution for you that I'm always using whenever I need to compute the instantaneous frequency of a discrete analytical signal and which works much better then the common approach using $tan(.)$ where you need to introduce thresholds and fiddle around with adding/subtracting $2\pi$. In addition to that it is even faster.

Given a signal $x(t)$ and the corresponding analytical signal $$\xi(t) = x(t) + i\mathcal{H}x(t)$$ with $\mathcal{H}$ being the hilbert transform operator the instantaneous frequency $\omega(t)$ can be computed as followed. We have $\omega(t)=\dot{\varphi}(t)$ with $\varphi(t)=\Im log\xi(t)$ being the phase of the analytical signal $\xi$ ($\Im$ denotes the imaginary part). Combining these equations we get
$$
\omega(t)=\Im\frac{d}{dt}log\xi(t)=\Im\frac{\dot{\xi}(t)}{\xi(t)}
$$
This can also be seen be as follows:
$$
  \Im\frac{\dot{\xi}}{\xi} = \Im\frac{ \dot{A}e^{i\varphi} + iAe^{i\varphi}\dot{\varphi} }{Ae^{i\varphi}}
                       = \Im\Big[ \frac{\dot{A}}{A} + i\dot{\varphi} \Big]
                       = \dot{\varphi}
$$
Discretizing this equation you can now compute the discrete instantaneous frequency from the discrete analytical signal:
$$ \omega[n]\approx\Im\frac{\xi[n+1]-\xi[n]}{\xi[n]} $$
or second order:
$$ \omega[n]\approx\Im\frac{\xi[n+1]-\xi[n-1]}{2\xi[n]} $$

Of course this is not the exact solution as it is still requires a discrete differentiation of the signal. But this could be done using standard techniques.


$x$ is a vector, composed of the individual elements $x_1$, $x_2$, etc.  The "T" is for transpose.  It is not a variable.


You can implement very efficient high order low pass filter using cumulative summing. The theory is as follows. A perfect integrator has the recursive definition

$$y[n] = y[n-1]+x[n]$$

If you take the delayed difference of the output as in
$$y'[n] = y[n] - y[n-d]$$

you are realizing a boxcar filter of length $d-1$, as can be seen from explicitly resolving the recursion for $y[n]$ in $y'[n]$:

$$ y'[n] = (y[n-1]+x[n]) - y[n-d]$$

$$ y'[n] = ((y[n-2] + x[n-1]) + x[n]) - y[n-d]$$

$$ \dots$$

$$ y'[n] = y[n-d] + \sum_{k=1}^{d-1}x[n-k] - y[n-d] = \sum_{k=1}^{d-1}x[n-k]$$

That means you have a low pass filter of effective order $d-1$ from just one addition and one subtraction per sample.

Of course, boxcars are only useful up to a certain point. But the same idea can be generalized  to a larger class of filters with a trade off between filter shape and computational effort. And that is done by concatenating several integrator-difference stages.

Two in series will generate a triangular impulse response, three a parabolic one, and so on. Each stage convolves the output with yet another boxcar. In the frequency domain the original sinc() shaped response is taken to the n-th power (or if you have different delay lengths for the different stages, multiplied by other differently scaled sincs).

But both response sequences in time and frequency domain behave very predictably, because the central limit theorem of statistics predicts that the time domain response quickly approaches a Gaussian distribution. And because the fourier transform of a Gaussian is also Gaussian, you'll get a very good approximation of a Gaussian frequency response low pass filter after only a few stages.

One note about numerical precision though. Because you are integrating, the content of the delay buffer can be become very large and the difference precision suffers. This can be resolved by making the integrator leaky, so that it doesn't accumulate a constant function up to infinity. The leaky integrator difference equation is

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

for some sufficiently small positive $\epsilon$.


By using array antenna and a MUSIC algorithm, how many targets could be detected?

not a foolish question at all.

for a uniform linear array.

there are a number of different answers in the literature that do not agree.

the first issue is that direction finding and beamforming are different and have different upper bounds on the number of signals that can be resolved.

for beamforming $N-1$ sources under the narrow band assumption is typically the upper bound. you might assume polarization might double that number.

again under the narrow band assumption, the rank of the cross spectrum matrix is an upper bound for direction finding.
One might reasonably assume $N$ for $N$ sensors but one can construct a virtual array from the negative coarray.

The virtual array is based subtracting outputs of pairs of real sensors.  One can form $N^2$ such pairs if $x_k-x_j$ is distinct from $x_j-x_k$.  If one were to believe this, the upper bound at a single frequency is $N^2$ sources.

The $N^2$ result is given in the peer reviewed literature by an eminent IEEE Fellow but I can’t say it is compatible with my intuition.


The definition as you quoted it is only correct for static nonlinearities. The output of a static nonlinearity only depends on the input function directly, and not on its integral or derivatives. Otherwise, the nonlinearity is called dynamic (i.e. having memory). A static nonlinearity can also have memory if it is multi-valued, i.e. if there is more than one possible output value for a given input value (e.g. a system with hysteresis). Finally, a static nonlinearity is indeed memoryless if it is single-valued, i.e. if its output is uniquely defined by its input value:

$$y(t) = F(x(t))$$

where $y=F(x)$ is an injective function (i.e. each $x$ is mapped to only one value $y$).

So in sum you have

dynamic nonlinearity: has memory
static nonlinearity

multi-valued: has memory
single-valued: memoryless