Finding the fundamental frequency from Autocorrelated data

Three suggestions:

Try searching only over your expected pitch range: from a lag or period a bit shorter than that of your highest expected pitch for the instrument in question, to a lag a bit longer than twice the period of the lowest pitch expected.
Make sure there is a sufficient dip between your highest and your 2nd highest peak, or else you might find the shoulder of the 1st peak, and thus completely miss finding any 2nd peak.
You might also want to try using parabolic interpolation to better estimate each pitch period if any are between autocorrelation lags.


Frequency response of numerical derivative

When I try this, the results looks as expected. So you really have to explain in more detail what exactly it is that you're doing, because it doesn't seem to be a property of the discrete-time derivative.

EDIT: Now that I see your code, I'm convinced that the problem will disappear if you use

dydt = zeros(1,sr);

to initialize the derivative vector.


Is "Windowed Fourier Transform" a synonym for "STFT"?

In Chapter 2.4 Previous Work of  The Short Time Fourier Transform and Local Signals, S. Okamura, 2011, one reads:

  The STFT is also known under many names such as the windowed Fourier
  Transform, the Gabor transform, and the local Fourier transform.

Later, one discovers that the definition may  slightly vary with different authors, but, at first glance, they are the same. Be careful though. For some persons, a "Windowed Fourier Transform" is a mere Fourier transform of  a windowed signal. Because there is no sliding nor shift involved. For some, "short" implies "finite support".

Leaving theses versions aside, my personal interpretation and use  is that a “Windowed Fourier Transform” is a special case of a slightly more generic  “STFT” acceptation, for two (mild) reasons:

"windowed" somehow conveys the idea that only one window is used. But you can easily use different windows for STFT, long ones for low frequency, shorter ones for high frequencies.
the concept behind  "what a window is"  is not evident. However, for most people, a window  often has a non-zero sum (or unit), and often is either increasing/decreasing (unimodal, flat being a trivial version) or monotonic (like exponential window or weighting) for causal implementations. For short-term, it only needs to be "short", or well-concentrated in time, and could be more complicated.


How to Tell How Likely a Signal Is Present in Another One (Variance Unknown)?

Though the Matched Filter is the best tool detection of a known signals under AWGN it should work well here as well.

To say something about the probabilities the question is, do you know something about the energy of the received signals?
If you do, you should easily say something about the probabilities.

Pay attention that if the assumption is a signal could be within a different signal you need to apply a generalized Matched Filter (Multi Decision).


How to apply a mel-filterbank to a signal?

Binning is an averaging operation on the (squared) magnitudes of the DFT. You would maybe have $256$ DFT bins but only around $20$ outputs of the filter bank. So you need to average groups of DFT bins to reduce the dimension from $256$ to $20$. For a mel-scaled filter bank, the averaging functions (kernels) are usually triangular, i.e. the center DFT bins get more weight than the rest. As an example, let the first kernel that is used for averaging the lowest DFT frequencies be stored in an array $W_1[k]$, $k=0,\ldots,K-1$, then the output of that filter bank channel is computed as

$$y_1=\sum_{k=0}^{K-1}W_1[k]R[k]\tag{1}$$

where $R[k]$ are the (squared) DFT magnitudes. In this way you compute all $20$ filter bank outputs. You just need to know (define) the weights $W_i[k]$ and the range of DFT indices that are averaged by the respective kernel.

This answer to a related question might also be helpful (also have a look at the comments).


Determining Fourier Series coefficient for Discrete time

In the formula displayed for $x[n]$, the exponentials $e^{jk\omega_0 n}$ should be $e^{-jk\omega_0 n}$. I'm assuming that the sum $\sum\limits_{n = \langle N\rangle}$ is extended over all integers $n$ from $0$ to $N-1$. Fix $r\in \{0,\ldots, N-1\}$. Then

$$\sum_{n = \langle N\rangle}x[n]e^{jr\omega_0 n} = \sum_{n = \langle N\rangle}\sum_{k = \langle N\rangle}a_ke^{-jk\omega_0 n}e^{jr\omega_0 n} = \sum_{k = \langle N\rangle} a_k \sum_{n = \langle N\rangle} e^{j(r-k)\omega_0 n}.\tag{*}$$

If $k = r$, then each of the terms $e^{j(r-k)\omega_0 n}$ equals $1$, in which case

$$\sum_{n = \langle N\rangle} e^{j(r-k)\omega_0 n} = N.$$

If $k \neq r$, let $z = e^{j(r-k)\omega_0}$. Then $z\neq 1$ and $z^N = e^{2\pi j(r-k)n} = 1$. Thus

$$\sum_{n = \langle N\rangle} e^{j(r-k)\omega_0 n} = \sum_{n = \langle N\rangle} z^n = \frac{1 - z^N}{1 - z} = 0.$$

In summary,

$$\sum_{n = \langle N\rangle} e^{j(r-k)\omega_0 n} = N\delta_{k,r}.$$

Therefore,

$$\sum_{k = \langle N\rangle} a_k \sum_{n = \langle N\rangle} e^{j(r-k)\omega_0 n} = \sum_{k = \langle N \rangle} a_k (N\delta_{k,r}) = N\sum_{k = \langle N\rangle} a_k \delta_{k,r} = Na_r.\tag{**}$$

Combining (*) and (**),

$$Na_r = \sum_{n = \langle N\rangle} x[n]e^{jr\omega_0 n} = \sum_{n = \langle N\rangle}x[n]e^{jr(2\pi/N)n}.$$

Hence

$$a_r = \frac{1}{N}\sum_{n = \langle N\rangle} x[n]e^{jr(2\pi/N)n}.$$


Estimate amplitude of signal between 45 Hz-60 Hz

What you describe is a very simple signal model: one tone in a constrained frequency range, and superimposed non-sinusoidal noise.

You can do a lot of things based on that. For example, you could do the FM radio / lock-in amplifier dance, and just use a PLL (with a low control loop filter bandwidth, to remove most of the noise!) to lock onto the tone; using the resulting tone to mix your signal, aggressively low-pass filter it and simply take the average amplitude of the resulting signal.

You could also use (other) parametric estimators.

So, to me, estimation would roughly work like this:

band pass filter so that you only retain the signal of interest; preserving a bit of roll-off for the filter, that'd be something like filtering to 40–65 Hz. That gives you a remaining information-containing signal bandwidth of 25 Hz, so decimate your signal in the filtering step down to a sampling rate of 50 Hz. The rest of the signal doesn't contain anything you care about, so you're not losing anything by doing so!
The resulting signal will only contain the little noise power that falls into that bandwidth, plus the signal we care about. So, that's really a simple signal model now! Apply a suitable parametric estimator.

W.r.t. what a suitable parametric estimator would be:

