How to achieve a smoothed frequency phase response with no outliers I have two sampled signals, which I then perform wavelet analysis on, by synthesising the Filterbank, perform the dft on each kernel, and the dft on the sampled signals, and multiply together.

Using this I can get successfully extract the magnitude of the signal, and either take an average or identify the time response of a frequency.

Now I would like to extract the phase.  Taking either an instance in time, or averaging the mean signal I can get the phase and plot it, using atan2(im / re) for each signal, and then subtracting one from the other. However it appears to be very noisy.  Although I can identify the phase trend, there are a lot of sporadic results that lie far out of the trend.

Researching into this a little, I have found there can be issues with rounding small floats, and it is good to set a threshold for the phase, so everything below the threshold is rounded to 0.  This returns similar results.

What I am witnessing looks like some kind of phase noise, reading this it is my understanding this occurs when the frequency is not an exact integer number of cycles for the dft size. This is where I start to get confused as I have my Filterbank frequencies and I have the dft size. How do these two relate?   Would I need to remove all results where the centre frequency of each filter does not meet an integer number of cycles for the sample length?

Am I on the right path here, or is there something else that would be causing my phase noise.  Failing this, is there any decent way to filter out outliers in a signal.

If I go about not performing Wavelet analysis, and do a straight forward dft on an incoming signal and extract the phase, I get similar results, so I think phase noise is along the right path,  please point me in the right direction!

If your signals are not perfectly integer periodic in the FFT aperture (width), the phase measurements of those narrow band signals will toggle between FFT result bins.  This phase toggling will look like phase noise, but is actually a windowing artifact.  It is due to any window edge discontinuity from any non-exactly-integer-in-aperture periodicity in the signal (plus the fact that all the basis vectors of a DFT are perfectly circular in aperture).

To stop this phase toggling (apparent phase noise), perform an FFTShift before the FFT.  The fftshift will rotate the edge discontinuity away from sample 0, and reference the DFT phase estimate to the center of your data window, where it is likely continuous (instead of to a phase discontinuity between the last sample and first sample of the data).

Added:

There are (at least) 2 ways to perform a time-domain fftshift.  The first is to circularly rotate the data vector 50% after windowing but before performing the FFT.  The 2nd way is to complex conjugate every other result bin (odd ones usually) after the FFT.  Due to the DFT shift property, these two methods produce identical results.  Note that these are in the opposite order from doing a frequency domain fftshift.


How do MATLAB and/or Python treat $2^n$ samples rule in FFT As far as I have read, an FFT requires that the number of original data points must be a power of 2.

I'm wondering whether the tools like MATLAB or Python which have FFT functions take care of this fact.

Especially in MATLAB and Python; imagine I have samples of data in a file containing 13000 samples. If I want to read this file and plot its FFT, how does the FFT functions of MATLAB or Python treat the file? Does the function truncates the samples and obtain the number of samples as $2^n$ before performing the algorithm? Or should I sample it myself as $2^n$ samples at the beginning?

Matlab doesn’t require the number of samples to be a power of 2. You can read more about it in Matlab FFT help page.
It will handle any length of signal and will not truncate or pad the length unless you specifically say so.
BTW in the help page it says that Matlab FFT is based on FFTW. In the FFTW documentation it says:

  FFTW is best at handling sizes of the form 2^a 3^b 5^c 7^d 11^e 13^f, where e+f is either 0 or 1, and the other exponents are arbitrary. Other sizes are computed by means of a slow, general-purpose algorithm (which nevertheless retains O(n log n) performance even for prime sizes).


Scaling Property in DFT In below image, we have scaling property of DFT, how the final equation is obtained from the above equation. That is how we are getting the scaling factor , $ \frac{1}{|ab|} $ in the final equation ? Please explain how we are getting $ \frac{1}{|ab|} $. Thank you

Equation 4.12 is given as

$\sum\limits_{m=0}^{M-1} \sum\limits_{n=0}^{N-1} f(am,bn)e^{-j \cdot 2\pi(\frac {k}{M} \cdot m + \frac {l}{N} \cdot n)} $  (let $m'=am, n'=bn$)

$=\sum\limits_{m'=0}^{(M-1)a} \sum\limits_{n'=0}^{(N-1)b} f(m',n')e^{-j \cdot 2\pi(\frac {k/a}{M} \cdot m' + \frac {l/b}{N} \cdot n')}$

($m'$, $n'$ with the incremental interval of $a$ and $b$ respectively)

Note that
$\sum\limits_{m'=0}^{M-1} \sum\limits_{n'=0}^{N-1} f(m',n')e^{-j \cdot 2\pi(\frac {k/a}{M} \cdot m' + \frac {l/b}{N} \cdot n')}$
$=F(k/a,l/b)$

Here $m'$, $n'$ with the incremental interval of $1$, and the points $m$ and $n$ are evenly distributed on the circle, according to the symmetry, you have

$\sum\limits_{m'=0}^{(M-1)a} \sum\limits_{n'=0}^{(N-1)b} f(m',n')e^{-j \cdot 2\pi(\frac {k/a}{M} \cdot m' + \frac {l/b}{N} \cdot n')} = \frac {1}{|ab|}F(k/a,l/b)$


BIBO Stability in Z-domain I'd really appreciate it if someone could please explain to me the condition for a LTI system to be BIBO stable, in z-domain.

I have a background in control, and in linear control for example, if we have the transfer function, we simply look at the poles. But in digital signal processing, sometime we look at the poles - they should be inside the unit circle in order to have BIBO stability - and soemtimes the professor says it depends on the ROC - ROC should include the unit circle -.

I am really confused. Could someone please explain to me the difference?

There is no principal difference between continuous-time and discrete-time systems when judging stability. The imaginary axis of the $s$-plane corresponds to the unit circle in the $z$-plane, and the region inside the unit circle in the $z$-plane corresponds to the left half-plane of the the $s$-plane. In both cases the region of convergence (ROC) is essential. Given only the pole locations, you generally cannot say whether a system is stable or not, unless the poles are on the imaginary axis (the unit circle), in which case the system is unstable. In the general case you need additional information, and that information is given by the ROC. Information about the ROC can also be given in the form of a causality constraint, which basically tells you that the ROC is a right half-plane, or - in discrete time - the region outside the pole with the largest radius.

As explained above, the ROC is essential for deciding if the system is stable or not. As you've mentioned in your question, if the ROC contains the imaginary axis (the unit circle) then the system is stable. The ROC of causal systems is to the right of the right-most pole (or outside the pole with the largest radius in discrete time). That's why all poles of a causal and stable system must lie in the left half-plane (inside the unit circle). For an anti-causal (left-sided) system, the opposite is true. For (non-causal) systems with two-sided impulse responses, the ROC is a strip (annulus) in the complex plane, limited by two or more poles.

It is important to understand that pole locations do not uniquely define a system. That's why they are not sufficient to decide if the system is stable or not. Take as an example a discrete-time system with two poles, one outside the unit circle and one inside. There are three systems (impulse responses) with these same poles. One causal unstable system, one anti-causal unstable system, and one non-causal (two-sided) stable system. Only the two-sided system is stable because its ROC is between the two poles and it includes the unit circle. The ROCs of the other two systems do not include the unit circle.

Also take a look at this answer for a concrete example of three impulse responses corresponding to the same transfer function with different ROCs.


Applying a filter on an audio signal with Python Let's say I have a filter described by its transfer function:

$$
H(\omega) = \frac{1}{1 + j\frac{\omega}{\omega_0}}
$$

And I want to apply this filter to an audio signal (a .wav file) using Python. My initial idea was this:

Split the signal into fixed-size buffers of ~5000 samples each
For each buffer, compute its Fourier transform using numpy.fft.rfft
Apply my filter to the coefficients of the Fourier transform: ft[i] *= H(freq[i])
Compute the inverse Fourier transform of the result using numpy.fft.irfft
Join the buffers together into the output file

It seems like I am doing something wrong, because if I set H to a constant factor ($\neq 1$), then the sound gets completely noisy and unrecognizable. I also tried with different values for the buffer's size.

Here is my code:

import wave
import numpy, numpy.fft as fft
from math import *

inputFile = wave.open('input.wav', 'r')
outputFile = wave.open('output.wav', 'w')

params = inputFile.getparams()
sampleRate = params.framerate
bufferSize = 5000

outputFile.setparams(params)
samples = []