If your SNR in this very reduced bandwidth is very good (and that's not unlikely!), calculate the analytic signal, and simply calculate the arctan of the ratio of successive samples of that. Apply some averaging/low-pass filtering to the result to get rid of remaining noise influences. It's not fast, it's not overly stable, but it's easy to grok (frequency is the angular velocity with which the signal runs along a origin circle in the imaginary plane). Frequency resolution: Defined by the accuracy of your arctan implementation (and of course limited by noise).
If your SNR is so-and-so, and you can live with limited resolution: Take $m$ samples of input. Concatenate $N-m$ zeros, so that you get an $N$-long zero-padded vector. Do the FFT of that. You'll get a PSD estimate with $N$ equidistant points within your bandwidth. Find the three strongest adjacent bins. You know the actual frequency of your tone is somewhere between these, on the maximum of a sinc running through these three bins. Since close to its maximum, a sinc is pretty similar to a quadratic function, use the highschool formulas to find the maximum of that sinc from the three points. Resolution is mainly limited by SNR, but this algorithm is somewhat more sensitive to noise than others.
Else, I'd recommend doing ESPRIT. A bit over the top, but gives you a frequency estimator that's pretty close to the Cramér-Rao Lower Boundary for short observations, and pretty stable in the presence of non-sinusoidal noise. Also, easy to adapt to more than one tone (simply change one parameter). The calculations involved (one vector-vector product, one quasi-pseudoinverse, two eigenvalue decompositions) are benign enough for your rates.


How to design a filter with a certain magnitude response

if you are only interested in designing a signal representing the wind turbine noise, you could just generate it from your magnitude. E.g., with your magnitude $A$ over frequency $f$ you can generate a random phase $\phi$ for each frequency and over a time $t$ and then add all together. Here a small code for illustration (I tried to use your varibles. However you need to check orientation etc.):

phi         = rand(size(fmag))*2*pi;
x           = (mag*sin(2*pi*fmag*t+repmat(phi,1,N)));

Hope this helps.


The result of an $N$-point DFT is an $N$ element complex vector.  For strictly real data, only the first $1+N/2$ bins matter, as the others are a redundant complex conjugate image.  For properly bandlimited data sampled at $F_s$, the frequency range is $0$ to $F_s/2$.  Divide $F_s/2$ by $N/2$, and you get $F_s/N$ as the frequency spacing of $1+N/2$ equally spaced DFT results, including both end points, DC and $F_s/2$.

So $F_s/N$ is the DFT result bin spacing of an $N$-point DFT.

However, the actual resolution can range from a fraction of the DFT bin spacing to $2$ or more DFT bins of separation, depending on the signal-to-noise ratio and what kind of resolution you want: frequency peak estimation, or peak separation (with a gap between peaks, typically of $3\textrm{ dB}$ or more).

In low noise (high signal-to-noise ratio), an isolated frequency peak can be resolved, by $\mathrm{Sinc}$ interpolation, to a much higher resolution than the DFT bin spacing, due to added information from the a-priori assumption (or a good guesstimate) of that signal-to-noise ratio and peak isolation.


Numerical Calculation of impulse response out of frequency response

I would suggest the frequency sampling design method, which results in a linear phase FIR filter. You sample the desired magnitude on an equidistant grid, add a linear phase term, and then you compute the impulse response by taking an inverse DFT. The resulting impulse response is usually windowed. I've explained the details in this answer, which also includes some Matlab/Octave code. Note that windowing is not done in that simple code fragment.


Figuring out how an RF pulse affects the MRI spins

Yes, it basically is that simple. If the resonance frequency of the spin falls in the frequency range that the pulse excites, the spin will be manipulated by the RF pulse. Reality usually kicks in when you have imperfect slice profiles (i.e. to achieve a perfect rectangular slice profile, you would need an infinitely long RF sinc-pulse).

Is it possible that you accidentally mixed up angular frequency and normal frequency? I am asking, because your frequency spread and RF bandwidth would excite $\frac{2000\,\textrm{Hz}}{639.7\,\frac{\textrm{Hz}}{\textrm{slice}}} \approx 3.12\,\textrm{slices}$, if I am not mistaken, and... when it comes to frequencies and a ratio is so close to $\pi$, and there is a "/2.0" in the original question, there could be something going on with wrong frequency units.


What does mean and variance do in Gaussian noise?

Does mean increases the brightness?

As you add a constant mean to each pixels of the image you will get an image that has an overall higher average image pixel value. Since images are generally encoded with the maximum pixel value being white, you would get an image which is closer to white, and hence perceptually appears brighter.

  Does variance increases the amount of noise?

That one deserve a little bit more background. Indeed, quoting from Wikipedia:

  In experimental sciences, noise can refer to any random fluctuations of data that hinders perception of an expected signal.

So in that context, even adding or subtracting a constant value to the entire image could be coined as noise (especially if this makes the image too bright or too dark so as to make it harder to discern objects in the image). If on the other hand we focus our attention to what is commonly described as a "noise" effect then this typically refer to some kind of perceived "dots" that occurs due to somewhat large variations in pixel values in close proximity. Since the variance is directly linked to the variation in pixel values, it is clear that an increase of the variance would increase the amount of noise (in the sense of perceived quantity of "dots").


It is a frequency dependent gain to compensate attenuation of the high frequencies in the signal band (< 0.6Hz). As a side effect it also increases the background noise (>0.6 Hz)


How does gaussian convution stacking work?

From this Wikipedia page, you can get:

So for your example, $\sigma_X = 1$ and $\sigma_Y = 2$. The result has a standard deviation of $\sigma_Z = \sqrt{\sigma^2_X + \sigma^2_Y} = \sqrt{1 + 4} = \sqrt{5}$.


DFT and perfect reconstruction of a square wave on a digital computer

First of all, it's a misunderstanding that the Gibbs phenomenon disappears if you use infinitely many Fourier series coefficients to reconstruct a discontinuous periodic function, such as a square wave. It doesn't. The reason is that generally the Fourier series doesn't converge point-wise, but it converges in the mean, i.e.,
$$\lim_{N\to\infty}\int_{0}^{T}\left|x(t)-\sum_{k=-N}^{N}c_ke^{j2\pi kt/T}\right|^2= 0\tag{1}$$
if $x(t)$ is a $T$-periodic function, and $c_k$ are its Fourier coefficients.
Taking the discrete Fourier transform (DFT) of a finite sequence of numbers just corresponds to a matrix multiplication:
$$\mathbf{y}=\mathbf{Ax}\tag{2}$$
and as long as the matrix $\mathbf{A}$ is invertible you can compute $\mathbf{x}$ from $\mathbf{y}$:
$$\mathbf{x}=\mathbf{A}^{-1}\mathbf{y}\tag{3}$$
This has nothing to do with the Fourier series of a continuous periodic function, and it has nothing to do with the Gibbs phenomenon.


How to change x-axis from samples to time?

You're getting an error because the two inputs to plt.plot need to be the same length. Notice that len(y)/48000 is just a number, so python tells us it of length 1.

Try plt.plot(range(len(e))/ 48000, e) **. What this does is change the values of your x axis to correspond to real time. If the sampling rate is $48000 \textrm{Hz}$, this means that $\textrm{sample}=1$ on your current axis now corresponds to 1/48000 of a second in the new axis, as this is the amount of time required to get one sample. Extend this to the whole axis - the time value on the new $x-$axis is simply the amount of time needed to take $y$ samples.

** if range() doesnt work, try numpy.arange() instead

Also, yes, these two lines

e -= e.min()
e /= e.max()

correspond to first shifting the signal so that its minimum value is now at 0, then dividing by the max. This is just a normalization, as you guessed.


What is the beneficient of taking the FFT of an ecg signal?

I'm not a doctor, I'm an engineer. However, I think I see the point of this.

An electrocardiograph is used as a diagnostic tool to assess the state of the heart and cardiovascular system. Several signals are measured, each of which can be used to make judgements about the health of the patient. Most of the signals measured in an ECG are periodic (in other words, they occur at a given frequency). The DFT/FFT is basically a way of analyzing the frequency content of a signal. So, it follows that by taking the DFT/FFT of an ECG, you can easily see those frequency components.

From the Wikipedia article on ECGs, here are some common components of an ECG and their associated frequencies:

P Wave: 0.08s, 12.5Hz

QRS Complex 0.09s, 11Hz

T Wave: 0.16s, 6.25Hz

I don't actually know the medical meaning of those features, but in any case, they're periodic, meaning you can analyze them using spectral analysis.


Sampling rate in Wifi 802.11 -> 20 MHz enough with a BW of 20 MHz?

The radio-frequency signal of 802.11n can have a bandwidth of 20 MHz or 40 MHz. But even if it is 20 MHz, it can be sampled at 20 MSa/s if an I/Q demodulator is used: at the receiver the radio-frequency signal is downconverted to the baseband. The I/Q demodulator has two lowpass output signals (Inhphase/Quadrature or real/imaginary part) with bandwidth 10 MHz each. Two samplers with 20 MSa/s each are then sufficient to fulfill the sampling criterion.

There are other receiver architectures where the received signal is first downconverted to an intermediate frequency (for example 10 MHz) and is then sampled at a higher rate (for example 40 MSa/s). In this case only one sampler is necessary.

The guard interval has only a negligible influence on the signal bandwidth. It is not included in the "Fourier transform period $T$", because it is only added after the transform. The symbol duration after the transmitter IFFT is therefore
$$
T = N/f_\mathrm s = 64/(20\,\text{MHz}) = 3.2\,\mathrm{\mu s}
$$


How to convert multi-band Energy Time Curve (ETC) to Impulse Response (IR)?

You can't really get the impulse response from the ETC. This is because ETC represents the energy of the signal, while the impulse response contains time information too.
One may argue that ETC also has some time information but it's not the same. You need phase to be able to reconstruct the impulse response and this is simply not present in the ETC.
What you could possibly try to do is to sum the energy of the frequency bands to acquire the total energy (broadband), calculate the square root of the result to acquire the amplitude and add random phase, or even impose the phase characteristics you want. Keep in mind that this may lead to unstable results or complex values in your time series (IR).
In order to impose your desired phase response, you could convolve your signal (square root of total energy) with the IFFT of a "manufactured" signal that will have constant amplitude spectrum and the desired phase.
Once more, please keep in mind that arbitrary phase curves may lead to results that "do not make sense". No matter what, the resulting sound may, or may not sound natural/real. Most probably (haven't tested though) parts of it, such as the late reverberation tale may sound natural (it exhibits random-like phase characteristics) while the fist parts of the IR may not sound convincing.
The main issue with the manufactured phase is that real spaces exhibit mixed-phase characteristics. Some parts of the IR (early reflections) have minimum phase (most often) while late reverberation shows maximum phase characteristics.


Choosing samples per symbol for modem pulse

The pulse rate is ultimately limited by the channel bandwidth. Assuming you're using raised cosine pulses with rolloff factor $\beta$, then $$B = \frac{(1+\beta)R_p}{2},$$ where $R_p$ is the pulse rate.

The number of samples per pulse is, in all libraries I know, independent of the pulse rate. Strictly speaking, you don't need more than $2B$ samples per second. However, some algorithms work better with more samples per pulse, and any signal plots you create will be clearer. In addition, normally you don't want your DAC to be working close to the theoretical limit.

In my experience (mainly with GNU Radio and USRP radios), 4 samples per pulse works very well with root raised-cosine pulses.

To further clarify, let us work through an example. Assume we want to transmit at a pulse rate $R_p=10,000$ pulses per second, with $\beta=0.5$. We know that we'll need a bandwidth of $7,500$ Hz. If we set our system to use 4 samples per pulse, then we'll have a sampling rate of $40,000$ samples per second. This is the amount of samples that our DSP system will need to process per second. If, instead, we go with 6 samples per pulse, then the sampling rate will increase to $60,000$ samples per second and the computational load of the system will increase accordingly.


Questions about the stability (and stationarity) of a system and state space representations

@w00zie: Hi: I'm putting this here just because I find it easier to write in here. This is not an adequate answer from a DSP point of view (and I would like to see one also) but let me just it explain more clearly since I babbled above and it's still not that clear.

Take the more standard AR(1):

$y_{t} = a y_{t-1} + \epsilon_{t}$.  ( denote this as first equation )
This can be re-written as

$y_{t}(1 - a L) = \epsilon_t$ where L is the backshift operator which is kind of the analogue to $z^{-1}$ in DSP notation. This expression can be re-written as

$y_{t} = \frac{\epsilon_{t}}{(1 - a L)}$ which can be re-written as

$y_{t} = \sum_{i=0}^{\infty} a^{i} \epsilon_{t-i}$.

But the only way that the infinite sum on the RHS can converge is if $|a| < 1$.
This is why that restriction arises for stability. I hope this helps even though it's still a time-series viewpoint rather than a DSP viewpoint.

Oh, and you asked why it looks so different from the standard AR(1) notationally.

You have: $s_t = 2  s_{t-2} +  3 w_{t-2}$.

If you let $y_t = s_t$ and $\epsilon_t = w_{t-2}$,  then that gives

$y_t = 2 * y_{t-2} +  3 \epsilon_{t}$. ( denote this as last equation ).

The equivalence of "first equation" and l"last equation" can be argued in two different ways:

1) you can think of last equation as an AR(2) with a zero $y_{t-1}$ coefficient. But the restriction on $a$ for stability doesn't change.

2) you can think of last equation as an AR(1) where the time step is doubled.

Note that the difference in the lags of the error terms, $\epsilon_{t}$ and $w_{t-2}$ is okay because $w_t$ is iid and uncorrelated with all other terms.
You are correct that, if $w_t$ was not iid and uncorrelated with everything else, then you couldn't use 1) or 2) as arguments for equivalence.

Note that I really shouldn't have used $y_{t}$ here since you use $y_{t}$ in a different context but my $y_t$ and your $y_{t}$ are not related.

Finally,  your $y_{t}$ is not stationary (i.e: not stable )  because it's just $s_t$ + noise and $s_t$ is not stationary.


Bose Chaudhuri Hocquenghem (BCH) vs Reed Solomon (RS) Error correcting codes and sound transmission channel modeling

What channel model fits data transmission over air best, I guess distortion generated by speakers should also be take under account?

Disclaimer: I'm really not an Audio person. There's a lot of audio DSP engineers on here that probably have a far better idea of what's happening on an acoustic channel. But:

Frequency selectiveness will be pretty important. It's unlikely your phone speakers can go below 40 Hz, or far above 12 kHz.
Audio typically might experience a bad case of multipath
prepare yourself for what would be a relatively solid amount of Doppler if your communication partners, or significant reflectors, or the medium, are moving

  Is there a general study comparing BCH vs RS in different channels?

probably, hundreds. But: it's really a thing that needs to be compared under specific channel models, and since there are freely available implementations of both, as soon as you have a realistic channel model for your specific use case down, simulation should give you more application-specific results than any study might ever give you.

I'd really go for channel sounding, in your case; send a known signal of suitable bandwidth, record the raw reception, try to extract a transfer function from that (assume LTI, and verify that assumption). Change the environment/channel, and do it again. Do that ten times (and write a paper), and come up with a stochastic channel model (e.g. "we have this general amplitude response for the ideal single channel, we have x~(PDF) multipath components, and their delays and weights follow the XYZ distribution", and write a paper).

Then write a channel simulator (at this point, I might, again, recommend GNU Radio, for which several channel sounders and channel simulators exist, but mainly for RF channels – but that's not inherent to its design, just to most of its userbase; obviously, write a paper), and test a few codes on that (and write a paper :))!


Harmonics below oscillation frequency in the brain

You explain in the comments that the signal is sampled at 50 KHz, then low pass filtered at 128 Hz with a 2nd order biquad filter and then resampled and stored at 256 Hz.

If you weren't prudent with choice of sampling rate and low pass filtering, then the harmonics of 10 Hz could create a component at 5 Hz through aliasing. For example this would occur if you used a sampling rate of 25 Hz: the third harmonic at 30 Hz would fold to the 5 Hz frequency. I explain aliasing in more detail at this post here:

Where should I set my anti-aliasing filter corner frequency for this signal?

That said, there are two prominent mechanisms that could create folding images in vicinity of 5 Hz due to aliasing during this sampling process.

The first would be due to interference components above the Nyquist frequency of 25 KHz prior to the sampling of the analog signal - specifically any frequency components at $N f_s \pm 5$ Hz would fold to 5 Hz in your digital signal where $N$ is any integer and $f_s$ here is 50 KHz. (So 50 KHz +/- 5 Hz, 100 KHz +/-5 Hz are example locations where interference can alias in). To avoid this be sure to use a good anti-aliasing filter and pay particular attention to its rejection around $N f_s \pm 5$ Hz with a rejection band around those frequencies that are within the frequency range that you want to maintain to be free of interference for your measurement.  Given this is many of orders larger than your 10 Hz interference, it is less likely this is an issue however should not be discounted depending on your measurement sensitivity desired.

The second mechanism, and likely dominant one in this case, is your resampling approach. Resampling to 256 Hz is a fractional rate conversion from 50 KHz and there could be many opportunities to create harmonic folding effects. The actual implementation needs to be reviewed in detail including significantly the performance of the resampling filter prior to decimation. Resampling to 256 Hz is identical to sampling at 256 Hz directly, so any frequency content above the Nyquist frequency of 128 Hz would fold into band. Specifically in this case the frequencies at $N 256 \pm 5$ Hz would land at 5 Hz in the digital domain.

Without seeing the actual implementation or ability to change it this would be difficult to predict although certainly can be measured. Induce a 10 Hz test signal with high harmonic content (impulses) and measure the results without your actual device you are measuring connected to assess your measurement noise floor in that condition.

Regardless of aliasing, to be sure that the interference itself does not contain 5 Hz components, I would also recommend analyzing a measurement without your true data source connected (such that you are only evaluating background noise). This would also give you a fidelity of your measurement system in showing you what sensitivity can be achieved, and the related confidence of your results.


Is there any method to find K value in KNN?

Basically too large k may tend to underfit, while smaller k may tend to overfit. You may try to run the classification with different k within the preset range with cross-validation. Observe and find out the k that gives you the least cross validation sum-of-square error.

Another possible metric besides the sum-of-square error is called Silhouette applied to interpretation and validation of clusters of data. The technique provides a succinct graphical representation of how well each object lies within its cluster.


DFT for signals of different length and sampling rate

Fft normalization is often an overlooked topic and many articles dealing with FFT just say to divide the fft ouput by N (divide by either N or N/2 depending on the actual
FFT algorithm used).
Real FFT  usually discards the upper part of the spectrum (due to FFT redundancy) so the normalization coefficient here is N/2, rather than N.

  So, my question is: what normalization coefficient for PSD should I
  use in my case?

I prefer the approach taken by Heinzel et al, and according to them, the normalized/scaled PSD  can be computed as

PSD = ( 2 * (yi * yi) ) / ( fs * S2 ); // i = 0…N/2  (1)

S2 = sum(wi * wi); //i = 0...N (2)

(yi is the complex result of the fft at bin i and wi is the value of the window function at point i )

The sum S2 is used instead of N in order to better estimate the scaling
factor, given a particular window (rectangular, Hann etc). Using no windows at all corresponds to a rectangular window and in the case of a rectangular window S2 is equal to N.

Heinzel et al also define a metric called ENBW (equivalent noise bandwidth), which is at the heart of the here defined scaling technique and for those interested, complete derivations (including 1 and 2) can be found here https://holometer.fnal.gov/GH_FFT.pdf

Note also that for a rectangular window, ENBW is equal to fsampling / N, which may explain why the results started to make sense for you once this normalization factor was applied.


The magnitude is given by $ {\left| X \left( f \right) \right|} $ hence it is 6.
Assuming all data is given in the graph, the function isn't complex it is only negative -> Phase is -180 degrees.
You should use its symmetry. Without over thinking it, It could be generated by a window on the center convolved with 2 deltas (Shifted) -> Sinc in time multiplied by 2 exponentials.


Generally Bayer mask (also known as CFA = color filter array) imagers are the most popular.

3Chip cameras cost considerably more - you are not only paying for 3x as many chips but the alignment in manufacturing is expensive and the range of cameras is much more limited. The largest are HD (1920x1080) although a lot of cameras actually have much smaller chips and rely on a deliberate half-pixel misalignment to interpolate each colour into a higher resolution - as with a CFA. Since most are aimed at the TV market they are usually video output (so you need a grabber card) and are a fixed TV frame rate

There are potentially aliasing issues with CFA cameras if you have a very small point image which may hit only one of the colour cells, this can be fixed with an anti-alliasing filter which blurs the image slightly, or by software.

One area in which 3chip wins is in low light where you need high colour quality. A 3chip camera detects essentially all the light in the scene in each colour. A CFA will detect only 25% of the Red and Blue and 50% of the green since the other photons are blocked by filters.


Your intuition is correct. The reason is that the rectangular pulse train has  discontinuities, whereas the rectified sine is continuous (but its derivative has discontinuities). The consequence of this is that the Fourier coefficients of the rectangular pulse train decay as $1/n$, whereas the coefficients of the rectified sine decay as $1/n^2$ (where $n$ is the index of the Fourier coefficients).


Bounded or unbounded in that context always refers to amplitude. E.g., BIBO (bounded-input-bounded-output) stability means that a bounded input signal $x(t)$ satisfying $|x(t)|\le A$ always results in a bounded output $y(t)$ satisfying $|y(t)|\le B$ with some finite positive constants $A$ and $B$.

So the Dirac delta impulse is not well-suited for testing stability of a system. E.g., an ideal integrator's response to an impulse at the input is a unit step, which is of course bounded. Yet, the ideal integrator is unstable because its response to a step is unbounded.

In sum, a Dirac delta impulse is unbounded.


Transfer function of linear phase filter & its frequency domain representation

You finally have changed your $H(z)$ and added a qualifier "linear phase" , so is the new answer:

Now assuming that you have already determined $H(z)$ from the problem specifications, so the rest of obtaining DFT $H[k]$ can be approached in a number of ways one of which (perhaps longer than necessary) is the following:

Using the linearity and shift properties of the Z-Transform, we can show that that if a given signal $x[n]$ is of the form $$x[n] = a\delta[n] + b\delta[n-1] + c\delta[n-k]$$ then its Z-Transform is of the form: $$X(z) = a + b z^{-1} + c z^{-k}$$

Applying this observation, in the reverse direction, to $H(z)$ in your problem we get a prototype for $h[n]$ as: $$h[n]=k(\delta[n] - {5\over6} \delta[n-1] - {19\over3} \delta[n-2] - {5\over6} \delta[n-3] + \delta[n-4])$$

Also as you state that $h[0]=2$, yields $k=2$, hence the signal $h[n]$ is: $$h[n]= 2\delta[n] - {5\over3} \delta[n-1] - {38\over3} \delta[n-2] - {5\over3} \delta[n-1] + 2\delta[n-4]$$

Obtaining DFT $H[k]$ of this $h[n]$ is actually a very similar procedure to obtaining $H(z)$, you shall simply plug $e^{j{2\pi \over N} k}$ in place of $z$ where you simply get this:

$$H[k]= 2 - {5\over3} e^{-j{2\pi \over N} k} - {38\over3} e^{-j{4\pi \over N} k} -{5\over3} e^{-j{6\pi \over N} k} + 2e^{-j{8\pi \over N} k}$$

where $k$ ranges from $0$ to $N-1$ and plugging $N=5$ produces:
$$H[k]= 2 - {5\over3} e^{-j{2\pi \over 5} k} - {38\over3} e^{-j{4\pi \over 5} k} -{5\over3} e^{-j{6\pi \over 5} k} + 2e^{-j{8\pi \over 5} k}$$

Note that eventhough it is perfectly valid to state that, under suitable conditions, $H[k]=H(z)$ evaluated at $z=e^{j{2\pi\over N} k}$ , I guess it does not make an answer...


There are several different ways to look at noise:

The probability density function (PDF) of a noise process, $x(t)$ say. It is sometimes written $p_x(x)$.  This tells you how frequently the the variable $x(t)$ takes on specific values.

One way to think of a PDF is as a (normalized) histogram of the possible values that $x(t)$ can take.  The PDF is normalized because $x(t)$ must take on some value between $-\infty$ and $+\infty$, so the integral
$$
\int_{-\infty}^{+\infty} p_x(x) dx = 1.
$$

For example, the figure below shows the PDF (in blue) and the histogram (in black) for a normal (a.k.a. Gaussian) random variable and for a uniform random variable.

Another way to look at noise is whether there is a relationship between the value of $x(t_1)$ and the value $x(t_2)$ (where $t_1 \not= t_2$). If there is no relationship, then the noise is said to be "white".  If there is a relationship, the noise is said to be "colored".

The way this is usually illustrated is using autocorrelation (autocovariance) of the signal.

Below are two plots of the autocorrelation of normally distributed noise.  In the first plot, the noise is white.  In the second plot, there is a high level of correlation between subsequent samples.

scilab CODE BELOE

// 15426

X_gauss = rand(1,1000,'normal');

X_uniform = rand(1,1000,'uniform');

figure(1)
clf
subplot(211)
histplot(100,X_gauss)
vals = [-3:0.1:3];
plot(vals,exp(-vals.^2/2)/sqrt(2*%pi))
subplot(212)
histplot(100,X_uniform)
plot([0 1], [1 1])

figure(2)

X_white = X_gauss;
X_colored = filter(ones(1,50),10,X_white);

figure(2)
clf
subplot(211)
plot(xcorr(X_white))
mtlb_axis([970 1030 -1000 1000])
subplot(212)
plot(xcorr(X_colored))
mtlb_axis([970 1030 -1000 1000])


Getting audio samples from power spectrum

No, the power spectrum only holds half of the information that the DFT'ed signal held, which held all the information that the original audio signal held.

Sadly, that half is important, and not reconstruction is possible. You might come up with an estimator that "guesses" the phase information you're losing by doing the magnitude square, but it will not recover the original signal in the general case. I'm not aware of any such estimators that do something useful.

In short: you can reverse the DFT, it's a reversible operation. You cannot revert calculating the magnitude square of complex numbers – that is a irreversible operation, and hence, you can't recover the original signal from a PSD estimate.

To illustrate this:

Imagine you record one hour (let's say from $t=-18000$ to $t=1800) of perfect silence, and in the middle of that is one second of a pure 100 Hz sine tone.

The DFT of that is pretty simple: In time domain, it's a sine (which is infinite) multiplied with a unit rectangular window. We know the spectrum of that; it's the convolution of the Fourier transform of the sine (which is dirac with "area" $\frac 12j$ at the positive frequency, and $-\frac12j$ at the negative frequency) convolved with the Fourier transform of a rectangular window (which is always some scaled $\text{sinc}$ function). Hence, we know that the DFT of that recording must be the sum of two shifted sincs, one of those being "mirrored" on the frequency axis.

Now, we estimate the PSD by just calculating the magnitude square point-wise, and we end up with something that looks like two positive lumps at the positive and negative frequencies. For obvious reasons, it's symmetric to the magnitude axis. We call the DFT $X(f)$, and that leaves us with the PSD estimate $|X(F)|^2$

We now decide that the 1s of tone shouldn't be at the $t=0$, but let's say at $t=900$; now, instead of recalculating the whole spectrum, we just apply the time-shift property of Fourier transform: a $\Delta t$ shift in time domain is a multiplication with $e^{j2\pi f \Delta t}$ in frequency domain. That way, we just need to multiply with this function in frequency domain.

And lo! Let $X(f)\cdot e^{j2\pi f \Delta t}$ be our spectrum, and $|X(f)\cdot e^{j2\pi f \Delta t}|^2$ be our PSD estimate. Observe that the modulation with the complex $e$xponential doesn't change magnitude square – hence, you can't tell time shifts from your PSD alone.


Question regarding using cardinality in pattern recognition

There's a lot of red flags in that paper that shouldn't have passed review. In fact, I'm a bit surprised – Mr. Ortiz-Catalan is a well-published researcher in a research group called "signals and systems" at a very prestigious university, and yet there's sentences like:

  Altering the dimension of the unit used for sampling (ADC
  resolution) to a high precision unit (for example, a double) would
  alienate the discrimination power of cardinality.

There's no ADC that give double precision floating point values. (Aside from "double precision" being a relative term, only defined within one programming language or CPU architecture). At any rate, the "D" in ADC stands for "digital", and that means "discrete in time and value, so

  This is because
  every sample value would be unique, so cardinality will always be
  equal to the number of samples in the time window (maximum
  value of cardinality).

is class-A hogwash, because sampling leads to discrete samples, and thus, the maximum number of unique values you can get is the cardinality of the alphabet of your ADC.

The same applies to:

  A simple operation like dividing the EMG sampled values
  by the gain of the amplifiers, in an attempt to use the original
  EMG amplitude in further computations, could increase the unit
  precision automatically by the processing software; that is, from
  the original acquisition values (for example, 14 bits) to a double,
  in the case of MATLAB (Massachusetts, United States).

For a fixed gain, you get as many unique values out of that division as you put in – division by a constant is a function (school math), for the same input value it should give the same output value; the data type doesn't matter here (it might, in a corner case, mean that the number of output values is smaller than the number of input values. Never more. And generally: Information can't be generated during signal processing – they cite Information Theory as basis for their research, so they must know the Data Processing Theorem, and gaining values would be equivalent to generating information).

I have to credit the author for the "I cite a software tool as coming from Massachusetts" joke – I liked that, it was cute :)