for i in range(params.nframes//bufferSize):
    signal = inputFile.readframes(bufferSize)
    ft = fft.rfft(signal)
    for i in range(1, ft.shape[0]):
        ft[i] *= 0.5 # This does not do what I want...
    ift = fft.irfft(ft)
    samples += list(map(lambda x: max(min(int(x), 255), 0), numpy.abs(ift)))

outputFile.setnframes(len(samples)//params.nchannels)
outputFile.writeframes(bytes(samples))

Thanks for your time and patience :)

Apart from what Marcus suggests in the comments, ft[i] *= H(freq[i]) is probably wrong.  The overlap-add and overlap-save methods should deal with this correctly, but read on!

The FFT does not implement linear convolution. It implements circular convolution, which can approximate linear convolution provided certain conditions are met. One of those conditions is that $$N_{\tt fft} \ge N_{\tt filter} + N_{\tt signal} - 1.$$

In your case, you've effectively set $$N_{\tt fft} = N_{\tt filter} = N_{\tt signal}$$ which will introduce time-aliasing --- the problem that happens when the "certain conditions" are not met.  See here for an answer to a similar question.


finding maximum frequency Assume we have:
$$x_c(t)=\sum_{k=-9}^{9} a_k e^{j(2000\pi kt)}$$

and Fourier coefficients out of $$|k|>9$$  is zero. So, it is band-limited.

The sampling frequency is:$$F_s=\frac {6} {10^{-6}} HZ$$ and $$T_s=\frac {10^{-3}} {6}$$

My question: what is the maximum frequency? how we can calculate it?

I know the answer that is $$F_n=\frac {9} {10^{-6}} HZ$$

But how?

What you need to know is that the signal $$s(t)=e^{j2\pi f t}$$ has frequency $f$. At any time $t_0$, $s(t_0$) is a point on the unit circle. As time advances, this point moves on the circle, and it completes $f$ turns in one second.

Knowing this, you can easily see that the maximum frequency in your signal $x_c(t)$ is 9000 Hz. That is because, for $k=9$, you can write the exponent as $j2\pi 9\cdot 1000 t$.

The sampling frequency you mention has no bearing on the signal's maximum frequency. However, if you intend to sample this signal, you need to make sure that the sampling frequency is at least twice the maximum frequency, or $f_s=18000$.


Rayleigh energy equation - applied to Fourier Transform of a sinusoid I saw details indicating "total energy of a Fourier-transformable signal equals the total area under the curve of squared amplitude spectrum of this signal."

I have no issue with that. But I have found pretty much no discussions about the consequence of applying the Rayleigh energy equation to the Fourier transform of sinusoidal signal - which of course has infinite energy. There are certainly discussions indicating that the energy equation can only be applied to finite energy signals, and therefore is not applicable to the Fourier transform of a sinusoidal time domain signal.

So - what is the bottom line? Does this mean that if one attempts to apply the energy equation to the Fourier transform of a sinusoid, then the result will be a finite value or finite number - and it will be a meaningless number, right?

That is - when applied to finite energy signals, the energy equation returns meaningful results. But when applied to infinite energy signals like a sinusoid, then it returns a value, but it will be a meaningless result, right?

An (infinite) sinusoidal signal does not really have a Fourier transform, since as you say it is not a finite-energy signal. It can be represented as a Fourier series, and there is a useful way to sort of write it as a "pseudo-Fourier transform" with Diract delta distributions, as well. However, if you try to take the magnitude of that spectrum and square it, it is not well-defined, since you would need the squared-magnitude of the Dirac deltas.


Cramér-Rao lower bound I have been trying to implement the Cramér-Rao lower bound from the paper - A reference-free time difference of arrival source localization using a passive sensor array (eq. 6 and eq. 7).

$$
\eqalignno{ &{\rm cov}(\hat{d}_{m,n},\hat{d}_{m,n}) \,\,=\,\, {3c^2\over \pi^{2}N\kappa^{2}}{\rho_{m}\rho_{n}\over {\rm SNR}}\left(2+{\rho_{m}\rho_{n}\over{\rm SNR}}\right)&\hbox{(6)}\cr &{\rm cov}(\hat{d}_{i,j},\hat{d}_{k,\ell})\,\,=\,\, {3c^2 \over \pi^{2}N\kappa^{2}}{\rho_{i}\rho_{j}\rho_{k}\rho_{l}\over {\rm SNR}^{2}}\phi_{i,j,k,l}&\hbox{(7)}}
$$

As can be seen the bound is dependent on length of the signal.

Whats the ideal value one must take for the length of the signal?

The bound was derived assuming that $N$(length) tends to infinity (reference 8).

So are we supposed to take huge value for $N$ like $10^{21}$?

Sounds like you have read reference 8 and understand that for an infinite $N$ you can one could theoretically achieve the CRLB. But clearly in practice you need to drop that idea and think about what your constraints are for your problem. That is, do you really need an estimate with such low variance and if you do then how long are you willing to wait? Since increasing $N$ means you need more and more samples of your signal to come but if you are doing offline processing then maybe you can afford this and can take a very large $N$.


Prewarping both resonant frequency $f_0$ and bandwidth (or $Q$) when using bilinear transform RBJ's Audio EQ cookbook takes into account only frequency prewarping when  case Q is used for bandwidth. Why not $Q$ prewarping as well with some of those filter types defined there?
At KVR I got info that BLT turns to kind of II (Impulse Invariant) when both prewarpings are in use and BLT looses some nice properties in that case.

Well, I tried the $Q$ prewarping and made few plots:

Right column plots uses EE's $Q$ definition just for the (original) cookbook. Did not check the situation of phases.

If my plots are correct then it looks like prewarping $Q$ might bring some improvement at least for peak filter magnitude. For LP filter, situation is different.

BTW, some doctoral thesis showed plots where both prewarpings were used and set to have same factor. How's that factoring done?

Any thoughts?

EDIT: Plots showing the end results for this query.

Clarifications can be found from RBJ's answer and comments related to that answer.

alright, there are a couple of different issues that may (or may not) need to be de-conflated.

i'm gonna try to keep the number of symbols minimized.

$dB_\text{gain}$ is the number of dB gain of the peak (for $dB_\text{gain} > 0$) or cut (for $dB_\text{gain} < 0$).  it appears to be 6 dB in the plots.

$$A^2 \triangleq 10^{dB_\text{gain}/20}$$

is the linear gain for the boost or cut.

$f_\text{s}$ is sample rate.

$f_0$ is frequency of the boost or cut in the same units as the sample rate.

the analog transfer function of a resonant second-order filter (a.k.a. "biquad" or "SOS") is

$$\begin{align}
\\
H(s) &= \frac{b_0 + b_1 s + b_2 s^2 }{1 + a_1 s + a_2 s^2} \\
\\
 &= \frac{b_0 + b_1 s + b_2 s^2}{1 + \frac{1}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2} \\
\\
\end{align} $$

don't mistake the $a_k, \ b_k$ coefficients of the analog prototype with those in the resulting digital filter from the recipe in the cookbook. this is how the "EE definition" of $Q$ is defined.

an analog BPF with passband gain of 0 dB, has transfer function:

$$ H(s) = \frac{ \frac{1}{Q}\frac{s}{2 \pi f_0} }{1 + \frac{1}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2} $$

note when $f = f_0$ then $H(j 2 \pi f_0) = 1$.  note also that there are two bandedges $f_+$ and $f_-$ such that $f_- < f_0 < f_+$ and

$$ f_+ = f_0 \cdot 2^{bw/2} $$
$$ f_- = f_0 \cdot 2^{-bw/2} $$

and

$$ |H(j 2 \pi f_-)|^2 = |H(j 2 \pi f_+)|^2 = \frac{1}{2} \quad . $$

we define those bandedges for the BPF to be the "half-power frequencies" (a.k.a. "-3 dB frequencies").  and the $bw$ parameter is the bandwidth expressed in octaves.  the higher bandedge $f_+ = 2^{bw} f_-$ is $bw$ octaves higher log frequency than the lower bandedge $f_-$.  turns out that this bandwidth $bw$ in octaves is related to $Q$ as follows:

$$\begin{align}
\\
\frac{1}{Q} &= \frac{2^{bw} - 1}{2^{bw/2}} \\
\\
 &= 2 \ \sinh \left( \frac{\ln(2)}{2} \ bw \right) \\
\end{align} $$

keeping the same definition of $Q$, the "bell-shaped" boost/cut parametric EQ takes that BPF, gives it some gain (with sign) and adds it to a wire.  it has transfer function:

$$\begin{align}
H(s) &= (A^2 - 1)\frac{ \frac{1}{Q}\frac{s}{2 \pi f_0} }{1 + \frac{1}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2} \ + 1 \\
\\
 &= \frac{1 + \frac{A^2}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2}{1 + \frac{1}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2}\\
\end{align}$$

this is the "traditional" analog parametric EQ.  note that
$$|H(j 2 \pi f_0)| = A^2 = 10^{dB_\text{gain}/20} \quad .$$
the problem is that, leaving $f_0$ and $Q$ constant, the curve for $dB_\text{gain} > 0$ is not a mirror image for $dB_\text{gain} < 0$, given the same number of dB, $|dB_\text{gain}|$, boost or cut.  the cut will be much skinnier than the boost.

some people might want the cut to exactly undo the boost given all other parameters being the same.  so in the cookbook (and in some other papers), a redefinition of $Q$ is made.  this redefinition is the substitution:

$$ Q \ \leftarrow \ A \cdot Q $$

so that makes the transfer function for the parametric EQ (with adjusted $Q$):

$$\begin{align}
 H(s) &= \frac{1 + \frac{A^2}{A \cdot Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2}{1 + \frac{1}{A \cdot Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2} \\
\\
  &= \frac{1 + \frac{A}{Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2}{1 + \frac{1}{A \ Q}\frac{s}{2 \pi f_0} + \left(\frac{s}{2 \pi f_0} \right)^2} \\
\end{align} $$

note that if $dB_\text{gain}$ is replaced with $-dB_\text{gain}$, then $A$ is replaced with $\frac{1}{A}$ and the numerator and denominator of $H(s)$ are essentially swapped, which causes the frequency response of the cut to mirror that of the boost.

now, using the same relationship between bandwidth $bw$ and $Q$:

$$ \frac{1}{Q} = 2 \ \sinh \left( \frac{\ln(2)}{2} \ bw \right) $$

then the bandedges

$$ f_+ = f_0 \cdot 2^{+bw/2} $$
$$ f_- = f_0 \cdot 2^{-bw/2} $$

satisfy this gain definition:

$$ |H(j 2 \pi f_-)| = |H(j 2 \pi f_+)| =  A = 10^{dB_\text{gain}/40} $$

which are the "mid-gain frequencies" having dB gain of $\frac{dB_\text{gain}}{2}$.

so the definition of bandedge gain is a bit different, but at least for the analog prototype, we're keeping the relationship between $bw$ and $Q$ the same.  higher $Q$ means tighter $bw$.

so, first, before we discuss "warping Q", let's be completely consistent about which Q to compare.  for this, i might recommend leaving the "EE definition" of Q behind to not confuse.

so the bilinear transform which compensates for the frequency warping at the "significant frequency" or the resonant frequency $f_0$ makes this substitution:

$$ \text{normalized }s \triangleq \frac{s}{2 \pi f_0} \ \leftarrow \ \frac{1}{\tan(\pi f_0/f_\text{s})} \ \frac{1 - z^{-1}}{1 + z^{-1}} $$

$H(z)$ is the resulting digital filter transfer function after making that substitution for normalized $s$.

the frequency warping function turns out to be:

$$ f_\text{analog} = \frac{f_0}{\tan\left(\pi \frac{f_0}{f_\text{s}} \right)} \ \tan\left( \pi \frac{f_\text{digital}}{f_\text{s}} \right)  $$
or
$$ f_\text{digital} = \frac{f_\text{s}}{\pi} \ \arctan \left( \tan\left(\pi \frac{f_0}{f_\text{s}} \right) \frac{f_\text{analog}}{f_0} \right) $$

(note that DC gets mapped to DC and $f_0$ gets mapped to $f_0$, but when $f_\text{analog} \to \infty$, then $f_\text{digital} \to \frac{f_\text{s}}{2}$.)

fixing (or "prewarping") $f_0$ does not compensate for a cramped bandwidth because, besides the resonant frequency $f_0$, so also are the bandedges $f_+$ and $f_-$ warped by the bilinear transform.    measured in octaves the bandwidth, $BW$, in the digital filter is:

$$\begin{align}
 BW &= \log_2\left( \frac{f_\text{s}}{\pi} \arctan \left( \tan\left(\pi \frac{f_0}{f_\text{s}} \right) \right) \frac{f_+}{f_0} \right) - \log_2\left( \frac{f_\text{s}}{\pi} \arctan \left( \tan\left(\pi \frac{f_0}{f_\text{s}} \right) \frac{f_-}{f_0} \right) \right)  \\
\\
&= \log_2\left(\arctan \left(\tan\left(\pi \frac{f_0}{f_\text{s}}\right) 2^{+bw/2} \right) \right) - \log_2\left(\arctan \left(\tan\left(\pi \frac{f_0}{f_\text{s}}\right) 2^{-bw/2} \right) \right)  \\
\\
&= \log_2\left(\frac{\arctan \left(\tan\left(\pi \frac{f_0}{f_\text{s}}\right) 2^{+bw/2} \right)}{\arctan \left(\tan\left(\pi \frac{f_0}{f_\text{s}}\right) 2^{-bw/2} \right)} \right)  \\
\\
 \end{align}$$

you can see that the mapping of $bw$ to $BW$ is an odd-symmetry function, so it goes through $0$ and has no even-order terms in a Maclaurin (a.k.a. Taylor series) expansion.  now if you fix $f_0$ and plot digital $BW$ vs. analog $bw$, you will see that $BW < bw$ and that is the bandwidth cramping done by frequency warping of the bilinear transform.  to uncramp the bandwidth, you would have to solve for $bw$ in terms of $BW$ and $f_0$ and $f_\text{s}$.  and that is (how shall we say?) a female canine.

if you read the cookbook a little, you will notice a slightly adjusted mapping between $BW$ and $Q$:

$$ \frac{1}{Q} = 2 \ \sinh \left( \frac{\ln(2)}{2} \ BW \frac{2 \pi f_0/f_\text{s}}{\sin(2 \pi f_0/f_\text{s})} \right) $$

so the digital filter bandwidth $BW$ was increased by a factor of $\frac{2 \pi f_0/f_\text{s}}{\sin(2 \pi f_0/f_\text{s})}$ as a first-order attempt to compensate for the bandwidth cramping.  this result comes from assuming a narrow bandwidth in the first place and evaluating the derivative:

$$ \frac{\partial \, BW}{\partial \, bw} \Bigg|_{bw = 0} = \frac{\sin(2 \pi f_0/f_\text{s})}{2 \pi f_0/f_\text{s}} $$

the first term of the Maclaurin series is

$$ BW \approx \frac{\sin(2 \pi f_0/f_\text{s})}{2 \pi f_0/f_\text{s}} bw $$

compensating (or "prewarping") that first-order cramping of the bandwidth (which is what i think you mean by "prewarping Q") is done in the cookbook and always has been.  but it is only a first-order compensation.  if you want to do it better, be my guest and try to compute the next third-order term for the Maclaurin expansion (and flip it around so you have an expression of the approximate analog $bw$ in terms of digital $BW$).  doing this with Maclaurin series expansion is the only way i know because i don't think you will invert the $bw \to BW$ mapping above directly.


How to estimate the poles of a linear combination of exponentials? Consider a signal that can be expressed as a linear combinations of exponentials:

$x(n) = \sum_{k=1}^K c_k z_k^n$

I would like to estimate the poles $z_k$ and I was told that I could estimate them by calculating the SVD of the Hankel matrix of the observations $\mathbf{H}$:

$\mathbf{H} = \left(\begin{array}{cccc}x(1) & x(2) & x(3) & \dots \\
x(2) & x(3) & x(4) & \dots \\
x(3) & x(4) & x(5) & \dots \\
\vdots & \vdots & \vdots &  \\\end{array}\right)$

But how exactly do I utilize the SVD? The poles will be the singular values? And is it possible to find $K$?

Thanks in advance!

My transform, EST, finds exactly what you are looking for.

The steps of EST are as follows.

Let N be the number of exponentials.

Find the LPC (Linear Prediction Coefficients) for the signal.

For the sake of simplicity I will assume that all the exponential bases (the z values) are complex. The size of the LPC has to be 2N.

You must have at least 4N samples.
Case 1 is when you have exactly 4N samples. Create a system of linear equations with the Hankel matrix and the last 2N samples as the results vector. The LPC vector is the solution to this system.

Case 2 is when you have more than 4N samples. In this case you find the LPC using the covariance method.

Convert the LPC vector to a prediction polynomial, exactly as LSP (Line Spectral Pairs) is doing.
Find the roots of the polynomial.

The roots of the polynomial are the exponential bases you are looking for.

If the signal is noise-free, then the exponential bases are precise. If not, it is an estimate.

Once the exponential bases are known, finding the coefficients (the c values), is simply a matter of linear equations.

You also asked if you can find the number of exponentials.
In the case of a noise-free signal: yes.
Otherwise, you should find a large enough number of exponentials and filter out some of them.
Please note that each exponential base holds a frequency and the growth/decay rate of the signal at that frequency and each coefficient holds an amplitude and a phase. These parameters can be used for filtering criteria.


Onset Detection Peak Picking I have a guitar music signal and I want to detect the onset points via energy.

This is the signal example:

This is the energy plot:

My question is what would be the best method defining the threshold and consequently finding the peaks. I reckon I will have to set a window length and a window step (I am thinking 1000 and 100 samples respectively on a 44.1KHz sampling rate) and then every loop defining a new threshold (adaptive/local threshold). Every loop I would have to find the window average and compare it to the threshold. I have these questions:

1) Does window length = 1000 samples and window step = 100 seems reasonable?

2) How the threshold would be defined each loop and in relation to what?

3) If a window is above length which samples in it would I consider peaks? I think it would be somewhere in the middle but what range around it?

Maybe, window middle - window step < peak samples < window middle + window step  ?

Thank you for your help.

compute RMS from audio signal to get power
do AGC (automatic gain control)
perform "discrete differentiation" (the simplest is 1st order: $y[i] = x[i] - x[i-1]$)
if the value is greater than certain threshold, it means we have an onset. You have to determinate the threshold experimentally or use adaptive algorithm. Obviously you also need some kind of lowpass filter to filter the noise.

You may want to swap steps 2 and 3.

If you can use software libraries, check out aubio.

Possible duplicate of https://stackoverflow.com/questions/294468/note-onset-detection


RF Downconversion using Heterodyne Principle after a RF down-conversion using heterodyning principle, will there be an information loss?
If we translate to a fixed IF frequency, is there is an intuitive way of explain that the information can be retrieved?
Is it correct to say that it is the relative change that matters during demodulation?
For example for AM it is the amplitude change and for FM its
frequency change and this does not change even if we down-convert to
a fixed IF frequency?

after a RF down-conversion using heterodyning principle, will there be an information loss?

No, if the original signal was band-limited, and the bandwidth of your IF processing is sufficiently large to capture that. (so, if you've built a sensible heterodyne receiver.)

  If we translate to a fixed IF frequency, is there is an intuitive way of explain that the information can be retrieved?

Yes. Shifting a frequency is just multiplications with $e^{j\Delta_f t}$, and that's an easily invertible operation (multiply with $e^{-j\Delta_f t}$).

If an operation is perfectly invertible, no information was lost.

  Is it correct to say that it is the relative change that matters during demodulation?

I don't understand that sentence, but I think you really mean "what matters is the bandpass content, not the carrier frequency", and that would be right. See: equivalent complex baseband.

  For example for AM it is the amplitude change and for FM its frequency change and this does not change even if we down-convert to a fixed IF frequency?

This depends on how you do the AM and FM demodulation – if they work the same no matter the carrier frequency, then clearly yes.


Properties of root raised cosine filter In communication systems, the transmitted signal is often passed through a root-raised cosine filter to avoid, in textbook language, intersymbol interference (ISI). An identical filter at the receiver end processes the signal so that the combined response of the two filters is a raised cosine filter.

The raised cosine filter has two properties: roll-off factor and oversampling rate.

I have a bipolar sequence ($\pm 1$ values) that I am supposed to process through a raised cosine filter. It is not clear to me what exactly is intersymbol interference here, and how can I visualize it. Here the symbol has only one bit. How should I choose the roll-off factor and the oversampling rate?

-ryan

Your bipolar sequence at one sample per symbol is equivalent to a perfectly synchronized sequence of symbols with no intersymbol interference (ISI). Intersymbol interference in this context refers to the tendency for symbols that precede and follow the current symbol to interfere at the receiver. That is, at the correct sampling instant for symbol $N$, the receiver also sees some of symbols $N-1$ and $N+1$ (and possibly more in each direction) mixed in. This causes an increase in symbol error rate.

In general, ISI is caused by some non-ideal channel impulse response between your transmitter and receiver. This can be caused by several phenomena, including:

A multipath communications channel (e.g. a wireless channel). In this case, the receiver observes a signal that can be modeled as your originally transmitted signal convolved with some impulse response.
The use of a non-ideal pulse shape. On the AWGN channel, the optimum receiver consists of a filter whose impulse response is matched to the transmitted pulse shape. The receiver therefore observes a stream of pulses where each has the shape of the originally transmitted pulse shape convolved with itself. The choice of the pulse shape is important to ensure that the receiver doesn't encounter any ISI after applying the matched filter.

This is why the root-raised cosine filter is often used. When convolved with itself, the result is a raised-cosine response, which has zero ISI (it is exactly zero for all multiples of the symbol period $T$, apart from $t=0$).

With that said, there are some practical considerations when designing your pulse shape. When you actually transmit it, you'll need to oversample it in order to meet the sampling theorem's reconstruction requirements. A typical oversampling rate might be 4.

The rolloff factor controls how much bandwidth the RRC-shaped signal will occupy. In general, the more excess bandwidth, the better you can expect your synchronization algorithm to perform (i.e. you would usually observe smaller error from your timing recovery subsystem with a larger rolloff factor). However, this comes at a cost of more bandwidth required to carry the signal.


What is the benefit of FFT when we're already sampling at an acceptable rate? As I'm learning about DSP I can understand the process of taking a continuous signal into a discrete signal. I understand an oscillator producing different types of signals, including a periodic signals (especially sinusoids), phase shifts and amplitudes etc.,

I'm following a coursera class, and reading DSP theory.. Rather than blindly using the DSP libraries out there I want to know exactly why something is happening. My 1st roadblock is this:

If we already have enough samples of amplitudes collected at a reasonably sample rate (using Nyquist-Shannon) Why does it matter that a signal (Sinusoid) can be broken down into other sinusoids? Aren't the samples enough to reconstruct the original signal/sound?
What is advantageous by knowing that the Fast Fourier Transform can decompose a signal into multiple signals?
The sinusoids that make up a single sinusoid also have other sinusoids that compose them. If a signal has say 5 sinusoids, those 5 sinusoids, according to DFT if I'm correct, can also each have sinusoids.. Isn't this recursive or infinite? At what point can we say we have "enough" sinusoids to represent a signal.

Thanks in advance for the clarification & patience.

Not all periodic signals are sinusoids.  Items 1 and 3 should not assume such.  Not all sampling is done just for reconstruction.  Look up other needs or uses for sampled signals.
A linear time-invariant system can behave differently to different frequency sinusoidal inputs (as provided by common solutions to certain low order differential equations), so it may be easier to figure how that system responds to other types of far more complicated inputs by breaking them down into the sum of simpler sinusoidal ones.
A DFT or FFT won't break down a pure (unmodulated) sinusoid that is integer periodic in aperture into other sinusoids (other than that a pure sinusoid of an arbitrary phase can be broken down into component cosine and sine sinusoids of that same frequency).  But an FFT can "break down" other types of signals.


Energy of audio signal From audio file [y,FS]=something.wav

So the y is typically N cross 2 matrix

I want calculate energy so i.e  sum(abs(y.^2))

So when I apply these I end with two answer for each column. So I want only one energy for whole signal what to do???

That really depends on what you need to use the energy for and how exactly you define "energy" in your application.

You can either take the sum of the energies or do a more advanced mode taking into account the inter-channel correlation and the spectral correlation of stereo playback. If you want anything like "physical" energy, you need a calibration.


Biquad and polarity flipping Assuming normalized biquad coefficients:

$$
H(z) = \frac{b_0 + b_1z^{-1} + b_2z^{-2}} {1 + a_1z^{-1} + a_2z^{-2}}
$$

And taking the definition of flipping the polarity as multiplying every tap in a signal by -1 i.e.

$$[0.2, 0.1, -0.1, -0.2] \to [-0.2, -0.1, 0.1, 0.2]$$

For a zero delay biquad is it sufficient to check that $b_0 >= 0$ to assume the biquad is not flipping the polarity accidentally (although it will be modifying the waveform) and if it were, a correction $H'(z)$ would be:

$$
H'(z) = \frac{-b_0 - b_1z^{-1} - b_2z^{-2}} {1 + a_1z^{-1} + a_2z^{-2}}
$$

This relates to audio.

A filter will not flip the polarity for all input signals. It might flip it for certain input frequencies, namely for the ones for which the filter's phase response equals $\pi+2k\pi$. The way to check if your filter flips the polarity of a relevant input signal (i.e., for a certain frequency), you have to evaluate the frequency response at that frequency. For a low pass filter you can check the filter's behavior at DC (i.e., frequency zero):

$$H(1)=\frac{b_0+b_1+b_2}{1+a_1+a_2}\tag{1}$$

You want this value to be positive. For a high pass filter you may want to check the frequency response at Nyquist:

$$H(-1)=\frac{b_0-b_1+b_2}{1-a_1+a_2}\tag{2}$$

And for a band pass filter you should check the frequency response at the center frequency.

In sum, a non-trivial filter will not flip the sign for all input frequencies, it might flip it for specific frequencies. Depending on the filter type, there are certain prominent frequencies (e.g., DC or Nyquist) that you need to check.

EDIT:

It is straightforward to show that checking the sign of $b_0$ is not sufficient. Take the coefficients

b = [0.13111  -0.26221  -0.13111];
a = [1.00000  -0.74779   0.27221];

We have $b_0>0$, but the frequency response at DC is

$$H(1)=-0.5$$

and so the sign of a constant input signal will be flipped (as soon as the transients have died out):


Building periodic sequences from finite-support sequences Given a discrete-time finite-support signal x[n]

$$x[n] = \left\{ {\begin{array}{*{20}{l}}
{{{( - 1)}^n}n}&{{\rm{ }}n = 1,2,3}\\
0&{{\rm{otherwise}}}
\end{array}} \right.$$

And consider also its periodic repetition is

$$y[n] = \sum\nolimits_{k =  - \infty }^\infty  {x[n + 7k]} $$

How we can generate the y[n] in Matlab or Python?

Here is piece of code to generate and plot x[n]. But for Y[n] I have no idea!

% finite support signal
n = -5:1:8;
x = zeros(size(n));
x( find(n==1 | n==2 | n==3) ) = [(-1)^1*1 (-1)^2*2 (-1)^3*3];
stem(n, x, 'filled')


Derivation of $ R_{N(t)}(\tau) $ from its $f_{N(t)}(\eta)$ How can we prove the auto-correlation function of white gaussian noise $\{ R_{N(t)}(\tau) \}$ is $\frac{N_0}{2} \delta(\tau)$ from its p.d.f in equation below?
$$
f_{N(t)}(\eta)=\frac{1}{\sqrt{2 \pi \sigma^2}}\exp\left(\frac{-(\eta-\mu)^2}{2 \sigma^2}\right)
$$

My Approach:

We know: $\frac{N_0}{2}=\sigma^2 \quad(\text{ or variance}) \tag{1}$
We also know:
$$
R_{X(t)}(\tau)= E[X(t)X(t+ \tau)]=\int_{y=- \infty}^{\infty} \int_{x= - \infty}^{\infty} x \space y \space f_{X(t) X(t+\tau)}(x,y)  \space dx  \space dy
$$
Where $f_{X(t) X(t+\tau)}(x,y) \text{ is joint p.d.f}$
$$\text{so , } R_{N(t)}(\tau)= E[N(t)N(t+ \tau)]=\int_{y=- \infty}^{\infty} \int_{x= - \infty}^{\infty} x \space y \space f_{N(t)N(t+\tau)}(x,y)  \space dx  \space dy$$
Now we know
$$f_{N(t)}(\eta)=\frac{1}{\sqrt{2 \pi \sigma^2}}\exp\left(\frac{-(\eta-\mu)^2}{2 \sigma^2}\right)$$
then, how can we find $f_{N(t)N(t+\tau)}(x,y)$ ?

As @MattL points out in a comment, a Gaussian pdf does not imply whiteness. Indeed, it can be argued that the assumption that the process is a continuous-time white noise process is contrary to the belief that the random variables constituting the process even have a pdf of any kind.  Continuous-time white noise is a mythical beast that is used to account for the empirical observation that the power spectral density $S_{\scriptstyle{\text{output}}}(f)$ of the noise at the output of a linear filter with transfer function $H(f)$ can be expressed as
$$S_{\scriptstyle{\text{output}}}(f) = \frac{N_0}{2}|H(f)|^2,-\infty < f < \infty,$$
a result that is perfectly consistent with the standard theory of wide-sense-stationary processes in linear systems which claims that
$$S_{\scriptstyle{\text{output}}(f)} = S_{\scriptstyle{\text{input}}(f)}|H(f)|^2, -\infty < f < \infty,$$
provided that we are willing to suspend belief and assume that the input noise has power spectral density $\frac{N_0}{2}, -\infty < f < \infty,$ which has infinite power. (This suspension of belief must be accompanied by a most solemn vow to not attempt to solve the world's energy crisis by plugging in the power grid into a white noise generator!). Formally, a continuous-time white noise process is defined to be a continuous-time zero-mean wide-sense stationary process whose power spectral density has constant value for all choices of frequency $f$. Continuous-time white noise can only be observed through some kind of device (whether linear or nonlinear) which perforce performs some kind of filtering or bandwidth reduction so that the observation has finite variance and the observed process has finite power. The observed process is not called a white noise process.

The inverse Fourier transform of $\frac{N_0}{2}, -\infty < f < \infty,$ is the autocorrelation function $\frac{N_0}{2}\delta(t)$ of the white noise process where $\delta(t)$ is the Dirac delta. It is a shibboleth often murmured on this site that the value $\delta(0)$ of the Dirac delta is $\infty$ with the usual consequences for this mispronunciation being administered by the cognoscenti of dsp.SE, but the point is that we cannot define any kind of pdf for the random variables constituting a white noise process, and very certainly not the Gaussian pdf that the OP wishes to use.  What is true that the noise at the output of the linear filter is Gaussian noise meaning that the random variables constituting the output noise process are (zero-mean) Gaussian random variables with variance
$$\sigma^2 = \int_{-\infty}^\infty S_{\scriptstyle{\text{output}}(f)}\, \mathrm df = \frac{N_0}{2}\int_{-\infty}^\infty |H(f)|^2\mathrm df.$$
This result is also consistent with the view that the input noise process is white Gaussian noise in the sense that it is a standard result in the theory of wide-sense-stationary process that linear filtering of a Gaussian process results in a Gaussian process. But it is not possible to ascribe a variance  (not even an infinite variance) to these alleged Gaussian random variables in the input process because for all $x$,
$$\lim_{\sigma\to\infty} \frac{\exp\left(-\frac{x^2}{2\sigma^2}\right)}{{\sqrt{2\pi}\sigma}} = 0.$$

  In summary, continuous-time white Gaussian noise is a myth used to provide a rationale for the empirical observation that the noise at the output of a linear filter with transfer function $H(f)$ is a zero-mean wide-sense-stationary (and hence strictly stationary) Gaussian random process with power spectral density
  $$S_{\scriptstyle{\text{output}}}(f) = \frac{N_0}{2}|H(f)|^2,-\infty < f < \infty.$$
  The only assumption about the random variables constituting the wide-sense-stationary input white noise process is that they all are zero-mean random Gaussian random variables but we cannot ascribe any variance to them or observe them directly or do anything else with them until we have processed them through some device of finite bandwidth.

More than what you probably want to know about continuous-time white Gaussian noise
can be found in the Appendix of this lecture
note
of mine.

In view of a veritable flood of comments from @markleeds on the original version of this answer, it is worth noting that on dsp.SE (and any place else where MATLAB is the lingua franca), there is also the concept of a discrete-time white noise process.  A discrete-time process is a countable collection of random variables $\{X[k]\colon k \in \mathbb Z\}$ (notice the change from $X(t)$ to $X[k]$: this is very commonly-used notational convention to distinguish the two cases) and a discrete-time white noise process is just a collection of zero-mean finite-variance independent identically distributed random variables $N[k]$, $k \in \mathbb Z$. The (discrete-time) autocorrelation function of the process is $$R_N[k] = E[N[m]N[m+k]] = \sigma^2\delta[k]$$ where $\delta[k]$ is the Kronecker delta (not the Dirac delta) which is defined as
$$\delta[k] = \begin{cases}1, & k=0,\\0, & k\neq 0.\end{cases}$$
Although the same symbol is used in both cases, we use $()$ and $[]$ as clues to tell us what $\delta$ means. Now, what, if any, is the relationship between samples of continuous-time white noise and discrete-time white noise? Well, as mentioned above, sampling white noise is not a physical reality -- one cannot draw blood from a beast that does not exist! -- and any attempt to do so must take into account the physical properties of the sampler which is a finite-bandwidth device with a switch that closes for a very short period of time $\varepsilon$ and then opens again. A simple model for what happens during this time is that what the sampler reads at time $t$ is the value of
$$\int_{t-\varepsilon}^{t}N(\tau)\, \mathrm d\tau$$ which
is a random variable  with mean $0$ and variance
\begin{align}
\sigma^2 &= E\left[\int_{t-\varepsilon}^{t}N(\tau)\, \mathrm d\tau\int_{t-\varepsilon}^{t}N(\lambda)\, \mathrm d\lambda\right]\\
&= E\left[\int_{t-\varepsilon}^{t}\int_{t-\varepsilon}^{t}N(\tau)N(\lambda)\, \mathrm d\tau \,\mathrm d\lambda\right]\\
&= \int_{t-\varepsilon}^{t}\int_{t-\varepsilon}^{t}E[N(\tau)N(\lambda)]\, \mathrm d\tau \,\mathrm d\lambda\\
&= \int_{t-\varepsilon}^{t}\int_{t-\varepsilon}^{t}\frac{N_0}{2}\delta(\tau-\lambda)]\, \mathrm d\tau \,\mathrm d\lambda\\
&= \int_{t-\varepsilon}^{t}\frac{N_0}{2}\,\mathrm d\lambda\\
&= \frac{N_0}{2}\varepsilon
\end{align}
If $N(t)$ is continuous-time white Gaussian noise with power spectral density $\frac{N_0}{2}$ and the samples are spaced at least $\varepsilon$ apart, say the $k$-th sample is at time $kT$ where $T\geq \varepsilon$, then we have that the $N[k]$ are independent zero-mean random variables with variance $\frac{N_0}{2}\varepsilon$, that is, _discrete-time white noise.  If $N(t)$ is a Gaussian process, then the discrete-time process is a Gaussian process in which the random variables are independent $\mathcal N(0,\frac{N_0}{2}\varepsilon$ random variables.


Combining two FIR band-pass filters I am using IPP for FIR filtering over 1D data. I want to use the same input and filter it using two band pass filters. This is how I do it:

output1 = doFIRBandPassFilter(input,low1,high1);
output2 = doFIRBandPassFilter(input,low2,high2);

I want to combine those outputs so that the end result will show BOTH (low1,high1) and (low2,high2) filtered data.

If output1 and output2 are arrays - should I add them? multiply them? convolve them?

As you've indicated there may be some mismatching in the filter bands resulting in undesirable results. If one of the filters is heavily attenuating the frequencies of the other filter, there's little you can do to reconcile this by just cascading the filters together.

Given you'd like to maintain the response of both filters independent of the other, you need to add the filter outputs together.

Here's some code that shows the result of one filter heavily attenuating the frequencies of the other. It also shows the equivalence of cascading two filters either by connecting the output to the input of the next stage or just convolving the filter taps.

b1 = fir2(2048, [0 0.1 0.15 0.3 0.35 1], [0 0 1 1 0 0]);
b2 = fir2(128, [0 0.4 0.45 0.6 0.65 1], [0 0 1 1 0 0]);

x = randn(1, 32768);

% Parallel sum
y1 = filter(b1, 1, x) + filter(b2, 1, x);

% Cascade
y2 = filter(b2, 1, filter(b1, 1, x));

% Equivalent cascade
y3 = filter(conv(b1, b2), 1, x);

figure, hold all
subplot(311), plot(20*log10(abs(fft(y1))))
subplot(312), plot(20*log10(abs(fft(y2))))
subplot(313), plot(20*log10(abs(fft(y3))))


Adding AWGN noise with certain power Assume I have a signal(BPSK) and it has a certain amplitude. I would like to add a noise with a certain power to it. Assume I am using coding rate Rc as well. One method I am thinking would be to consider

EbNo="value" Eb=sum(signal.^2)/(length(signal)

sigma^2=Eb/(2*Rc*EbNo);
noise=sqrt(sigma^2)*randn(1,length(signal));
rec=signal=noise;

However it seems power of my noise does not match what I expect. Does anyone have any other suggestion or algorithm or do you guys know if there is anything wrong with what I wrote.

This is solved slightly differently if you're simulating a continuous-time or discrete-time AWGN channel. I'll assume you're simulating a discrete-time channel; that is, you're simulating the output of the receiver's matched filter sampled at correct time instants.

First consider the case with no noise. In this scenario, the matched filter's output is an array of numbers, each of which is equal to $+\sqrt{E_b}$ or $-\sqrt{E_b}$. The average energy per sample is then equal to $E_b$.

Now consider the case where only white noise is present at the matched filter's input. The sampled filter output is a sequence of random values with variance $\sigma_n^2$. The signal-to-noise ratio is defined as $E_b/\sigma_n^2$. By convention, $\sigma_n^2$ is denoted $N_0/2$, which means that $\text{SNR}=2E_b/N_0$. (This is because when the filter input is white noise with power spectral density $N_0/2$ and the matched filter's impuse response has energy equal to 1, then the variance of the sampled filtered noise is also $N_0/2$. This is a technical detail that you may omit for now and come back to when you're more comfortable with the theory).

Then, $$\sigma_n^2=\frac{N_0}{2}=\frac{E_b}{SNR}.$$ Using this equation you can easily simulate any SNR you want. A word of caution: it is common, when plotting the bit error rate vs the SNR, to use $E_b/N_0$. Be aware of this when comparing your results to those published elsewhere.


Unable to obtain correct correct fourier transform for homework Okay so I have this question for my homework which I am asked to find the fourier transform of a signal. However I am unable to obtain the correct answer, the question and my working is shown below.

Question( The answer is the equation of Xp(f) ):

My working(Look from b) onwards):

As can be seen, I am unable to get the answer. But I don't feel that I am wrong.

Thanks.

If $x(t)$ is the (non-periodic) function in the interval $[-2,2]$, then the periodic function $x_p(t)$ is given by

$$x_p(t)=x(t)\star \sum_{n=-\infty}^{\infty}\delta(t-6n)\tag{1}$$

where $\star$ denotes convolution, and, consequently, its Fourier transform is

$$X_p(f)=X(f)\cdot\frac16\sum_{k=-\infty}^{\infty}\delta\left(f-\frac{k}{6}\right)=
\frac16\sum_{k=-\infty}^{\infty}X\left(\frac{k}{6}\right)\delta\left(f-\frac{k}{6}\right)\tag{2}$$

where $X(f)$ is the Fourier transform of $x(t)$. It remains to find $X(f)$. The function $x(t)$ can be written as

$$x(t)=\sin\left(\frac{\pi t}{2}\right)\cdot\left[\text{rect}\left(\frac{t-1}{2}\right)-\text{rect}\left(\frac{t+1}{2}\right)\right]\tag{3}$$

With the Fourier transform pairs

$$\begin{align}\sin\left(\frac{\pi t}{2}\right)&\Longleftrightarrow \frac{1}{2j}\left[\delta\left(f-\frac14\right)-\delta\left(f+\frac14\right)\right]\\\text{rect}\left(\frac{t\pm 1}{2}\right)&\Longleftrightarrow 2e^{\pm j2\pi f}\text{sinc}(2f)\end{align}\tag{4}$$

we obtain from $(3)$

$$\begin{align}X(f)&=\frac{1}{2j}\left[\delta\left(f-\frac14\right)-\delta\left(f+\frac14\right)\right]\star 2\;\text{sinc}(2f)(e^{-j2\pi f}-e^{j2\pi f})\\&=2\left[\delta\left(f+\frac14\right)-\delta\left(f-\frac14\right)\right]\star\text{sinc}(2f)\sin(2\pi f)\\
&=2\left[\text{sinc}\left(2f+\frac12\right)\sin\left(2\pi f+\frac{\pi}{2}\right)-\text{sinc}\left(2f-\frac12\right)\sin\left(2\pi f-\frac{\pi}{2}\right)\right]\end{align}\tag{5}$$

The final expression for $X_p(f)$ results from combining $(5)$ and $(2)$.


Selected bin oscilation in STFT I am looking for a noisy signal using STFT. My window length is 128 and I am using 75% overlap. I am using a Hanning window before running the FFT process. I am using Quinn's 2nd estimator for interpolation after the fact.
What I observe is the detected signal jumping frequently between bins. In the picture below the red is the detected frequency (around 190Hz) the green is the detected bin (sampling rate is 1000Hz so bin size is 7.8Hz)

What I would like to know is - is this likely to be an artifact of the method I am using, or is it simply that the frequency is varying a lot. What would be a suitable way of coming up with a more reliable prediction without losing time resolution?
I am using a low pass filter at 20Hz to smooth the output, but it still does not remove this oscillation.

Here is a full FFT of the signal over a much long time period:

I guess that indicates that its simply the signal causing this.

I switched to using the normalized power spectrum and changed my output filtering from two-pass to one-pass. This appears to have made a massive difference:

It's difficult to tell from the picture but the response is actually much smoother and you can see that the selected bin is not jumping around in the same way. The two pass filtering I was using in particular caused a lot of ringing after sudden impulses and a one-pass filter was generally much, much better.


Is there a recommended amount of attenuation for rectifying edge artefacts? I just noticed that Spotify fades in/out the music when you play/pause - is this to avoid clicking?
This attenuation is noticeable; is there a minimum amount of attenuation required to fix clicking (while not introducing a noticeable fade). Also, is there a type of attenuation that's best, linear, parabolic, etc?
On youtube, there's no noticeable attenuation, yet there's no clicking; could this be very small attenuation?

My guess is that for Spotify, it is not just about avoiding clicks, but also to improve the DJing kind of experience. The fade length seems to be a user setting in Spotify. Based on my testing with 20 Hz and 440 Hz test tones and Virtual Audio Cable, YouTube doesn't fade in or out or do any other kind of click removal when pausing, here with 20 Hz test tone:

The worst-case test signal to be faded in or out would be a full-scale constant-valued signal (for 16-bit audio that would be 32767, 32767, 32767, 32767, 32767, 32767, ...), because it has no audible sound of its own, and because it will still allow to hear the fade as loud as possible. Then we'd need the best possible fade shape, but I don't know what it is, so here is a sub-optimal fade in:
$$f(x) = \begin{cases}0&\text{if } x < 0,\\
\frac{1}{2}-\frac{1}{2}\cos(\pi x)&\text{if } 0 \le x \le 1,\\
1&\text{if } x > 1.\end{cases}\tag{1}$$

The following "floatbeat" program plays this fade twice a second (once in and once out per second):
f = function(x){return 0.5-0.5*cos(x*PI)}, fadeLength=0.03, loopLength=1, sampleRate=44100, p=(t/sampleRate)%loopLength, q=p-loopLength/2, (p < fadeLength? f(p/fadeLength): (q < 0)? 1: (q < fadeLength? 1-f(q/fadeLength): 0))

Try running it and adjust fadeLength. At least for me, using good headphones at normal music listening volume, a fade length 0.03 s is only barely audible as a thump, and increasing it to 0.04 s makes the thump disappear. I think 0.03 s is pretty safe.
But is 0.03 s already so long that it sounds like a fade? Let's apply the fades to a 400 Hz sine tone, in floatbeat:
f = function(x){return 0.5-0.5*cos(x*PI)}, fadeLength=0.03, loopLength=1, freq=400, sampleRate=44100, p=(t/sampleRate)%loopLength, q=p-loopLength/2, (p < fadeLength? f(p/fadeLength): (q < 0)? 1: (q < fadeLength? 1-f(q/fadeLength): 0))*sin(2*PI*t*freq/sampleRate)

Running this, I can't hear a click in fade in or fade out. I know there is a fade from the way it sounds and because I do a lot of music editing, but it does not register as a gradual process, more like a smooth event.
Let's try some other fade shapes of length 0.03 s too, with my comments comparing each to Eq. 1:

Repeating Eq. 1 here for easy clicking.
A linear fade $f(x) = x$ sounds much clickier to me.
Smoothstep $f(x) = 3x^2-2x^3$ sounds maybe tiny bit worse to me.
Asymmetrical polynpomial $f(x) = 6x^2 - 8x^3  + 3x^4$ sounds significantly thumpier.
Smoothstep $f(x) = 10x^3-15x^4+6x^5$ sounds about as good, maybe even a bit better. Less mid-freq but more bassy thump.

Calculating the Fourier transform $F(f)$, where $f$ denotes frequency, of each fade of length 0.03 s appears to confirm what I heard: (For this, I calculated the Fourier transform of the derivative of each fade as function of time $t = \frac{x}{0.03s}$, and compensated for derivation by multiplying the result by a factor $-\frac{i}{f}$.)

I also included the pure tone hearing threshold curve from ISO 226:2003. Some caveats: I do not know the listening volume level so I just shifted the curve in dB scale until it matched my subjective hearing experience with the fade curves peaking above the threshold at roughly those frequencies at which I heard them. Also, the threshold curve has been measured for pure sinusoidal tones whereas here we are dealing with transient sounds. The level comparison between the two cannot be well defined. It would be better to do the analysis with something like a compressive gamma-chirp filter bank which better models hearing. I used closed-back headphones, but still my listening environment has some fan noise and such, so not perfect silence, so there could be some auditory masking going on. Despite all the caveats, I think the threshold curve is useful to show here as it seems to match what I heard. It also hints to that there may be a better fade shape that moves some of the energy from the first spectral side lobe to the main lobe, using window function terminology. Perhaps a better fade shape, one that barely touches the hearing threshold at multiple early lobes, could be obtained by integrating a suitably parameterized Ultraspherical window. Or perhaps there are ways to trace the hearing threshold for all lobes.
I'd choose between 5th order smoothstep and Eq. 1, or 3rd order smoothstep if you want to save some processing power. The choice is not too critical.


How to find the DTFT of unit step with $(-1)^n$ scalar multiplied? What are the DTFTs of the following two signals?

$$x[n] = e^{j \pi n}\left\{u[n] - u[n-8]\right\}\quad\text{and}\quad h[n] = (-1)^n\left\{u[n] - u[n-4]\right\}$$

I am trying to find $X(\omega)$ and $H(\omega)$ and the scalar $(-1)^{n}$ is throwing me off.

First note that $e^{j\pi n}=(-1)^n$. Furthermore, note that $u[n]-u[n-N]$ is a rectangular impulse of width $N$.

Now, you can calculate
$$
\textrm{DTFT}\left\{u[n]-u[n-N])\right\} = \sum_{n=0}^{N-1}e^{-j2\pi nf}=\frac{1-e^{-j2\pi Nf}}{1-e^{-j2\pi f}}
$$
which is the rule for the geometric sum.
Furthermore, you will find

$$
\textrm{DTFT}\left\{\exp{(j\pi n)}\right\}=\delta\left(f-\frac{1}{2}\right)
$$

Finally, the result follows from the convolution theorem, i.e. multiplication in time becomes convolution in frequency:

$$
\textrm{DTFT}\left\{\exp({j\pi n})(u[n]-u[n-N])\right\}=\delta\left(f-\frac{1}{2}\right)\star\frac{1-e^{-j2\pi Nf}}{1-e^{-j2\pi f}}=\frac{1-e^{-j2\pi N(f-\frac{1}{2})}}{1-e^{-j2\pi (f-\frac{1}{2})}}
$$


Determining if the given system is LTI Problem

  Given the compound system below, with the input $x(t)=\operatorname{sinc(t)}$, the output of A is $y(t)=\operatorname{sinc(2t)}$ and the output of B is $z(t)=\operatorname{sinc(t)}$, determine which of the A, B systems are Linear time invariant (LTI).

  a) A
  b) B
  c) Both
  d) Neither