Having cleared up that ADCs only give discrete values, let's move forward.

Every ADC I encountered in EMG applications so far is used with a proper anti-aliasing filter up front. Needing to do that is mathematical/physical fact – if you don't band-limit the analog signal before sampling it, you get artifacts, namely aliasing in your signal, and also, especially for "peaky" signals such as EMG tends to be, instead of noticing every peak, you'd, in the best case, just get a non-1 probability of observation.

So here we go, feeding in band-limited signals into an ADC with a finite set of values.

Now, I mainly work in communication technology, so I'm much more familiar with the mathematical tools meant for analysis of harmonic signals, but we can still work with this pretty well:

Nyquist's sampling theorem says that, if the original signal is sufficiently band-limited, we can reconstruct the original signal from a periodical observation of instantaneous values. That periodicality is the sampling rate, and in the real-valued case, that needs to be just a bit higher than twice the highest frequency in the signal (i.e. the signal's bandwidth).

Nyquist doesn't require any synchronization between the signal and the sampling.

Thus, let's do the following thought experiment:

Sample a signal with bandwidth $B$, at times $t_i = i T_s$, $T_s < \frac1{2B}$.
Following Nyquist, we now have a complete representation of the signal
We find that if we configured our sampling system optimally, the extreme values of the sampled signals are at least close to the extreme values that the ADC can produce – otherwise, we'd be wasting dynamic range and could very well better represent the signal by amplifying it a bit before, reducing the quantization errors. This should have tipped us of that cardinality is extremely sensitive to the physical voltage amplitude – and thus a very hard-to-control meausure, considering that depends on a multitude of factors like electrode placement, contact quality, temperature of the system...
Now, the information how many different values where taken by the signal is identical to how many were NOT taken – since we know how many values our ADC can produce.
The observational windows were 200ms @ 2000Hz sampling rate, i.e. 400 samples. Hence, there were exactly 400 different cardinalities that could be realized here: 1 unique value (constant) or up to 400 different values. We can, with a high degree of certainty, say that certain cardinalities are more likely than others – for example, the constant signal is very unlikely to happen, and 400 different values isn't much more likely, either, since the signal is discrete and supposed to be correlated.
The paper finds that decreasing the sampling rate to $\frac12$ of the orignal 2 kHz only reduces the quality by 0.2 % – that is surprising, as half of the signal goes missing! We must pick one of three hypotheses:

A) the signal digitized was sufficiently bandlimited in the first place for the lower sampling rate, and in that case oversampling was done, which should reduce the meaningfulness of value set cardinality (as oversampled signals will just take more intermediate signals, thus shifting the density functions for that measure closer together higher up in of the range of possible cardinalities). The interesting part is here that additive noise will effectively lead to dithering – ie. an increase of cardinality, regardless of the original signal.
B) the signal was undersampled when it was sampled with 1 kHz, by a factor of $u = \frac{2B}{f_\text{sample}}$.
C) Most of the energy in the signal was in the lower half (or, in fact, the lower fourth, as 500 Hz sampling rate still seems to work nearly as good – which at some point should make one wonder what we're really looking at), but not quite all – and thus, the fast changes that (almost) always lead to unique values are rare and correlate highly with the classification.

Intuitively, C) works quite nicely: anything that changed by more than 1 quantization step between two successive samples at the 1 kHz rate will almost always have one intermediate value (excluding strong noise for now) sampled at 2 kHz. But: Knowing typical EMG spectra, there's less enery at higher frequencies (i.e. "faster jumps"), and thus, it's very unlikely that strong peaks and intermediate values from these will be duplicate within an observation window (especially if that is only 400 samples long).
Thus, the counting of unique values is but an aliasing detector – or, if you want so, a detector for energy in the upper half of the spectrum. A high-pass filter!

I personally have my problem with papers that offer $p$-values but don't even mention how they calculated those. Maybe you could nicely ask to get their algorithms/software. This should be no problem – the paper boldly boasts "open access" insignia, and that ideologically implies that they're proud to share not only their results but also their methodology. I do believe that the paper is on to something – it seems to have demonstrated a classifier that outperforms others. However, these others might have been badly parameterized, and might not perform that much worse, especially on properly preprocessed data.


How to implement a $1$-bit quantizer in Simulink?

You can see this example of Sigma-Delta A/D Conversion.


First a comment on your noise-generation process. The Matlab function randn() generates Gaussian noise with zero mean and mean power 1. So steps 2 and 3 are obsolete.

If you'd like to achieve a given SNR, then creating the noise signal with the wanted power is as simple as

w = w .* sqrt(P_n)

Where P_n is the mean noise power and can be calculated by the equation you've already found yourself. Actually you're doing the exact same thing in your step 4, from which you can also derive that

P_n = W0 * fs

assuming that fs is the sampling frequency of your system and thus equal to the bandwidth and using Pw = 1.


Complex Cepstral Transform on Audio in MATLAB

Natural log in matlab is 'log(x)' and log to the base 10 is 'log10(x)'

More information on Mathworks


obtaning discrete-time system impulse response $h[n]$ from step response $s[n]$

Considering an LTI system, its impulse response $h[n]$ and the step response $s[n]$ are related by $$h[n] = s[n]-s[n-1]$$ and $$s[n] = u[n] \star h[n] = \sum_{k=-\infty}^{n} h[k] = \sum_{k=0}^{\infty} h[n-k] $$

Given your definition of $s[n]$ therefore its $h[n]$ is
$$h[n] = \left[5 - 4 \left( \frac{4}{5} \right)^n \right] u[n] -\left[5 - 4 \left( \frac{4}{5} \right)^{n-1} \right] u[n-1]$$

After simplifications it would become:
$$h[n] = 5\delta[n] + 5 \left( \frac{4}{5} \right)^{n} u[n-1] - 4 \left( \frac{4}{5} \right)^{n} u[n]$$ further simplifies to

$$ \boxed{ h[n] = 5 \left( \frac{4}{5} \right)^{n} u[n] - 4 \left( \frac{4}{5} \right)^{n} u[n] = \left( \frac{4}{5} \right)^{n} u[n] }$$

Looking at the signal flow graph and using Z-transform variables you may achieve the same result by considering:

$$Y(z) = \frac{4}{5} Y(z) z^{-1} + X(z)$$
$$Y(z)( 1 - \frac{4}{5} z^{-1}) = X(z)$$
$$ H(z) = \frac{Y(z)}{X(z)} = \frac{1}{1 - \frac{4}{5} z^{-1}} $$

From which the inverse Z-transform (with the causal ROC) yields:
$$ \boxed{ h[n] = \left( \frac{4}{5} \right)^n u[n] }$$ which is the same result obtained from the time-domain approach.


Take the chirp signal

$$
x(t) = e^{j\pi\alpha t^2} = e^{j\phi(t)}
$$

Where $\alpha$ is the chirp rate of the signal. You can see that the function describing the phase, $\phi(t)$, is in the form of a quadratic equation. So you can plot

$$
\phi(t) = \pi\alpha t^2
$$

and see the quadratic nature of the phase at any point in time. The actual phase-change between each sample depends on what your sample rate is.


Find the length of the impulse response for the given output and input

You correctly figured out that the system must have $3$ zeros. So you just have to define a corresponding transfer function $H(z)$ that satisfies
$$H(z_i)=0,\qquad z_i=e^{j\omega_i},\quad i\in\{0,1,2\}$$
In other words, what is the minimum number of coefficients of a polynomial with $3$ zeros? Note that you might need to provide two solutions: one general solution, and one with the restriction that the polynomial coefficients be real-valued (as is required for a real-valued system).


What is the maximum allowed transition band of the low-pass filter used in the reconstruction kernel for a CD with maximum speech frequency of 4 kHz?

If you don't want DAC noise near 40-44 kHz to alias down into your audio, you want the narrower/lower transition band reconstruction filter.


Sudden loudness changes producing unwanted noises?

You are changing the volume in discrete steps, applied at random points in the audio signal.

If the signal happens to be close to zero, you won't hear much of anything.

If the signal is close to a peak, then you will get a "step" change in the signal.  This audible.

When you roll the volume up and down that way, you get a kind of "zipper" sounding noise.

What you want to do is to apply volume changes only when the signal is close to zero.  The "zippers" still occur, but are much less annoying - and may be completely inaudible if you change volume only at times when the signal is at zero.

Looking at the code you referred to in a comment, it appears you aren't changing volume when the amplitude is close to zero.

You are changing volume when the signal is less than zero.  What you want to do is to change volume when the absolute value of the signal is less than some small value.

Start with like 0.1, then try 0.01.  You should notice far less "zipper noise" that way.

Sort of like this:

if (abs(amplitude) < 0.1d) {
   currentVolume = targetVolume;
}


Machine learning for denoising MRI images

My personal feeling is that you should do each things separately and compare the results.  For example, take your MRI dataset and denoise using "standard algorithm 1", "standard algorithm 2" and "neural network algorithm 1".  I would keep things reasonably simple unless you have a good justification of doing "standard algorithm 2 and neural network algorithm 1".   Doing two in a row is going to be harder to justify in a paper.

(I think the other question that is difficult is how to tell if one denoising algorithm is "better" than another.  I have written a couple papers on standard image processing denoising in MRI and it is always a difficult question to answer.  Also, if the output of the denoising is a dataset that is to be interpreted by a radiologist then other consideration would have to happen vs denoising as a pre-step to some scientific quantification.)

One thing you could do is look through papers, for example at https://arxiv.org/ to see if there are other neural network denoising examples that can help guide you in parameter selection.