My approach:

Using two facts below (if I got them right):

A and B are inverse of each other (because the overall output is
equal to the input)
If an LTI system is reversible, then the reverse is also LTI.

after a quick scanning I crossed a, b as the answered. I concluded that either both systems are LTI or neither of them are.

For the next step, using the given information I derived the system A:

\begin{align}
x(t)&=\frac {\sin(\pi t)}{\pi t}\\
y(t)&=\frac {\sin(2\pi t)}{\pi t}=\frac{2\sin(\pi t)\cos(\pi t)}{\pi t}=\frac{2\sin(\pi t)\sin(\pi t+\frac{\pi}{2})}{\pi t}=2\pi t x(t)x\left(t+\frac{\pi}{2}\right)
\end{align}

now with the $t$ as the coefficient of the system we can say $y$ is not time invariant and therfore not LTI. so I chose d as the correct answer.
But the text book I'm reading proves that b is the answer, using another method.

Questions:

So my questions are:

Is my approach/facts flawed?
If so, what is the correct approach to solve this problem.
OR possibly did the book get it wrong?

Any help would be highly apprecieted!

This is how I would approach the problem.

System $A$ cannot be LTI because its output has frequencies not present in its input. The input spectrum is a brickwall from $f=-1/2$ to $f=1/2$, but the output spectrum is a brickwall from $f=-1$ to $f=1$.

It's impossible to arrive at a conclusion for system B. It could be an ideal LTI low-pass filter, but it could also be something else (for example, it could be a signal generator).

Note that we cannot say that A and B are inverses of each other. That statement requires $x(t)=B[A[x(t)]]$ to be true for all $x(t)$. Here, we only know that is true for one specific input.


Where 2-D DFT is necessary in image processing applications instead of 2-D DCT? As we know the 2-D DCT is a real-valued kernel and less computational complexity and it used in several image processing applications like image compression etc. 2-D DFT is a complex kernel and high computational complexity, So in which image processing application of 2-D DFT is more advantageous than 2-D DCT?
Why DFT is instead of DCT in image processing applications?
Any suggestions would be greatly appreciated!
Thanks
Sukumar.

Any place you would like to filter an image with anti symmetric filter.


What is the essential bandwidth with the given function? Let the function be

$$
f(t) = e^{-3t} u(t -1)
$$
The criterion is set at 70% of the total signal energy.

I know how to find total energy, but not how to find the the $W$ (or frequency) at which the signal is a certain amount.

Find the Fourier Transform $F(\omega)$ and perform the following integration.$\frac{1}{2\pi}\int_{-W}^{W}|F(\omega)|^2\mathrm{d}\omega=0.7E$, where $E$ is the total Energy of the function given by: $E=\frac{1}{2\pi}\int_{-\infty}^{+\infty}|F(\omega)|^{2}\mathrm{d}\omega$.The value of $W$ is the 70% Bandwidth.