You can use the Laplace transform, but can also simply use convolution in the time domain. In any case, you will need the system's impulse response $h(t)$. Let the input signal $x(t)$ satisfy $x(t)=0$ for $t<0$, and $x(t+T)=x(t)$ for $t>0$ and $T>0$, as required. Furthermore, let $f(t)$ be the first period of $x(t)$:

$$f(t)=\begin{cases}x(t),&0<t<T\\0,&\text{otherwise}\end{cases}\tag{1}$$

Then

$$x(t)=\sum_{n=0}^{\infty}f(t-nT)=f(t)\star\sum_{n=0}^{\infty}\delta(t-nT)\tag{2}$$

where $\star$ denotes convolution, and $\delta(t)$ is the Dirac delta impulse. If $g(t)$ denotes the convolution $h(t)\star f(t)$, i.e., the system's response to $f(t)$, then the output signal can be written as

$$\begin{align}y(t)&=h(t)\star x(t)\\&=h(t)\star f(t)\star\sum_{n=0}^{\infty}\delta(t-nT)\\&=g(t)\star \sum_{n=0}^{\infty}\delta(t-nT)\\&=\sum_{n=0}^{\infty}g(t-nT)\tag{3}\end{align} $$

According to $(3)$, the output signal can be written as a sum of shifted responses to the finite length signal $f(t)$, which corresponds to the first period of the input signal.

The Laplace transform of the output signal can also be written in terms of the transforms of the input signal and of the function $f(t)$. Note that

$$f(t)=x(t)-x(t-T)\tag{4}$$

The Laplace transform of $(4)$ is

$$F(s)=X(s)(1-e^{-st})\tag{5}$$

If $H(s)$ is the system's transfer function, i.e. the Laplace transform of the impulse response $h(t)$, then the output is given by

$$Y(s)=H(s)X(s)=H(s)\frac{F(s)}{1-e^{-st}}\tag{6}$$


How to prove $cos(t) + cos(\pi t)$ is non periodic function? Also can I represent this signal using fourier series?

$cos(t)$ and $cos(\pi t)$ are both periodic with periods $2\pi$ and $2$ respectively. To find the period of the sum we need to find an integers $n,m$ such that $\pi n = 2 m$ or $\pi / 2 = m /n $, which is not possible since $\pi$ is irrational.

I think the Fourier series uses specific frequencies to form the signal. You can't necessarily put a frequency component at irrational frequencies.

Fourier series work fine with irrational frequencies too, but the frequencies are all harmonics, that is, integer multiples, of the fundamental frequency.


I think this code would be helpful, as displaying intensity part gives you the spatial characteristics of an image.

%Write a MATLAB function which converts an RGB to HSI. Display the Hue
%image, Saturation image and the Intensity image. Histogram equalize the H,
%S and I images and the HSI image. Then convert the HSI image back to RGB.