Model Validation After Estimation for System Identification Task (Assistance with MATLAB Code) QUESTION: I want to determine how well the estimated model fits to the future new data. How do I validate the estimated model...what is the procedure? After system identification, how to do model validation with the channel estimates? Please help.

In this example, I am using Least Mean Squares for system identification.
The system has coefficients / impulse response given by $\mathbf{h}$ which are estimated. The input signal to the system is $d$. The noisy output is $x$ where, the noise $e$ is Additive White Gaussian. I am trying to estimate the coefficients -- system identification using Least Mean Square and estimating the input $d$ as well. I am using an equalizer for this.

Thank you.

$x_n = r_n + e_n = \mathbf{h}^T \mathbf{d} + e_n$ where $\mathbf{h}$ represents the coefficients, $d$ is the input and $e_n$ is Additive White Gaussian Noise.

%% Channel and noise level
h = [0.9 0.3 -0.1]; % Channel
SNRr = 10;              % Noise Level

%% Input/Output data
N = 1000;               % Number of samples
Bits = 2;               % Number of bits for modulation (2-bit for Binary modulation)
data = randi([0 1],1,N);        % Random signal
d = real(pskmod(data,Bits));    % BPSK Modulated signal (desired/output)
r = filter(h,1,d);              % Signal after passing through channel
x = awgn(r, SNRr);              % Noisy Signal after channel (given/input)

%% LMS parameters
epoch = 10;        % Number of epochs (training repetation)
eta = 1e-3;         % Learning rate / step size
order=10;           % Order of the equalizer

U = zeros(1,order); % Input frame
W = zeros(1,order); % Initial Weigths

%% Algorithm
To get this we have to find error between input and output. As I understood `x` is an input for our model and `y` is an output.

    %% Algorithm
    for k = 1 : epoch
        for n = 1 : N
            U(1,2:end) = U(1,1:end-1);  % Sliding window
            U(1,1) = x(n);              % Present Input

            y(n) = (W)*U';             % Calculating output of LMS
            e(n) = x(n) - y(n);           % Instantaneous error
            W = W +  eta * e(n) * U ;  % Weight update rule of LMS
            J(k,n) = e(n) * (e(n))';        % Instantaneous square error
        end
    end

UPDATE : Please find below the picture that shows the step at which I am stuck. In the picture, the plus sign is the addition of Additive White Gaussian noise. I think the outcome of system identification of the unknown system is that the system should respond by producing the same output when given the same input. Or to predict. I don't want the MSE.

I agree with AlexTP's idea of simply plotting the instantaneous error between the model's output and the system's output over time. The exact implementation is a bit tricky though, because you are not only estimating the impulse response (using LMS) but also doing an "inverse" filter operation to estimate the unknown inputs. One way is to use a standard machine learning approach outlined in this answer.

Split your known input data $\mathbf d$ into two halves - training and testing (or validation$^*$).
Using the training half, run LMS and estimate the channel coefficients $\hat{\mathbf{h}}$.
Now use your test set and apply them to the noisy system to obtain outputs $x_n = \mathbf h^T \mathbf d + e_n$ and also apply them to the filter learned in step 2 to get $\hat x_n = \mathbf{\hat h}^T\mathbf d$. Compute the error$^\dagger$ between $\hat x_n$ and $x_n$ and make sure it's acceptable/small.$^\&$
If you are happy with the performance in step 3, your filter $\hat{\mathbf h}$ is ready for real world where you use the inverse of the filter $\hat{\mathbf h}$ to estimate the unknown inputs.$^{\ddagger}$

$^*$ Since you aren't doing model selection, you don't need a validation set, or in other words, your test and validation set are the same. If you do care about model selection (eg. choosing the best filter order) you can split your data into 3 pieces - training, testing and validation.

$^\dagger$ I don't think there's a single best way of computing this error. You could simply plot $\hat x_n - x_n$ vs $n$ and make sure it stays close to zero. Or perhaps you could aggregate them as a mean absolute error or mean squared error by summing over all $n$.

$^\&$ Alternatively, you can use the inverse filter of $\hat{\mathbf h}$ and estimate the inputs and compare them with the known inputs. Again, there's no single best way of showing error and would depend on your application.

$^{\ddagger}$ Since this phase has unknown inputs you can't really compute a performance metric for your model. Or as this answer puts it "you can only speculate about the quality of your model output using the results of your validation phase".


Recommended Order of Performing Denoising, Deblurring and Super Resolution on an Image In the normal workflow of image processing, there are multiple operations, increase image contrast, image denoising, image deblur and image super-resolution. What are the sequence/order of performing these different operations. Are there any overlap among these operations?

If your image is modeled as an image which is noisy, blurry and heavily decimated the optimal thing to do is estimate the image given that model.

The model is well defined in @Laurent Duval's answer.
I'd remark that in most real world cases the blurring is spatially variant hence it can't be modeled by convolution (Well, it is a generalized convolution).

In practice it is too hard (Was?) to solve such things and like in most cases we take the greedy approach: Step by Step.

Since most algorithms are SNR dependent, it makes sense first to handle the noise by Denoising. Pay attention that if we assume the image is blurred we have spatially correlated noise which is more tricky to deal with.

Then I'd handle the Super Resolution / Deblurring.
Since those are closely related I'd certainly solve them in one step.
Yet again, in case it is important to do them separately, start with Deblurring then Super Resolution.


How to Deduce a Linear System's Impulse Response from a Set of Input and Output Signals? I want to know how to solve those types of problems.. is it by inspection ?

Consider the linear system below. When the inputs to the system $x_1[n]$, $x_2[n]$ and $x_3[n]$, the responses of the systems are $y_1[n]$, $y_2[n]$ and $y_3[n]$ as shown.

Determine whether the system is time invariant or not. Just your answer.
What is the impulse response?

EDIT:  Assuming a general case where the given inputs don't contain a scaled impulse like $x_2[n]$

There seems to be an image which is gone now and hence I might be missing something.

In order to declare if the system is invariant you should see if a delay of the input yields only a delay in the output.
In you case, does the input $x_1[n-m]$ yields $y_1[n-m]$ and etc...

As written by @DilipSarwate, By looking at the response for $ {x}_{2} - {x}_{1} $ one could see the system is neither causal nor time invariant:

consider that linearity requires that the response to $ {x}_{2} − {x}_{1} $ be
$ T \left( {x}_{2} − {x}_{1} \right) = T \left( {x}_{2} \right) − T \left( {x}_{2} \right) = {y}_{2} − {y}_{1} $.
Now, $ {x}_{2} − {x}_{1} $ is a pulse of amplitude $ 2 $ at $ n = 2 $ and the response to it is nonzero at $ n = 0, 1, 2 $ and so the system is not causal, and since $ {x}_{2} − {x}_{1} $ is just $ {x}_{2} $ delayed by 2 units while the response is not $ {y}_{2} $ delayed by 2 time units, the system is time varying, not time invariant.

If the input signals are band limited and their bandwidth is less than your system you won't be able to restore the impulse response.
You will be able only to get the response in the frequencies the input has energy in.
This could be done by frequency analysis of the input and the output.
If your system is indeed LTI the connection between input and output is given by convolution with the impulse response.
Convolution is multiplication in the frequency domain, hence you could easily get the impulse response (Again, only at frequencies the input has energy in).

Yet since this system isn't time invariant it doesn't have a single impulse response. It has a response to an impulse at $ n = 0 $ but it might have a different one for an impulse at $ n = 1 $.
Update
This is a nice case to show the commutative property of convolution.
Since $ y \left[ n \right] = \left( h \ast x \right) \left[ n \right] = \left( x \ast h \right) \left[ n \right] $ you can switch the roles between the "System" and the "Input Signal".
So this question is equivalent to Given known output and system response, how could one estimate the input of the system?
Well, this is a known problem called Deconvolution and in the case of no noise it is simple to solve.
As written above, one way to do so it write the problem in matrix form.


What is the cost of Hamming coding? We know that (7,4) hamming code has 4 message bits and 7 code bits and can detect any single bit error.  The (12,4) repetition code has 4 message bits and 12 code bits and can detect any single bit error.  I tend to think that since in hamming code we are using less number of parity bits, there should be some disadvantage vis a vis repetition code.  Is it true and if yes, what is the downside of hamming code? (One known disadvantage is that of encoding and decoding complexity.)

This is an important subject that (IMO) is poorly handled in many introductory textbooks. To start, assume that an uncoded system spends $E_b$ joules per bit and transmits at rate $R_p$ information bits per second.

How can you compare this system fairly with the Hamming (7,4) code? Making them operate at the same bit rate implies a bandwidth expansion of 7/4, which is not always possible or desirable. Making them operate at the same energy per bit means that the energy of each Hamming coded bit should be $E_c=4E_b/7$.

So, you could say: under the same bandwidth and energy constraints, the cost of using a Hamming (7,4) code is a reduction in the data rate by a factor of 7/4. The benefit is a slight reduction in bit-error probabilty (but quite a bit less than that predicted by the binary symmetric channel).

Regarding the repetition code: under the same assumptions, the repetion code has a rate reduction of a factor of 3, while offering zero gain in probability of error!

Under the assumption that $E_c=E_b$, then the repetition code has quite a bit to offer in terms of bit-error rate. However, the cost is that you're multiplying the energy consumption by 3 while reducing your data rate by 3. To be fair, you'd need to compare the repetion code with an uncoded system whose $E_b$ is also increased by a factor of 3.

To summarize: comparing codes is not trivial, but not terribly hard either, as long as you make all your assumptions clear.


How to calculate resolution of radar? I am working of FMCW radar. Below are all the parameters of FMCW radar. How can I Calculate the resolution of my radar

Parameters            Values
Center frequency (fc)     70GHz
Bandwidth (BW)            40 GHz
Sampling frequency (fs)   10 MHz
Sweep time            0.2 ms
Number of periods    10 period

The range resolution is normally some fraction bandwidth of the pulse. The fraction will depend on how you measure it i.e. null to null, 3dB points, 6 dB points etc. Also the windowing (Hamming, Kaiser) used will reduce the resolution i.e. increase the number. For a chirp pulse with the bandwidth you stated the null to null width is given by

$$ \delta r = \frac{c}{2BW}=\frac{3\cdot 10^8}{2\cdot  40\cdot 10^9 }=3.75 \cdot 10^{-3} \textrm{ m } $$