function rgbtohsi(x)
F=imread('C:\Documents and Settings\admin\Desktop\TM\tanya\img3.tif ');
F=im2double(F);
r=F(:,:,1);
g=F(:,:,2);
b=F(:,:,3);
th=acos((0.5*((r-g)+(r-b)))./((sqrt((r-g).^2+(r-b).*(g-b)))+eps));
H=th;
H(b>g)=2*pi-H(b>g);
H=H/(2*pi);
S=1-3.*(min(min(r,g),b))./(r+g+b+eps);
I=(r+g+b)/3;
hsi=cat(3,H,S,I);
HE=H*2*pi;
HE=histeq(HE);
HE=HE/(2*pi);
SE=histeq(S);
IE=histeq(I);
choice=input('1:RGB to HSI\n2:Display Hue, Saturation and Intensity Images\n3:HSI to
    RGB\n4:Hue-Equalization\n5:Saturation-Equalization\n6:Intensity-
    Equalization\n7:HSI-   Equalization\n Enter your choice :');
switch choice
case 1
    figure,imshow(F),title('RGB Image');
    figure, imshow(hsi),title('HSI Image');
case 2
    figure,imshow(F),title('RGB Image');
    figure, imshow(H),title('Hue Image');
    figure, imshow(S),title('Saturation Image');
    figure, imshow(I),title('Intensity Image');
case 3
    C=hsitorgb(hsi);
    figure,imshow(hsi),title('HSI Image');
    figure, imshow(C),title('RGB Image');
case 4
    RV=cat(3,HE,S,I);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-Hue Equalized');
case 5
    RV=cat(3,H,SE,I);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-Saturation Equalized');
case 6
    RV=cat(3,H,S,IE);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-Intensity Equalized');
case 7
    RV=cat(3,HE,SE,IE);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-HSI Equalized');
    F1=C;
    otherwise
    display('Wrong choice');
end
end
    function C=hsitorgb(hsi)
    HV=hsi(:,:,1)*2*pi;
    SV=hsi(:,:,2);
    IV=hsi(:,:,3);
    R=zeros(size(HV));
    G=zeros(size(HV));
    B=zeros(size(HV));
%RG Sector
    id=find((0<=HV)& (HV<2*pi/3));
    B(id)=IV(id).*(1-SV(id));
    R(id)=IV(id).*(1+SV(id).*cos(HV(id))./cos(pi/3-HV(id)));
    G(id)=3*IV(id)-(R(id)+B(id));

%BG Sector
    id=find((2*pi/3<=HV)& (HV<4*pi/3));
    R(id)=IV(id).*(1-SV(id));
    G(id)=IV(id).*(1+SV(id).*cos(HV(id)-2*pi/3)./cos(pi-HV(id)));
    B(id)=3*IV(id)-(R(id)+G(id));
%BR Sector
    id=find((4*pi/3<=HV)& (HV<2*pi));
    G(id)=IV(id).*(1-SV(id));
    B(id)=IV(id).*(1+SV(id).*cos(HV(id)-4*pi/3)./cos(5*pi/3-HV(id)));
    R(id)=3*IV(id)-(G(id)+B(id));
    C=cat(3,R,G,B);
    C=max(min(C,1),0);
end


Plotting magnitude and phase for frequency spectrum

for plotting magnitude spectrum versus frequency you can use the following script:

t = 0:0.001:0.2;
s = sin(2*pi*100*t)+cos(pi*500*t);
S = abs(fftshift(fft(s)));

Fs = 1000;
df = 1000/(length(s));
f = -Fs/2:df:Fs/2-df;

plot(f,S);

Here I used your first signal and plotted its magnitude spectrum. You must notice for your signal min choosable Nyquist frequency is 500 Hz But we use it two times of it(which is 1000Hz), because we also want to see the part before f =0 in the frequency spectrum. Using this we find df(delta f) frequency spacing. frequency vector must be multiples of df. By this way spectrum will be drawn correctly.

Here is the output:


What's the meaning of a complex zero/pole?

We usually talk of $j\omega$ when we're also interested in the Laplace transform of a signal / system, but want to just talk about the frequency response.

The physical meaning of the imaginary part is that it refers to purely sinusoidal signals and are constant "amplitude". The real part refers to signals for which the "amplitude" decays or grows exponentially.


Approximating Lorentzian Fourier transform with FFT

First of all, take a look at this very related question and its answer(s). Second, your analytic solution is wrong, so it's no surprise that you don't see the correct result.

The time domain function

$$x(t)=e^{j2\pi f_0t}e^{-Rt}\tag{1}$$

has the following Fourier transform:

$$X(f)=\frac{1}{R}\frac{1-j\frac{2\pi(f-f_0)}{R}}{1+\left(\frac{2\pi(f-f_0)}{R}\right)^2}\tag{2}$$

Also, scaling should be multiplicative and not additive (unless you're in the logarithmic domain), and you should also not only look at the real parts of complex-valued functions.

As explained in this answer, approximating the CTFT by the DFT usually results in two types of errors: the truncation error (due to truncation of the time domain function), and the aliasing error (due to sampling the time domain function). These errors can be made small by choosing a large sampling frequency and a large DFT length.

Below is an example code in Octave/Matlab showing how the CTFT of the given function can be approximated by the DFT:

F0 = 0.3;
R = 0.04;

Fs = 30;            % sampling frequency
Ts = 1/Fs;
Tmax = 100;         % length of time domain signal
N = round(Tmax/Ts); % DFT length

t = (0:N-1) * Tmax / (N-1);
N2 = round(N/2);
f = (0:N2-1)/N*Fs;

fun = exp(1i*2*pi*F0*t) .* exp(-R*t);

% analytic expression for CTFT
z = 2*pi*(f - F0) / R;
FTfun = (1 - 1i*z) ./ (1 + z.^2) / R;

% DFT approximation of CTFT
FTfun2 = Ts * fft(fun,N);
FTfun2 = FTfun2(1:N2);

plot(f,20*log10(abs(FTfun)),f,20*log10(abs(FTfun2)))
    title('magnitudes (dB)'), xlabel('f'), legend('CTFT','DFT'), grid on
    axis([0,Fs/2,-40,30])


The Sum-area table (SAT) was introduced in computer graphics, and as far as I know, was popularized in computer vision and image processing under the name Integral image.

Apart from a larger genericity of the SAT to objects of dimensions different to that of images, the concepts are the same to me.


First-order low pass filter with bounded derivative

I will show you two bounds on the derivative of the output signal, both based on the system's impulse response. First, by taking the Laplace transform of the time-domain equation and rearranging you get the Laplace transform of the output signal in terms of the transform of the input signal and of the system's transfer function:

$$Y(s)=X(s)H(s)=X(s)\frac{1/\alpha}{s+1/\alpha}$$

The inverse Laplace transform of $H(s)$ is the system's impulse response

$$h(t)=\frac{1}{\alpha}e^{-t/\alpha}u(t)$$

where $u(t)$ is the step function. In the following I will assume that $\alpha>0$ is satisfied, because otherwise the system is not stable.

The transform of the derivative $y^{\prime}(t)$ is

$$sY(s)=sX(s)H(s)=[sX(s)]H(s)=X(s)[sH(s)]\tag{1}$$

from which it can be seen that the derivative $y^{\prime}(t)$ can be obtained either from filtering the derivative of the input signal $x^{\prime}(t)$ with $h(t)$, or by filtering $x(t)$ with a system with transfer function $sH(s)$. The impulse response of this system is

$$h^{\prime}(t)=\frac{1}{\alpha}\delta(t)-\frac{1}{\alpha^2}e^{-t/\alpha}u(t)\tag{2}$$

From these two interpretations we can now derive two bounds on $|y^{\prime}(t)|$. First note that the magnitude of the output of a linear time-invariant system can be bounded as follows:

$$|y(t)|=\left|\int_{-\infty}^{\infty}h(\tau)x(t-\tau)d\tau\right|\le
\int_{-\infty}^{\infty}|h(\tau)||x(t-\tau)|d\tau\le
|x(t)|_{max}\int_{-\infty}^{\infty}|h(\tau)|d\tau\tag{3}$$

Since $y^{\prime}(t)$ can be obtained from filtering $x^{\prime}(t)$ with $h(t)$ we get the following bound

$$|y^{\prime}(t)|\le|x^{\prime}(t)|_{max}\int_{-\infty}^{\infty}|h(\tau)|d\tau=
|x^{\prime}(t)|_{max}\frac{1}{\alpha}\int_{0}^{\infty}e^{-\tau/\alpha}d\tau=|x^{\prime}(t)|_{max}\tag{4}$$

This is the first bound on $|y^{\prime}(t)|$, which is of course only useful if one knows an upper bound on $|x^{\prime}(t)|$.

The other bound is obtained by noting that $y^{\prime}(t)$ can also be interpreted as filtering $x(t)$ with the impulse response $h^{\prime}(t)$ given by (2):

$$y^{\prime}(t)=(x*h^{\prime})(t)=\frac{1}{\alpha}x(t)-\frac{1}{\alpha^2}\int_0^{\infty}e^{-\tau/\alpha}x(t-\tau)d\tau$$

This leads to the following bound

$$|y^{\prime}(t)|\le\frac{1}{\alpha}|x(t)|_{max}+\frac{1}{\alpha^2}|x(t)|_{max}\int_0^{\infty}e^{-\tau/\alpha}d\tau=\frac{2}{\alpha}|x(t)|_{max}\tag{5}$$

Combining (4) and (5) we obtain

$$|y^{\prime}(t)|\le\min\left\{|x^{\prime}(t)|_{max},\frac{2}{\alpha}|x(t)|_{max}\right\}$$

Obviously, both bounds depend on properties of the input signal $x(t)$, and only one of them depends on the system parameter $\alpha$. In sum, without knowing anything about the input signal $x(t)$, there is no way to choose $\alpha$ such that $y^{\prime}(t)$ is bounded for all input signals.


How to sketch a frequency spectrum for an AM signal?

Since a signal's spectrum is complex, it can't be drawn in a single two-dimensional plot. What is usually done is to do two plots; one is the magnitude spectrum and the other is the phase spectrum.

When looking at signals, most of the time, the magnitude spectrum contains all the information that is needed and the phase spectrum is ignored. When looking at filters, the phase spectrum is important because one usually wants the phase response to be linear in the filter's passband.

So, what you need is the magnitude of $V(f)$:

$$|V(f)| = 25 \left[ \delta(f - 500k)+ \delta(f+500k)\right] + 5\left[\delta(f-999.5k) + \delta(f+999.5k) \right] + 5\left[\delta(f-500.5k)+\delta(f+500.5k) \right].$$

The magnitude spectrum is real, so it can be drawn in a regular 2-D plot. Note that all values of the magnitude spectrum are positive.


There is no real rule, but you would usually choose around $10N$ frequency points, where $N$ is the desired filter order. Distribute the points equidistantly over pass band(s) and stop band(s). This does not mean that they are equidistant in the interval $[0,\pi]$ because usually there is at least one "don't care" region where no desired response is specified (such as in the transition band between pass band and stop band).


The typical thing to do is the low-rank approximation on separate channels. Assume that $C$ is a channel of the RGB image $I$:

rank = 10;
[U,S,V] = svd(C);
L = U(:,1:rank) * S(1:rank, 1:rank) * V(:, 1:rank)';

Now, L should be the compressed image. If you do this operation and compose the channels back, you should get a compressed RGB image.

However, such a method in my opinion is only good for mathematical understanding. It is not very practical due to the fact that it depends on image dimensions too much and cannot generate high SNR results. Or better stated, as you decrease rank, you drastically reduce the quality of the reconstruction. For this reason, one has to retain around ~50-100 components. Better methods doing this blockwise, or formulating the compression as an energy functional exist. Image denoising is a huge area of its own.

You might want to check some results with varying rank in this presentation.

Also, a similar question has already been asked; please follow the Stack Exchange discussion here.


There are some issues with the code.
1. There seems to be a problem number of taps used for frequency. For filtering 400Hz sine wave at samplerate 48KHz I get better results if I use tapslen more than 559.
2. ippsFIRGenBandstop apparently is implemented as band pass + band stop + band pass. By doing this IPP tends to add some ripple at the higher end (around 44,000 and above). See here:

When using this filter it is recommended to use some kind of a FIR filter design tool like TFilter: http://t-filter.appspot.com/fir/index.html
Or matlab FIR filter design toolbox.


The result of DTW is a matrix $N\times M$, where $N$ is the size of the first MFCC-dataset in comparison and $M$ is the size of the second one. Now to make a decision of how close two datasets are to each other you need to see at the normalized bottom-right value value of the DTW matrix. If it's less than some pre-defined threshold, then you can say that two datasets are of a similar nature. In order to compute the threshold you need to run DTW multiple times and plot a ROC curve. The following sample DTW implementation in Matlab should help:

function cost = dtw (A, B)
  [rowsA colsA] = size (A);
  [rowsB colsB] = size (B);

  D = inf (rowsA + 1, rowsB + 1);
  D(1) = 0;

  for i = 2:rowsA + 1
    for j = 2:rowsB + 1
      d = norm ( A(i - 1, :) - B(j - 1, :) );
      D(i, j) = d + min ([ D(i - 1, j - 1) D(i - 1, j) D(i, j - 1) ]);
    endfor
  endfor

  cost = D(end) / (rowsA + rowsB);
endfunction


Correction of signal after High-pass RC-filter

Yes, this seems like a normal thing to do. However, note that the signal-to-noise ratio will be degraded in the 1.3 Hz to 15 Hz band. Your signal will be smaller but the quantization noise of the ADC will stay the same.


Currently I am writing a script that will look at the time/size = compression rate. Anything that is below 300kb/s I want to delete.

That's a very bad idea.

Just because something compresses well, because it fits the signal model of the compressor well, doesn't mean it's bad quality. Usually, quite the contrary.

You give an example in your comments yourself:

  Im talking about tracks that have lot of clipping

Clipping is a nonlinearity, introducing a multitude of tones into the spectrum that shouldn't be there. Since psychoacoustic models and thus lossy compressors usually need to quantize the spectrum in some way, having clipping in your audio means harder to compress audio meaning larger file size.

Generally, if something is hard to compress, a variable-rate codec will increase the file size. If something is easy to compress without much loss, then the file size will stay down. I'm pretty sure you can record a pretty perfect tuning fork tone in high quality and compress it very nicely.

Also, unless we're talking early MP3 encoders and such, modern codecs above ca. 160 kb/s are really hard for a human to tell from the uncompressed audio. That's the whole idea of having a good lossy compressor: Just lose the information that's irrelevant.

Note that I'm not saying you can't generally say that for example an MP3 with 64 kb/s will sound bad in almost all cases. It's just that you don't need any intelligence to do that, but can often directly read that from the file's metadata.

  Like looking at a spectral analysis to determine what is low vs high quality.

Nope, you're looking at the reconstructed signal after decompression. Unless you have extensive knowledge of the spectrum of your original audio, there's little information about audible losses in that. Again, both musicians, sound engineers and audio codecs strive to do the same: produce sound that sounds good to the human ear.

Also, remember that clipping (or soft clipping, and other similar nonlinearities), as the example you chose, is actually a pretty commonly used tool in production of music these days. Who are you to tell someone that the recording of their neighbor's cat walking on a keyboard is quality-wise better than let's say a Kraftwerk song that uses clipping intentionally? How is a crisp-sounding KPop-girlband song with strongly emphasized treble quality-wise superior to the muddy sound of the 1968 original release of Jimi Hendrix' All along the watchtower? There's whole musical genres that depend on sound mechanics that would signify a lack of recording quality in others. Even within the same genre, quality can't simply be deducted from a rigid description of the signal – compare Shostakovich's The Bolt to other ballet suites of that era.

"Good audio quality" in the end is an extremely subjective matter, and I'm afraid without a large-scale database of manually tagged "good" and "bad" examples according to your perception, you won't be able to implement any classifier, be it through neural nets or through more classic statistical approaches.


Amplitude Estimating Using a Windowed DFT

I have several blog articles that solve this problem exactly in a theoretical sense, and quite accurately in an implementattion.

3. DFT Pure Tone Frequency Formulas

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

4. DFT Phase and Amplitude

Phase and Amplitude Calculation for a Pure Real Tone in a DFT: Method 1
Phase and Amplitude Calculation for a Pure Complex Tone in a DFT
Phase and Amplitude Calculation for a Pure Complex Tone in a DFT using Multiple Bins

The rest of the articles are pretty good too. ;-)

I posted one this morning:

Angle Addition Formulas from Euler's Formula


$h_0(t)$ is inverse Fourier transform of $H(jw)$. Your formula is OK, you can continue your calculation to practice your math manipulation, why not.
To check the result, you can remark that $H(jw)$ is rotated $sinc$ function. And $sinc$ must be the Fourier transform of a rectangular box function.


In my experience, sound card sampling clocks vary by as much as 10 Hz from one card to another, even between computers of the same make and model. I have never seen much jitter, but I haven't especifically looked for it.

Note, however, that it is essentially impossible to know the true sampling rate, unless you have access to an atomic clock. What you can do is approximate the difference relative to a clock that you designate as your reference.

A simple experiment is to generate a single-tone signal, then receive it and find its frequency. Say your reference is the transmitter and it's operating at 48,000 samples per second. Feed the transmitter a sequence 0, 1, 0, -1, 0, .... This will produce a sine wave of nominal frequency 12,000 Hz. In the receiver, you sample this signal for a relatively long time (to have good frequency resolution) and use the DFT to determine its frequency to be, say, $f_2$. Then, you can calculate $$\Delta=\frac{12000}{f_2}$$ which is the factor by which the receiver's clock is faster than the transmitter's. Of course, this is only one of many possible approaches.


Channel weighting of Mono channel in EBU R128

Use a channel weight of 1.0 for your mono signal.

The weighting of 1.5 presented in ITU-R BS.1770 addresses the fact that a signal coming from behind is perceived louder by a human than if it is coming from the front (Reason: To better detect an enemy sneaking at you from behind). With a single mono signal this fact doesn't need to be addressed.

The 3 LU difference you mentioned also makes perfect sense from a physical point of view. One speaker in front of you compared to two (equal) speakers only provides half the energy and therefore the SPL will be lower by 3dB.

I'm honored that you use my LUFS Meter to compare to, but I'd recommend to also use the test files from ITU and EBU. The LUFS Meter manual contains the links to these test files.

Are you aware that the sources of the measurement algorithm of the LUFS Meter as well as some development notes are available? (Sadly I'm not allowed to post a third link... They are at docs/developmentNotes/)

Disclaimer: I'm the developer at Klangfreund.


What does it mean for a signal to alias to another?

Sampling of a signal leads to repetition in the frequency domain. Typically, the frequency domain is considered from $-\pi/2$ up to $\pi/2$ or from $0$ to $\pi$, where $\pi$ corresponds to the sampling frequency.

For a real-valued signal, the spectrum is also symmetric around $n\pi/2$ for any integer $n$. For your example where a sine wave with frequency $7$kHz is sampled at $6$ kHz, this means that the sampled signal contains the frequencies $1+6\cdot n$ kHz and the symmetric/mirrored frequencies $5+6\cdot n$ kHz, e.g., 1 kHz, 7 kHz, 13 kHz and 5 kHz, 11 kHz (and also the negative frequency parts.

You speak of 'aliasing over another signal' when the aliases overlap with other signal content. For example, when you have a signal consisting of the sine wave of 2 kHz and 7 kHz, which you sample with a rate of 6 kHz, then you observe the frequencies 2 kHz; 4 kHz and 1 kHz; 5 kHz. These frequencies do not overlap. When the original signal contains 1 kHz and 7 kHz and you sample at a rate of 6 kHz, then the '7 kHz signal aliases over the 1 kHz signal'.

Of course this principle works also with other signals, such as narrowband, low passed and high passed signals.


FFT subset output yields computational efficiency

Assuming your question is about bypassing the later stages of the FFT if not all samples are needed, this may give you some insight:

Knowing that the FFT alorithm is a successive decimations into even and odd samples so that smaller FFT's can be computed and thus leading to the great efficiency of the algorithm; here is a view of what the digital spectrums look like prior to combining in the final stages:

Everything between the pciture above showing the FFT decimation and the picture below showing the final FFT combining is covered in this post, so I won't repeat it here, but please review it if the plot below is not immediately clear: Is this signal perfect reconstructable?

From reading the linked post and the plots shown above, it will hopefully be clear that unless the signal is combined in the final stage, the spectrum will be superimposed with other regions of the spectrum.  So if we work back one stage from the final output, as shown in the figure below, we can hopefully associate with the figure above to see that the upper N/2 pt DFT contains ALL the spectrum (folded on itself), and the lower N/2 pt DFT also contains ALL the spectrum, however separable due to the phase relationship described if combined with the appropriate phase rotations as shown.

So that said, we can see the reason for the combining (and that the signal would be irrecoverable if we did not have access to both prior DFT's). However as long as both DFT's are available for a given point of interest, we can use the mapped output from both DFT's for that point without calculating the rest of the final stage. This of course works back in the same pattern through the implementation of the smaller DFTs (since they two are made up of two DFTs, which are each made up of 2 DFT's etc.

We can see from this plot how the final stage works its way back through these successive DFT's to the input. This plot also gives immediate insight into the possible reduction if only a subset of the points were needed (and shows how quickly any savings if diminished- in this simple case, 1 point is touching 21 of the 48 nodes shown!).  I think it is also interesting (perhaps a little mind-blowing) that if you work your way back to the first stage, each 2 pt DFT output for the same reasons described also contains ALL the spectrum, but in this case folded over itself multiple times (depending on how many stages).  Hence every 2 pt DFT at the input is connected to every output in the final DFT.


MATLAB Rayleigh fading and white gaussian noise

The Rayleigh fading channel equation you provided comes from the property that given two independent zero-mean Gaussian random variables with equal variance $X \sim N(0,\sigma^2)$ and $Y \sim N(0,\sigma^2)$, the random variable $R = \sqrt{X^2+Y^2}$ is Rayleigh distributed (see for example wikipedia). In the code you provided, the real and imaginary components (generated by independent calls to randn) generally meet this condition (or at least approximates it quite well for reasonable pseudo-random generator) and the magnitude of h would thus have a Rayleigh distribution.

In addition, it is generally assumed that the signals' power is preserved on average, that is $E\{R^2\} = 1$. Now given how we defined $R$:
$$
\begin{align}
  E\{R^2\} &= E\{X^2+Y^2\} \\
           &= E\{X^2\}+E\{Y^2\} &\mbox{(linearity of expectation)} \\
           &= 2E\{X^2\} &\mbox{($X$ and $Y$ identically distributed)} \\
           &= 2\sigma^2 &\mbox{($X$ and $Y$ zero-mean Gaussian)}
\end{align}
$$

So the signal power is preserved on average when $\sigma = 1/\sqrt{2}$. This is also the scaling factor that must be used if starting with unit variance Gaussian pseudo-random variables (as is the case with MATLAB's randn).

Similarly, given a complex-valued Gaussian noise $n$ defined as $n = n_{\scriptsize \mbox{real}} + i\cdot n_{\scriptsize \mbox{imag}}$, where both real and imaginary components are Gaussian distributed with variance $\sigma_n^2$, the noise power is
$E\{|n|^2\} = E\{n_{\scriptsize \mbox{real}}^2 + n_{\scriptsize \mbox{imag}}^2\}$. By a similar argument as above, it follows that $E\{|n|^2\} = 2\sigma_n^2$. To obtain noise with unitary power, we thus need $\sigma_n = 1/\sqrt{2}$.

Note that the notion of unitary power only makes sense once a measurement unit has been established. Quite often the noise power is defined relative to the signal power. In that case, a noise power of 1 (or 0dB) in signal power's units would mean that the power of the noise is the same as that of the signal.


Any simple Digital All-Pass Filter Design Tools?

This is an excerpt from the FilterScript user's guide, which should provide you with enough information on how to design the filter. Please refer to the other product resources on the tool's homepage for more information.


Frequency response of system with feedback and feedforward co-efficients?

Yes, you can always use the FFT to compute the frequency response of a discrete-time system at discrete frequencies. For an IIR filter you get the following transfer function:

$$H(z)=\frac{\sum_{m=0}^{M}b_mz^{-m}}{\sum_{n=0}^{N}a_nz^{-n}}\tag{1}$$

Usually, $a_0$ is normalized to 1. In order to evaluate $H(z)$ on the unit circle $z=e^{j\theta}$ you can compute the FFTs of the numerator and of the denominator of (1) separately and then divide to get the frequency response of the filter.

If you use Matlab or Octave, you can simply use the command $\tt{freqz}$.


Basic question: Why is the output of a system the convolution between the impulse response and the input?

Because it is the response of the system when an unit impulse (delta function) is applied to the input of the system. So basically, If you  multiply this output impulse response of the system to each of the input samples and then add them all, you get the overall output of the system. And this is only true, If the system agrees with the linearity properties such as homogeneity, additivity and shift invariance. For more details and to understand the mathematics behind this, you can go through this link: http://www.dspguide.com/ch6.htm

I tried to put this as comment instead of answer, since its not detailed enough, however, due to my lack of reputation, I could not comment. Anyways, the link I have provided is very useful. Hope it helps!


How is Point Spread Function (PSF) related to Image Priors in Blind Deconvolution?

A "point spread function" denotes how a single point from the source, or a punctual object (mathematically a Dirac) would spread on the observed image through the imaging system.

In blind deconvolution, both image priors and PSF priors are useful to unmix those intricate systems (even more with noise). They are different, as the PSF is somehow invariant to the imaging system (there exist space variant PSFs or blurs, let us skip that for the moment), while different images can be acquired.

PSF priors may be related to some parametric shape (e.g. a Gaussian of unknown width), a bounded support in space or frequency, smoothness, etc, while data priors are often more generic to some image class, related to sparsity, bounded variation, etc.


Look into the MUSAN corpus from the Kaldi people (Dan Povey, et al):

Corpus audio/text
Check the arxiv paper

  MUSAN is a new corpus of music, speech, and noise. This dataset is
  suitable for training models for voice activity detection (VAD) and
  music/speech discrimination. Our corpus is released under a flexible
  Creative Commons license. The dataset consists of music from several
  genres, speech from twelve languages, and a wide assortment of
  technical and non-technical noises. We demonstrate use of this corpus
  for music/speech discrimination on Broadcast news and VAD for speaker
  identification.


Most reliable way to decode PAM2 signal with unknown signal levels?

I suggest as a simple approach start with the threshold you chose and then bin the samples that are above and below threshold and take the average of each of those separately and then adjust your threshold to be between those two values. Mid way would be optimum when the two symbols are equiprobable in the presence of identical balanced noise on all samples (meaning the distribution going positive is the same as the distribution going negative). In the case on more 1’s than 0’s or vice versa with known or computed SNR, the threshold decision would be determined by a weighted average of the two symbol means using the SNR’s to be closer to the symbol that occurs less likely. (To optimize bit error rate).

Given the OP's comment under the question of using complex IQ samples with Gaussian noise, the 0 bit would be Rayleigh distributed and the 1 bit would be Ricean. The approach above could determine the estimated variance for each of the symbols (and therefore SNR from mean and variance) and then if the probability of occurrence is known for each symbol, a threshold can be determined that would minimize the error rate by equating the net areas of each distribution beyond the threshold weighted by the probability of occurrence.


How to calculate EVM in %age of an Equalized Constellation in 16QAM?

The Error Vector is the Euclidean distance from the actual sample at the optimum timing location in each symbol to the actual symbol location in a reference constellation (as the distance to closest decision boundary, just prior to decision). The measurement metric Error Vector Magnitude (EVM) is computed as an rms quantity over multiple error vectors, where the waveform and the decision boundaries are scaled to either the rms amplitude or a maximum amplitude in the constellation, depending on which standard is used. We can therefore use one sample for every symbol to compute a sufficient number of errors to derive a statistical quantity.
Below shows a single Error Vector for 16 QAM, and the error vector magnitude would be the root-mean-square of all error vectors as given by the formula above (once the constellation and the received signal are properly normalized as described).  It also must be mentioned that the normalization by the rms of the waveform assumes the symbols in the waveform are equally distributed, but if not normalization can occur through known training sequences. The idea is that the error vector magnitude should be zero if all the symbols land on the exact locations in the constellation.

To compute EVM normalize the constellation (by either the rms, the peak level of an outermost point, or the rms of the magnitudes depending on which definition is used). The waveform samples just prior to decision (one sample per symbol) are similarly normalized which would then minimize the computed EVM (this is consistent with what the receiver would do anyway just prior to decision to minimize error rate, so it is very easy to compute).
Once properly scaled the EVM computation in its general form is just the formula for the standard deviation of samples from a zero mean random process and given as a percentage quantity; so the average error magnitude given as a percentage of a normalized constellation:
$$EVM = \frac{\sqrt{\frac{\sum_n|t_n-\hat t_n|^2}{N}}}{S} \times 100\%$$
(Where $S$ is given by whatever normalization is used for the constellation)
Notice at Wikipedia en.wikipedia.org/wiki/Error_vector_magnitude it is defined as "The ideal signal amplitude reference can either be the maximum ideal signal amplitude of the constellation, or it can be the root mean square (RMS) average amplitude of all possible ideal signal amplitude values in the constellation. ". In this paper they define the normalization as the rms of the constellation: eprints.soton.ac.uk/263112/1/paper_101.pdf , while here the outermost point is used: ieee802.org/16/tg1/contrib/802161c-01_28.pdf .
Point being, if an EVM computation result is being provided or used, it is important to provide the definition and equation that is being used for that particular computation.


Question regarding estimation of signal tone parameters (frequency, amplitude, and phase) using Macleods algorithm

I'm not sure about Macleod, but generally windowing is not good for that style of estimator. See this paper for details about similar estimators.

Generally the more data, the better for all parameters. Provided you use an estimator without bias, it should do better with more data.

It depends on what  you mean by "resolution".

If you mean the standard deviation of the estimate, then that will be governed (at higher SNRs) by the Cramer-Rao lower bound for frequency estimate $\hat{f}_0$, which is:
$$
{\rm var}(\hat{f}_0) \ge \frac{12}{(2\pi)^2 \eta N(N^2-1)}
$$
where the signal model is
$$
x[n] = A \cos(2\pi f_0 n + \phi) + w[n]
$$
the amplitude is $A$ and phase $\phi$ (which are also both unknown), the number of samples is $N$, $\eta = A^2/(2\sigma^2)$ is the SNR and $w[n]$ is white, Gaussian noise of variance $\sigma^2$.

See equations (3.41) of Kay's Fundamentals of Statistical Signal Processing, Volume 1: Estimation Theory.


According to this paper, the grayscale dilation of image $I$ by a non-flat structuring element $S$ is defined as follows:

$[I ⊕ S](x,y) = \displaystyle\max_{(s,t)∈S}\{I(x-s,y-t)+S(s,t)\}$

Since the origin of the structuring element is in the top right corner, we have that $s∈[0,1]$ and $t∈[-2,0]$ excluding the point $(s,t)=(0,-1)$. So using your example of element $I(0,1)$ we have the following:

$[I ⊕ S](0,1)=\max\{I(0,3)+S(0,-2),\ I(-1,3)+S(1,-2),\ I(-1,2)+S(1,-1),\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ I(0,1)+S(0,0),\ I(-1,1)+S(1,0)\}$

$=\max\{0+0,0+0,0+2,2+0,0+(-1)\}=\max\{0,0,2,2,-1\}=2$

Which is the correct value for element $(0,1)$ in the dilated matrix. You can then simply apply the formula to the other $(x,y)$ value pairs of interest to obtain the final result.

The definition of erosion (and the other morphological operations) can also be found in that paper and applied in the same manner.


How to Design a Pre-emphasis filter knowing the gain, a zero and a pole?

Traditionally, these filters were implemented as analog filters, and they are different from the usual discrete-time first-order FIR pre-emphasis filter used in speech processing.

The transfer function of such an analog pre-emphasis filter is

$$H(s)=g\frac{1+s\tau_1}{1+s\tau_2},\qquad\tau_1\gt\tau_2,\tag{1}$$

where $g$ is the DC gain, and $\tau_1$ and $\tau_2$ are two time constants determining the locations of the zero and the pole:

$$s_0=-\frac{1}{\tau_1}\\s_{\infty}=-\frac{1}{\tau_2}\tag{2}$$

If the pole (or the zero) are given in terms of frequencies $f_i$, then the following relation holds:

$$\tau_i=\frac{1}{2\pi f_i},\qquad i\in\{1,2\}\tag{3}$$

So for a given frequency $f_i$ you can use $(3)$ to determine the corresponding time constant, from which you obtain the desired transfer function $(1)$.

If you want to implement that filter in discrete time, the most common option is to use the bilinear transform:

$$s=\frac{2}{T}\frac{z-1}{z+1}\tag{4}$$

where $T$ is the sampling period, i.e., the inverse of the sampling frequency. Note that $(2)$ is the bilinear transform without pre-warping, which means that the frequency axis will be warped, so neither the pole frequency nor the zero frequency of the analog filter will be exactly realized in the discrete-time domain. You can pre-warp the analog frequencies (or time constants) to make sure that the resulting discrete-time filter has the desired pole and zero frequencies. In order to realize a (pole or zero) frequency $f_i$ in the discrete-time domain, you need to use the following time constant in the analog filter transfer function $(1)$:

$$\tau_i=\frac{T}{2\tan\left(\frac{\pi f_i}{f_s}\right)}\tag{5}$$

where $f_s=1/T$ is the sampling frequency.

You can find more information on pre-warping in this answer.


The primary problem in your time domain analysis is the assumption that if both $x(t)$ and $x(t-d)$ have the same Nyquist sampling rate so will $x(t)+x(t-d)$.

You can see that the summation may alter the signal in such a way that resulting signal's bandwidth may be different than the individual ones.

As an example consider the composite signal $x(t) = x_a(t) + x_b(t)$ where $x_a(t)$ is a low-pass signal with individual Nyquist sampling rate of $w_a$ and $x_b(t) = \sin(w_0 t)$, is a high frequency sine wave with a Nyquist sampling rate of $2w_0$. Therefore both $x(t)$ and $x(t-d)$ will have the same Nyquist sampling rate of $2 w_0$.

Then we have $x(t-d) = x_a(t-d) + \sin(w_0(t-d)) = x_a(t-d) + \sin(w_0 t - w_0 d)$. Therefore the sum of those two composite signals, $x(t) + x(t-d)$, will be $$x(t)+x(t-d) = x_a(t) + x_a(t-d) + \sin(w_0 t) + \sin(w_0 t - w_0 d)$$

Now if $d$ is chosen such that $w_0 d = \pi$ then $\sin(w_0 t) + \sin(w_0 t - \pi) = 0$, for all $t$, therefore $x(t) + x(t-d) = x_a(t) + x_a(t-d)$. There will be cancellation of the high frequency terms so that the Nyquist sampling rate of the sum will be $w_a$ (assuming the sum of low-pass signals will not further alter it). Hence the sum will have a different Nyquist sampling rate eventhough $x(t)$ and $x(t-d)$ individually have the same Nyquist sampling rate.


trying to generate a sweep sine wave that increases its frequency till an upper limit and then starts decreasing its frequency to the lower limit

I think that in your simulation, you only calculate the instantaneous frequency, but then fail to compute the phase for the sin() function correctly.

This is how you could correct this issue:

% parameters of linear up/down-sweep
freq1 = 20;
freq2 = 200;
fs = 44100;

% length of a sweep in one direction
endTime = 0.1;

% number of samples of a sweep in one sweep direction
N = endTime * fs;

% instantaneous frequency at each point in time ...
% ... first increasing for N points, then decreasing for N points
inst_f = [linspace(freq1, freq2, N) linspace(freq2, freq1, N)];

% since (in continuous time) instantaneous frequency is derivative of
% phase, we have to compute the "integral" to get the phase for sin().
phi = 2 * pi * cumsum(inst_f) / fs;
sweep = sin(phi);

% plotting
t = linspace(0, endTime * 2, length(s));
plot(t, sweep);

% wavwrite(sweep, fs, 32, '~/sweep_delete.wav');


What's the shortest possible code that correct all the single bit errors?

In order to correct 1-bit errors, we need to have a Hamming distance (the number of bits that must be flipped) of 3 between all the symbols: when we flip one bit at random, then flipping another bit by random must either get us back to the original symbol (if the same bit was flipped twice) or an invalid symbol.
With only two symbols we can encode as 000 and 111 because this gives us the required Hamming distance, so clearly we can encode 3 symbols with 6 bits. Can we can do better though?
If we use 4 bits to encode one symbol as 0000 then 3 flips result in 0111, 1011, 1101 or 1110. These all differ from each other by no more than 2 bits, so clearly 4 bits is not enough.
On the other hand, if we use 5 bits to encode one symbol as 00000 then 3 flips give us 00111, 01011, 01101, 01110, 11100, 11010 and 11001. If we then encode the 2nd symbol as 00111, we're left with the following valid choices for the 3rd symbol: 11100, 11010 or 11001. These are all distance 3 from the first symbol and distance 4 from the second symbol.
So one of the many possible encodings for 3 symbols using 5 bits is:
a = 00000,
b = 00111,
c = 11100.


It's more like a soft answer (I am happy to update it later), but Alex Acero explained the technology behind aniomoji on this years ICASSP 2018. Here is the link. Basically, they are using so-called Blendshapes to estimate the positions of the muscles and then they are mapping them to the actual animoji model.


The actual frequency in Hertz corresponding to DFT index $i$ with a DFT length of $N$ is

$$f=\frac{i}{N}f_s,\quad i=0,1,\ldots N-1\tag{1}$$

where $f_s$ is the sampling frequency. So if indices are numbered starting from $0$ (and index number $19$ is indeed the $20^{th}$ frequency bin, as assumed in (1)), then your answer is correct.

As for part 2 of the question, the answer is also given by (1), so

$$f_s=\frac{Nf}{i}$$

which gives $f_s=289.7\text{ kHz}$ for the given values of $i$, $N$, and $f$.


What does the word 'bin' mean in the context of FFT's?

A finite length DFT transform has a finite number of basis vectors, which are all single frequency complex sinusoids.  If the results weren't binned then you would need an infinitely long DFT or FFT result to represent any continuous spectrum.  However, Parseval's theorem says that all the energy in a finite input vector ends up in the finite FFT result vector.  So any single between-bin-center-in-frequency pure sinusoid (or between DFT basis vector frequency, windowed by any finite length FFT) ends up binned, spread mostly into FFT result bins very near in frequency.

"very near" because:

Each FFT result bin is not rectangular, but a filter in the shape of Sinc function (more precisely, a periodic Sinc function or Dirichlet kernel), since that is how a pure sinusoid correlates with any one DFT complex basis vector when swept across the frequency range from 0 to Fs/2.  The result is also conjugate mirrored in the full length FFT result for strictly real input.

Another interpretation is that there exist no in-between bin frequencies, that any cut-off sinusoid is really a repetitively discontinuous periodic waveform.  But that mythical input is not representative of how FFTs are actually used in many forms of DSP signal analysis (windowing, peak estimation, low offset overlapping windows, phase vocoders, etc.). Although knowing about this circularity is important when using an FFT for fast convolution.


Computing of speech signals modulation entropy

Ok, I've found what i was looking for.

This paper describes how to measure entropy in signals:

Entropy Based Voice Activity Detection in Very Noisy Conditions, by Philippe Renevey and Andrzej Drygajlo. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.123.6098&rep=rep1&type=pdf

And this question covers the part of coding side of question:
https://stackoverflow.com/questions/18962990/calculate-entropy-of-a-file


Differentiation of discrete sinusoidal signal

Hilbert transform has the property of harmonic correspondence (Vakman 1996), meaning that the Hilbert transform of a cosine wave is a sine wave etc. So the Hilbert transform is a 90 degree phase shift. In fact the Hilbert transform in the Fourier domain looks quite similar to the differential operator:

$$\mathcal{F}\left[\tfrac{d}{dt} f(t)\right](u) = -iw \mathcal{F}[f(t)](u)$$
$$\mathcal{F}\left[\mathcal{H}[f](t)\right](u) = -i\frac{w}{|w|} \mathcal{F}[f(t)](u)$$

We can see that for a sinusoid of frequency $w_0$, the derivative is the Hilbert transform scaled by $w_0$. e.g. Hilbert transform of $\cos(n t) = -\sin(n t)$, derivative of $\cos(n t) = -n \sin(nt)$. So careful choice of alpha in your formula will give you the phase shift, while compensating for the scaling by frequency, which is the same as the Hilbert transform.

If you have a non-sinuosidal signal, then the 90 degree phase shift must be applied to every sinusoidal component. This is not the same as differentiation, and you couldn't use your formula. e.g. Hilbert transform of square wave (from wikipedia):


How to process video with DFT/FFT - getting the frequency of green color changes?

If you need to track changes in time, then the most natural and obvious way is to form one temporal (green data) vector per pixel that you want to observe, then take its DFT (it's way 2, but with more than 2 frames).

The actual vector length and DFT function that you should use will depend on the actual case you are working on.

Suppose you work from a recorder video. Then, you can take the 1D DFT of the vector formed with the value of one pixel of interest in the whole sequence.
On the other hand, if you work from a live feed, then you need to use a FIFO data vector and probably apply a windowed 1D Fourier transform to avoid artifacts, see short-time Fourier transform.


Inter Symbol Interference (ISI) / Equalization in Digital Signal Processing

If the channel just adds noise and does not cause any distortion (as in the first figure), you can design the transmit filter $g(t)$ such that, when combined with its matched filter at the receiver, the Nyquist criterion for zero intersymbol interference (ISI) is satisfied. If the channel adds linear amplitude and phase distortion (as in the second figure), then generally the Nyquist criterion will not be satisfied and you will observe ISI, which will increase the error probability.

A linear equalizer is designed to (at least partly) compensate the effect of the channel $h(t)$. There are different optimality criteria, such as the zero-forcing criterion you mentioned, which eliminates ISI completely, but which does not take into account the noise. A zero-forcing equalizer basically inverts the channel response, if this is possible at all (with a causal and stable filter). If the noise is not negligible, it is usually better to use a different criterion such as the minimum mean square error (MMSE) criterion, which jointly minimizes the effects of ISI and noise.

Since many channel responses are unknown and/or time-varying, i.e., the impulse response $h(t)$ changes over time, equalizers must be adaptive. In practice, most adaptive equalizers use a sort of least mean squares (LMS) adaptation algorithm because it is relatively simple and robust.

This topic is quite broad and complex, so I would recommend that you read a good introductory text on digital communications, instead of gathering unrelated pieces of information from the internet.

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