You can reference the following links:

High Range Resolution Techniques
Understanding Millimetre Wave FMCW Radars,
Radar Basics


Learning roadmap for the audio signal processing beginner I would like to start learning audio signal processing. There are numerous books online and academic papers all of which seem to skip the fundamentals of the topic.

I would like to know a rough roadmap, so to speak, to follow in order to successfully learn audio signal processing.

I have read that calculus is the first step before starting with signal analysis.

It feels to me like audio signal analysis is only one part of the overall knowledge needed. Where other topics are music theory, audio engineering and programming.

If I may ask people with knowledge on this area to suggest possible steps towards understanding how to analyse and manipulate/create audio signals.

I recommend having a look at Prof. Julius O. Smith III's Physical Audio Signal Processing. It's available online, or can be purchased through Amazon's print-on-demand service.

In particular, the description in the Book Series Overview might be worthwhile.


How to generate IQ component for transmitting 2FSK I am trying to generate IQ components for 2FSK modulation.
I have started off from this thread FSK and IQ modulation , but I don't quite understand the IQ TX components plotted in the example there by GordonFreeman.

In my view, in a 2FSK modulation, the angle of I+jQ shall vary linearly with time (constant frequency +/- Δf), drawing a circle either counter-clockwise with frequency Δf if TX data is equal to 1, or clockwise with frequency Δf if the TX data is equal to -1 (given e.g. a random -1/+1 sequence). In this convention, the difference between high and low freqs in 2FSK is 2*Δf.

This is well explained from about 07:00 in the following video: https://www.youtube.com/watch?v=5GGD99Qi1PA (can be played at faster speed if convenient).

If I understand correctly, for 2FSK modulation:

I(k)=cos(2*pi*freq(k)*t)

Q(k)=sin(2*pi*freq(k)*t)

with e.g. freq(k)=Δf*k where transmitted symbol k=-1/+1 if we have a binary data stream of +1/-1. In this sense, I(k) will be just a cosine wave, while Q(k) will be a sine wave which will be pi shifted every time there is a data change. An example are below IQ components I generated by a pyscript:

The above components are generated when bitrate is maximum(4*Δf).

Does this make sense?

=====

Thanks for your comments Dan.

Some more info-  I started off from a 2FSK modulated signal of this form:

modulated(t)=cos(2*pi*(fcarr+freq(k))*t)

where fcarr is the carrier frequency, freq(k) is either +Δf or -Δf as discussed depending of the transmitted data.

Then deriving the below I/Q components
(those that fed into a quadrature modulator would lead to the above modulated signal, I omit the full derivation but that can be derived either in freq domain https://en.wikipedia.org/wiki/Single-sideband_modulation or more simply in time domain with trigo formulas https://en.wikipedia.org/wiki/Quadrature_modulation):

I(k)=cos(2*pi*freq(k)*t)

Q(k)=sin(2*pi*freq(k)*t)

Therefore when freq(k) is changing from +Δf to -Δf, the I(k) which is a cosine will remain unchanged (if we multiply by -1 the arg of a cosine, the cosine won't change), while Q(k) which is a sine will change of a 180 degrees since sin(-x)=-sin(x).
In this sense, the rotation of the I/Q vector will be changing from counter clockwise to clockwise, but there will be a sudden phase shift of 180 degrees.

If maximum bitrate = 4*Δf is choosen, then duration of 1 symbol is 1/fourth of Δf period, or 90 degrees in the I/Q vector diagram.

So I/Q vector will start from (1,0) i.e. 0 degrees, moving counterclock wise to (0,1) i.e. 90 degrees, that will be within the duration of symbol "1" , then there will be a symbol change to "-1", that will produce, as mentioned above, just a change of I component from sin(2*piΔft) to -sin(2*piΔft). This will make the I/Q vector moving from (0,1) to (0,-1) and the vector will move now clockwise to (-1,0) for the duration of symbol "-1".

From this, the above I/Q plot was generated.

Your I & Q plots look wrong to me.  Instead of trying to sort out what you may have done wrong, let me just give you a working example instead.

The following GNU Octave (Matlab clone) script shows how to generate the I & Q of 2-FSK with square pulses and a modulation index of $1/2$ (which should correspond to a bit rate of $4\cdot\Delta f$).

The "magic" happens on the lines following the comment "Frequency modulate the baseband packet"

% Keep Octave happy
pkg load signal;

% Bit rate
Rb = 9600;

% Modulation index; 1/2 is Minimum Shift Keying.
modulation_index = 1/2;

% Optional frequency shift for demonstration of modulation on a
% carrier within -Fs/2 to +Fs/2
freq_shift_hz = 0;

% Samples per symbol and sample rate out of the modulation process
sps = 10;
Fs = Rb * sps;

% pulse filter taps - square pulses
taps_pf = ones(sps, 1)/sps;

% Create a random bit string
nbits = 16;
message_bits = round(rand(nbits,1));
nrz_message_bits = (message_bits - 0.5)*2;

% Upsample and pulse shape the packet
% N.B. Octave's filter function sets up a history of 0's for us.
x = sps*[upsample(nrz_message_bits, sps); zeros(1, 1)];
packet_baseband = filter(taps_pf, [1], x);

% Throw in a frequency shift (as a level shift of the unscaled baseband)
% for demonstration purposes
freq_shift = freq_shift_hz / (Rb/2 * modulation_index);
packet_baseband = packet_baseband + freq_shift;

% Frequency modulate the baseband packet
fm_gain = pi/(Fs/2) * Rb/2 * modulation_index;
x = packet_baseband * fm_gain;
phase = cumsum(x); % phase is integral of frequency
packet_modulated = exp(1i*mod(phase, 2*pi));

figure(1);
stem(message_bits);
title('Message Bits');
xlabel('Bit number');
ylabel('Value');
grid on;

figure(2);
N1 = length(packet_baseband);
t1 = [0:N1-1]/Fs;
plot(t1, packet_baseband * Rb/2 * modulation_index, 'x-');
title('Upsampled, Pulse Shaped, Scaled, and Level Shifted NRZ Bits');
xlabel('Time (seconds)');
ylabel('Frequency Deviation (Hz)');
grid on;

figure(3);
N2 = length(phase);
t2 = [0:N2-1]/Fs;
plot(t2, phase, '.-');
title('Accumulated Phase');
xlabel('Time (seconds)');
ylabel('Accumulated Phase (radians)');
grid on;

figure(4);
N = length(packet_modulated);
t = [0:N-1]/Fs;
plot(t, real(packet_modulated), t, imag(packet_modulated));
title('I & Q of 2-FSK Modulated Signal');
xlabel('Time (seconds)');
ylabel('Amplitude');
grid on;


Balance Fader in Digital Amplifiers What is the purpose of Balance / Fader in a digital Audio amplifier?
In a signal flow what is the most apt position to place balance and fader(for example after before a compressor unit or at the starting of the signal flow or at the end of signal flow) and why?

What is the purpose of Balance / Fader in a digital Audio amplifier?

In a stereo amplifier the balance controls allows to set different gains for the left and right channel. That would be useful if the listener is not located right on on the center line between the speakers or if the acoustic loading the two speakers is very different. The balance control can be used to reposition the image in front of the listener.
A "fader" is more common a car. In addition to being able to balance between left and right a fader can also balance between front and rear as well. This control is mostly outdated: in a well tuned car, a fader will do more damage than good and many cars these days already have "front" "driver" "rear" etc. modes which are individually optimized for the use case.

In a signal flow what is the most apt position to place balance and fader(for example after before a compressor unit or at the starting of the signal flow or at the end of signal flow) and why ?

That depends a lot on what other functions the amplifiers has. It's mostly "in the middle", after the "flavoring" controls (different modes, tone controls, dynamic loudness, dynamic range adjustments, ...) but before the speaker management (cross over, EQ, limiter, bass manager, etc.).


Is there an asymmetric crypto analog to direct sequence spread spectrum? In section 12.1.1.2 of Digital Communications: Fundamentals and Applications, Bernard Sklar makes this claim about spread spectrum signals:

  Therefore, not only can the spread-spectrum signal be made difficult
  to jam, but additionally, the signal's very existence may be rendered
  difficult to perceive. To anyone who does not possess a synchronized
  replica of the spreading signal, the spread-spectrum signal will seem
  "buried in the noise."

This seems analogous to symmetric cryptography; there is a symmetric "key" (the spreading signal) that allows one to both transmit and receive messages, and without the key one can neither transmit nor receive a message. Furthermore, without the key one cannot even detect the presence of a message (or, at the very least, it makes it more difficult to do so).

This makes me wonder if there's an asymmetric cryptography analog to direct sequence spread spectrum. I assume there would be a public and private "key" pair, and one could use the public key to transmit a message, but only the private key could be used to receive that message. And without the private key, it would be very difficult to even detect the presence of a coded message (i.e. only the private key could be used to "de-spread" the message).

Does such a technique exist?

(This does not answer your question which is a good one, so I may move this to a comment below your question)

Warning! Typical Direct Sequence Spread Spectrum that uses linear spreading codes (based on linear feedback shift registers implemented from irreducible and primitive polynomials) are NOT a good cryptographic technique. See this post where the approach to break such codes is detailed:

Minimum Output Samples needed to crack a "Gold Code" Generator (Dual LFSR)

This is because although the signal will indeed be "buried in the noise" at certain reception points, it is very rare to not be able to get closer to the transmitter (either spatially or with highly directional antennas) such that a positive SNR can be established enough to get the relatively small number of consecutive sequences to crack the code.

Use of non-linear codes together with data layer encryption are recommended. I wanted to mention this due to the common misconception that DSSS alone is a good communication security technique.


Channel Vocoder producing output with "click" sounds I'm trying to make a channel vocoder that takes two inputs, one a frequency rich carrier (a musical sound) and the other a modulator (vocals). Applying operations involved in the channel vocoder generates an output that makes the vocals sing in the tone of the carrier. (Think Songify.)

Presently I'm splitting both the signals into windows of a specified size and also overlapping two consecutive windows by a certain factor. FFTs of the windows of the carrier and modulator are then computed. The fft is further split into bands of equal length. The average of magnitudes of modulator's fft in that particular

for band_no = 0 to band_count:
        for i = 0 to band_length:
                avg+=magnitudeOfModulatorFFT[band_no*band_length + i];

I am then multiplying "avg" to each of carrierFFT's frequency value lying within the band, and then taking the inverse FFT. Although this produces the desired output, it also includes a clicking sound. Even after multiplying both the signals by a window function.

Any way I could solve this?

I think it could be one of two issues.

What is your frame overlap? if the frame overlap is something less than 40% or so you will begin to get undesired artefacts when you perform overlap add synthesis at the end. Try an overlap of 50-80% (although more computationally expensive, it provides greater frequency resolution).

The second consideration which goes hand in hand with the above suggestion. Are you using any kind of analysis window on your framed signal? Perhaps try using a tapered window such has hamming, Blackmann or Hanning. This will help minimise spectral leakage between frames. Hope this helps.


Discrete Fourier Transform on Impulses? So I'm used to DFT being on an equation like this:

x(n) = cos(n), ∀n , with window n = 0,1,2

I can clearly plot discrete values and take a DFT to find Xk

But I found this equation on a past exam that's like this

x(n) = δ(n+2) +δ(n+1)−δ(n−1) with window of length 5.

This appears to be an impulse response to me, and I don't understand how one would take a DFT of it?

For digital signals, you represent the $s(n) = \delta(n)$ dirac impulse by

$$s[n] = \delta_\text{discrete}[n] = \begin{cases}
0& \forall n \ne 0\text,\\
1 & n = 0
\end{cases}\tag{*}$$

which is obvious if you think about the dirac impulse of having the function (in continuous math) to give you

$$\int\limits_{-\infty}^\infty f(x)\delta_\text{continuous}(x) dx = f(x) $$

and in the discrete case, your integral collapses to the sum of

$$\sum\limits_{-\infty}^\infty f[x]\delta_\text{discrete}[x] dx \overset!= f[x] \tag{$\dagger$}$$

The only function that fulfills $(\dagger)$ is $\delta$ from $(\text *)$.


Why should an image be blurred using a Gaussian Kernel before downsampling? I recently read that before downsampling an image, it should be blurred using a Gaussian Kernel. This way, the downsampled image is better than just picking a single pixel out of a NxN block or averaging over the block. After searching in this site as well as google, I didn't get any exact answer.

But there were questions on how to select $\sigma$ for blurring. Reading the answers on those posts, I learned that downsampling has to be done in accordance with sampling theorem. Downsampling without blurring causes aliasing effects.

Can someone please explain why image has to be blurred before downsampling? I mean what is the exact relation to sampling theorem. What happens when an image is downsampled without blurring? I mean what are these aliasing effects? How can I notice them in the downsampled image?
Why is Gaussian blurring better than Averaging over a block?

If you can give some examples of the images, I would be a lot grateful. Ans, I would appreciate all kind of answers, partial, intuitive, rigorous, anything.

Thank You!

An image "should not be blurred using a Gaussian Kernel" in general. This can be a safe bet for a lot of basic image processing needs, and a smoothing is almost mandatory when you want to control the information lost by the downsampling.
Blurring is (often, not always) another word for low-pass filtering. When an image contains high-frequency content (fast variations), downsampling without blurring can produce visually weird, annoying aliasing artifacts or moiré patterns. There is an example on wikipedia Aliasing: the original

and the downsampled, represented at the same size:

The ripples on the bottom right are low-frequency artifacts generated by a careless brute-force downsampling. A blurring would attenuate image sharpness, dim the borders between bricks, and reduce the apparent aliasing aspects.
The choice of appropriate blurring filters has a long history in image processing. Gaussian shapes have long been considered somehow optimal for different reasons for "theoretical" continuous images. Plus, it is both decreasing and symmetric in  the space and the frequency domains. In the time domain, this means that faraway pixels have less influence. In the frequency domain, frequencies are reduced monotonously from low to high.
Since most images are discretized, reality is somewhat different. Since the Gaussian convolution  used to be computationally expensive, early approximating filters were designed with short support, borrowed for instance from Pascal triangle. Later, fast recursive implementations were designed (Deriche, Shen, etc.)
I guess question 1) is answered. For question 2) simple averaging gives an equal weight to all pixels in the window. Hence, faraway pixels are given equal importance with respect to closer pixels, which is not optimal in regions where images exhibit weak stationarity, like trends, edges and textures.


Why doesn't this complex multiplication in the frequency domain produce my expected phase shift? I know how to change the phase of a complex number by multiplying by $\cos \theta + i \sin \theta$.  And I understand that the phase of a sine wave is reflected in its Fourier transform.  So, I am trying to phase-shift a signal by changing the phase of its Fourier transform..

This works for "synthetic" Fourier transforms, but when I try to FFT a signal, apply the phase change, and invert the FFT, I don't get the expected result.  I am using MATLAB in the examples below.

Fs = 1000;
Tmax = 2;
L = Tmax * Fs;

For example, I'll build a synthetic Fourier sequence.  The frequency intervals are 0.5 Hz, so a 10Hz component is at the 20th position after DC.  (I'm just spitballing here; I know there's better ways to pick frequency bins than hand-jamming 21.)

Ysynth = zeros(1, L);
Ysynth(21) = 1000;

And rotate it by 45°:

Ysynth = Ysynth * exp(j * pi/4);

It has a 45° phase on the 10-Hz component:

polarscatter(angle(Ysynth), abs(Ysynth))

Now take its inverse FFT:

Xsynth = ifft(Ysynth);

So far so good.  The generated signal has a phase shift of 45°.  MATLAB does complain about the presence of an imaginary part when I plot it; I think this is because I didn't bother with the negative frequency component.

t = (0:L-1)/Fs;
plot(t, Xsynth)

Now, I would expect the same principle to apply to a frequency spectrum taken from an actual signal.  But I do not get similar results.  Quickly:

X = cos(2*pi * 10 * t);
Y = fft(X);
Yshift = Y * exp(j * pi/4);
Xshift = ifft(Yshift);

There is an imaginary component in Xshift, which I do not expect, and again MATLAB complains about.

plot(t, [Xshift; X])

There is no phase shift here, and the amplitude (of the real part) is different.

I must be misunderstanding something about phase representation in FFTs.  Why doesn't my transformation produce a phase-shifted version of the original?

If you want a phase-shifted strictly-real result, then you have to make sure the data you feed to the IFFT is conjugate symmetric.  So make sure to reflect the complex conjugate of any phase changes you make.  e.g. Make the upper or negative half of your IFFT input vector mirror the complex conjugate of the bottom or positive half of the input vector.

Added: (Also, "mirror" usually means the array index goes in the opposite direction.)


I trigger recording when power of current sample is X dB greater than last Y samples, what is the terminology here? Doing time domain processing of incoming audio. I sample the power level every 0.25s and keeping a running average of intensity over 5 samples. When the power level of the current sample is 10dB greater than the moving average, we start recording.

What are the correct names for these numbers? I also use the same method to trigger recording stop. Is this the same terminology?

I have the code and will use your suggestion to name the variables correctly.

Thank you.

"Running average" is fine. You can also call it the moving average.

I suppose you could call the variable associated with the 10dB number the signal detection threshold.


Group delay of non differentiable unwrapped phase As we know group delay as a function of frequency, computed as differentiation of  unwrapped phase, but what happen  at those frequencies where unwrapped phase is non differentiable?how to calculate delay at that particular frequency?

If there are phase discontinuities after all $2\pi$-jumps have been removed, then these discontinuities are usually caused by zeros of the frequency response. The phase jumps by $\pi$ at these frequencies, and the group delay doesn't exist, or, if one prefers, is non-finite. Note that the group delay is meaningless anyway at frequencies where the frequency response is zero.


The notion of white and colored noise The autocorrelation function of AWGN is an impulse response at t=0, and its fourier transform or its spectral power density has a flat infinite BW

Now any receiver will have bandpass or low pass filters that will limit the noise BW

Does that mean that in reality we are actually dealing with "colored" noise and considering noise to be "white" within the filter BW is an assumption made for convenience purpose?

Yes.

However, if you are dealing with a discrete-time signal, if your A/D is proceeded by an ideal low-pass filter with cutoff frequency $F_s/2$, with unity passband gain, the output power spectral density before sampling is:
$$
S_y(f) = |H(f)|^2 S_x(f) = \frac{N_0}{2} \operatorname{rect}(f/F_s)
$$
This means that the autocorrelation of the filtered noise will be:
$$
R_y(\tau) = \mathscr F^{-1}\{S_y(f)\} = \frac{N_0}{2T_s}\operatorname{sinc}(\tau/T_s)
$$
where $T_s = 1/F_s$. If the noise is zero-mean, Gaussian, and sampled at time instants $t_n = n T_s$, then each sample will be uncorrelated, and since they are Gaussian, will also be independent.   In this case, in the digital domain, the power spectrum will be flat for all frequencies.

This is not the only way to generate noise that is white in the discrete-time domain, but is one way, albeit a physically unrealizable one.


Block diagram of filter Can you help me with this question. Can you explain where i should start to solve the problem.

The block diagram you have posed is composed of serial cascade of three IIR sections each having its own input output relationship. The first two are second order stages and the last one is first order.

If you know about Z-transforms, you could easily find the algebraic relation between those signal variables in terms of their Z-transforms.

I would denote the intermediate summation nodes as $w[n]$ and $v[n]$ together with the input $x[n]$ and output $y[n]$ the Z-domain relations become:

$W(z) (1 - a_1 z^{-1} - a_2 z^{-2}) = X(z)$
$V(z) (1 - b_1 z^{-1} - b_2 z^{-2}) = W(z)$
$Y(z) (1 - c_1 z^{-1} ) = V(z)$

where I have denoted the system coefficients as $a,b,c$ respectively. (Do not mix this $a,b$ with the standard LCCDE expression coefficients $a,b$.)

Then simple algebra tells us the relation between $X(z)$ and $Y(z)$ as:

$$Y(z) = 2 \left( \frac{1}{1 - c_1 z^{-1}} ~ \frac{1}{1 - b_1 z^{-1} - b_2 z^{-2}} ~\frac{1}{1 - a_1 z^{-1} - a_2 z^{-2}} \right) X(z)$$

By first expanding the product polynomials and then mapping back from Z-domain to time domain you can express $y[n]$ as a function of $x[n]$ (and feedback)


sparse representation for image denoising When I read papers on image denoising, I always encounter sparse representation. For image denoising, we try to separate image signal from noise. It is assumed that signal is correlated and noise is uncorrelated. Sparse representation represents one signal as a linear combination of a small number of dictionary elements. It hope to use as few as non-zero coefficients to represent signal. In my view, it is appropriate for image compression. How can it is related to image denoising?

Here are two interpretations of why sparse representations are useful for denoising:

1) Your clean image is well represented by your dictionary (usually a wavelet or 2D-DCT dictionary), but the noise isn't. This means that when calculating your sparse representation from the noisy signal, you will only capture the clean signal, and not the noise.

2) Enforcing the representation vector to be sparse (i.e. with only a few non-zero elements), means that you will eventually threshold out all the other low-energy noisy components.

As you said sparse representations can be pretty useful for compression. Compression is merely about keeping the relevant information, while discarding the irrelevant one. Well denoising is about keeping the relevant part of your signal (the one that fits your model, or your dictionary), while discarding the irrelevant one (the noise).


Detect keypoints in images without corners I am trying to achieve a pattern matching between an image and a scene in which the image is. I am having problems with some images that does not contain corners, since SURF is not able to find keypoints.

Does anybody know any algorithm or method for retrieving keypoints in images without corners? In other words: how can I find an image inside another when this image does not contain corners?

I am using Matlab. Thank you in advance,

The computer vision system toolbox in Matlab includes Feature Detection, Extraction, and Matching. The line, edges, histogram of oriented gradients are all included other than the corner features. And there is also an example on point feature matching in the toolbox (the names are called corners but they are actually not).


FFT of a mirrored array according to the definition of the Fourier transformation we have for a real function:

if $F{{(a(t))}}$ = $A(f)$
then $F{{(-a(-t))}}$ = $-A$*$(f)$, with * for conjugate

import numpy as np
a= np.array([1, 2, 3,4,5,6,7,8,9,10])
a_mirrored = np.array(list(-a[::-1]))

A_f = np.fft.fft(a)
A_f_star = -np.conjugate(np.fft.fft(am))

So why is A_f different from A_f_star ? Can anyone tell please?

This is occurring because you have not time-reversed the signal properly. The time reversal intended by your identity is for $X[n] = X[(N-n)_N]$ where the notation $(\cdot)_N$ means that you should take the result modulo $N$. In other words $X[0] = X[(N)_N] = X[0]$ doesn't change position.

For your example, you start with the finite-length signal $a = [1, 2, 3, 4, ... 9, 10]$. The true time reversal of $a$ is $a_r = [ 1, 10, 9, 8, ... 3, 2 ]$. But the python code results in simply reversing the order of the components without respect to their indices, thus the result is the signal $a_m = [ 10, 9, 8, 7, ... 2, 1 ]$.

The reason the first entry is the same in both cases is because it is simply the sum of the elements which is equal for all three signals $a$, $a_r$ and $a_m$.

Here is a function that will time reverse your signal.

def time_reverse(x):
    xrev = [ x[0] ]
    xrev.extend(x[-1:0:-1])
    return xrev

So, you can add this function definition to the top of the file and replace the line

a_mirrored = np.array(list(-a[::-1]))

with the line

a_mirrored = -np.array(time_reverse(a))

Then your plots should match up fine.