How to demodulate an AFSK signal in software I am trying to transmit binary data from one device to another over an audio channel (speaker/mic). I use AFSK (Audio Frequency Shift Keying) as in Packet Radio, with $1200 \text{ Baud}$ and two frequencies $f_{mark} = 1200 \text{ Hz}$ and $f_{space} = 2200 \text{ Hz}$. I played around a bit in Ruby and my first implementation simply imitates a classic incoherent demodulator, which works fine, so far.

The problem is, i am trying to port this to a mobile platform where performance is a concern and my current solution is too slow. I have found numerous ways to demodulate AFSK in software:

Sliding DFT (FFT)
Sliding Görtzel Filter
Phase Locked Loop
Zero Crossing

What would be the way to go? There are just too many options to choose from. I am sure there are even more options available. Perhaps there exist even better solutions than the ones i named above? Does someone even have a code examples for me? I am concerned with

Performance (should run on mobile Platform, say an iOS or Android device)
Stability (should be able to handle some noise)

Any suggestions and hints are greatly appreciated!

I think that you could get the best performance in terms of demodulator bit-error rate (BER) with a phase-locked loop.  You need it to be fast, though.  I think your best bet for a fast algorithm that still performs reasonably well is zero crossing.

On a side note, I would like to suggest that you change the 2200 Hz to 2400 Hz.  A naive implementation of the 1200/2200 Hz scheme would yield discontinuities, as seen about two thirds into the plot below, where the 2200 Hz transitions to 1200 Hz.

In order to minimize the bandwidth that you're using and avoid discontinuities that will distort the signal you'll need to make the phase continuous.  Even if you make the transmitter phase continuous, though, there will still be the issue that the 2200 Hz symbols will not always have the same number of zero crossings due to the different phases.   Usually they will have four zero crossings, but sometimes they will have three.  The 1200 Hz symbols on the other hand, will always have two zero crossings because the baud rate divides evenly into the FSK frequency.

You can solve both of these problems by changing the 2200 Hz to 2400 Hz.  Then the symbols will always start and end at 0 degrees (thus automatically making them phase continuous), and they will always have the same number of zero crossings- two and four.


Determining the Filter Coefficients of an FIR Filter I need to find the filter coefficients of an FIR filter that will block sinusoids of frequency $200\ \rm Hz$ if the sinusoid is sampled at $1.2\ \rm kHz$.

I feel like this is a fairly simple problem, but I'm not exactly sure how to go about it.

I am thinking I need to convert the frequencies to rad/s, then use these values in the transfer function, but am not sure how to find this transfer function.

EDIT:

I have to choose from one of the possible options:

$$ \{1, 1, 1\} $$
$$ \{1, -1, 1\} $$
$$ \{1, 0, 1\} $$
$$ \{1, \sqrt{2}, 1\} $$

or None of the Above.

From the coefficients above, I know the following equations can be created:

\begin{align}
H\left(e^{j\omega}\right)&=1+e^{-j\omega}+e^{-2j\omega}\\
H\left(e^{j\omega}\right)&=1-e^{-j\omega}+e^{-2j\omega}\\
H\left(e^{j\omega}\right)&=1+e^{-2j\omega}\\
H\left(e^{j\omega}\right)&=1+\sqrt{2}e^{-j\omega}+e^{-2j\omega}
\end{align}

The general formula for the frequency response of a causal second-order FIR filter with a pair of complex zeros on the unit circle is

$$H(e^{j\omega})=G(1-e^{j\omega_0}z^{-1})(1-e^{-j\omega_0}z^{-1})\tag{1}$$

where $\omega_0$ is the (normalized) frequency in radians where the zero occurs, and $G$ is a (real-valued) gain constant. If we choose $G=1$ and expand $(1)$ we obtain

$$H(e^{j\omega})=1-(e^{j\omega_0}+e^{-j\omega_0})z^{-1}+z^{-2}=1-2\cos(\omega_0)z^{-1}+z^{-2}\tag{2}$$

In your example we have

$$\omega_0=2\pi\frac{f_0}{f_s}=2\pi\frac{200}{1200}=\frac{\pi}{3}\tag{3}$$

where $f_0$ is the frequency of the zero, and $f_s$ is the sampling frequency. Plugging $(3)$ into $(2)$ with $\cos(\pi/3)=\frac12$ gives

$$H(e^{j\omega})=1-z^{-1}+z^{-2}\tag{4}$$

which corresponds to the impulse response $\{1,-1,1\}$.


Why is the ideal low pass filter not achievable? I find myself to be quite confused by some concepts of DFT. Suppose I have some spectrum-normalized signal (which I simply receive without any prior information) $$A\cos\left(0.5n\right)\ +\ B\cos\left(0.7n\right)\ +\ C\cos\left(1.2n\right)\ +\ D$$
with A,B,C,D being some constants
Now I find that the DFT of the same signal corresponds to a different set of basis functions given by $e^{\frac{2\pi jn}{N}}$ for N being the number of data points.
Now I can retrieve the signal by multiplying each "complex frequency component" with its "complex amplitudes" and adding them. Now my question is as follows:
If I apply an ideal low pass filter, I should get an infinite sequence of data points because the sinc function would have infinite points in it and convolution of it with an infinite set of data points would be again infinite and double-sided. (I hope this is a correct guess??).
But if I simply had a mechanism to be able to separate this unknown sequence into the above form, (say some alternate transform that would give me this decomposition instead of the exponential decomposition) then will I be able to design the "Ideal low pass filter" in that case just by knowing the above decomposition? I feel it can be proved (though the uniqueness I'm not sure) that such a decomposition will exist because it would be a set of N polynomials in the cosines of the found natural frequencies, like for the mth data point$$y\left(m\right)\ =\ \sum_{i=1}^{N}\cos\left(m\alpha_{i}\right)$$
where $\alpha_{i}$ is the ith frequency. Here the frequencies are normalized so the +$2\pi$ for each frequency idea cannot be used to refute the uniqueness
Also if I knew that the signal had exactly N points, then the signal could be represented at most by N real frequencies. Then why is the ideal low pass filter not realizable?
I sincerely apologize if the question was very naive, it could potentially be so.

General Background: Representing a signal numerically in a computer requires the signal to be discrete. That means it's periodic in the other domain. When you use something like a DFT you need the signal to be discrete in both domains: That means the signals are ALSO periodic in both domains.

If I apply an ideal low pass filter, I should get an infinite sequence of data points because the sinc function ...

If you want to do this numerically in a computer, it's NOT a sinc function. The lowpass filter MUST be periodic in both domains and so what you get in the time domain is a periodic repetition of a sinc function with a time domain aliasing.
It's still infinite it time. In fact it MUST be infinite in time since it's periodic. As long as all signals involved are actually periodic with a shared Period, than things work "as expected" .
At this point the term "low pass" filter becomes complicated. Since the spectrum is periodic, it will pass arbitrarily high frequencies. So it needs to be interpreted: it removes higher frequencies in the base band or something like this.

$y\left(m\right)\ =\ \displaystyle\sum_{i=1}^{N}\cos\left(m\alpha_{i}\right)$

Your missing phase or sine term here but you are re-inventing the DFT or cosine transform. You run into exactly the same problem. How to choose N ?

Also if I knew that the signal had exactly N points

That's a contradiction. The signal that you described in your original equation is infinitely long. As soon as you make it finite in time, you are windowing it and that changes the spectrum and the math behind considerably.

then the signal could be represented at most by N real frequencies

Correct. That's for example the Discrete Cosine Transform, which is a close cousin of the DFT. But it makes no material difference here.

Then why is the ideal low pass filter not realizable?

It is realizable as long as everything periodic "the right way". The periodicity helps dealing with the infinite length in time: you need only to compute what's happening in one period and you have all information that can be had since all other periods are the same.
However most real world signals are NOT periodic "the right way", so ideal low pass filtering is rarely feasible.


Why is the bandwidth for this QPSK system larger than expected? Below is my gnuradio-companion graph.

As you can see the sample rate is 100k, I am modulating with QPSK, and am filtering the waveform (RRC) with the "Excess BW" parameter, at 0.35.  The spectrum is shown below:

The bandwidth seems to be about 2*160kHz.  I would expect it to be about the baud rate, but it is far from that.  What am I doing wrong?

In gnuradio, the throttle block only determines the speed at which your computer is processing samples - it does not determine what the actual sample or baud rates are. As an example - consider that you generate random PSK symbols at 1 sample per symbol, and then upsample by 2 with a root raised cosine w/ 35% exccess BW so we're now at 2 samples per symbol (this is what the constellation modulator block does in your flowgraph). What is the absolute baudrate that you're running at? Answer: Whatever you want it to be. If I hook that I/Q stream up to an RF front end, I could theoretically transmit at any baud rate using those exact same samples. E.g. if the DAC is running at 1MSPS (assuming this is an I/Q DAC for simplicity like you might consider the ad9361) then my baud rate would be 1MSPS/(2 samples/symbol) = 500k symbols/sec. Assuming my PC could keep up (and I increased my throttle block so that my PC isn't slowing me down from producing samples), I could instead run the DAC at 10MSPS to transmit at 5M symbols/sec.
So - your question will probably be - how do I simulate a certain baud rate without a DAC/front-end? Well - it's whatever you choose it to be. In your flowgraph, you set the FFT sink to think that the sample rate is 500 kHz. The data you have hooked up to that sink is coming in at 2 samples/symbol, but factoring in the excess bw, the relative bandwidth will be about (1.35/2.0)*500e3 = 337.5 kHz which is exactly what the spectrum is showing you. If you did further upsampling and then looked at the FFT, the signal would appear narrower and vice versa.


System identification/ Filter estimation to mimic frequency equalizer of audio with Scipy At the current problem I'm working on, I have two signals: One "original" signal that contains audio (voice). The second signal is the same audio file but edited with a frequency equalizer, for example to make all frequencies between $649-677\textrm{ Hz}$ louder.

What I'm trying to achieve now is to "teach" a SCIPY filter to mimic this behavior. Therefore I created a function with two butterworth-filters and try to curve_fit the $x$-data to the $y$-data. The intuition is, that like this the system learns the optimal frequency intervals to mimic the equalizer. Obviously it doesn't work as supposed (no curve fit achieved most of the time. If I achieve a curve fit and apply the filter to the original data, it didn't change at all). Any suggestions?

import numpy as np
from scipy.optimize import curve_fit
from scipy import signal

xdata = voice_1
ydata = voice_2

# Function to optimize, supposed that two bandpass filters are enough
def func(x, p1,p2,p3,p4):
    b, a = signal.butter(2,[p1,p2],btype="band")
    x = signal.lfilter(b,a,x)
    b, a = signal.butter(2,[p3,p4],btype="band")
    x = signal.lfilter(b,a,x)
    return x

# Optimize function to achieve a curve-fit and therefore get the right frequency response in p1,p2,p3,p4
popt, pcov = curve_fit(func, xdata, ydata,p0=(0.5,0.5,0.5,0.5))

What I also tried was creating a xsignal containing of white-noise, creating a ysignal where I changed some frequency levels and calculating a transfer function H as FFT(ysignal)/FFT(xsignal). As soon as I apply this function to other data than white noise the results seem to be wrong as well.

Thanks for all your comments. I finally found the answer this article and implemented it in python:https://www.dsprelated.com/freebooks/filters/Time_Domain_Filter_Estimation.html It works pretty well, if you first find the transfer function h(t) with raw and equalized white noise and then apply the function as a filter to any audio data:

import numpy as np
import scipy

def learn_equalizer(signal_x, signal_y, detail_level):
    # Create Toeplitz Matrix
    padding = np.zeros(detail_level, signal_x.dtype)
    first_col = np.r_[signal_x, padding]
    first_row = np.r_[signal_x[0], padding]
    X = scipy.linalg.toeplitz(first_col, first_row)

    # Implement the formula for finding the optimal solution numerically
    X = X[0:len(X)-detail_level,]
    h = np.dot(X.T,X)
    h = scipy.linalg.inv(h)
    h = np.dot(h,X.T)
    equalizer = np.dot(h,signal_y)

    return equalizer

Finally looking at the curves let's consider two plots: One with the raw and equalized signal, to show their differences. The other plot shows the equalized signal and the "mimicked" equalizer. Here you see how good the learned filter is, when applied to real sounds after being trained on white noise.

Signal comparison:


Is CMA equalization applicable for OFDM I have a question regarding the constant modulus algorithm (CMA), as you know it's one of the known algorithms used for blind equalization. The cost function of the CMA is build blindly, it means based on received signal only. The cost function of that algorithm is:

$J = \frac12\left[\left\lvert w^T*y(n)\right\rvert^2 - R_{cm}\right]^2$

where $R_{cm}$ is the constant of CMS algorithm which is usually equals 1. $w^T$ is the equalizer coefficients and $y(n)$ is the received signal.

My question, is that algorithm applicable for OFDM system?

If so,  and what's  the constant $R_{cm}$ will be?

No.

OFDM isn't constant modulus (i.e. constant envelope) in time domain, if you look at it as one system. It's quite the opposite; it's known for its high PAPR (which you probably know!). This is the case, by the way, even if you use a CM modulation for the individual subcarriers.

So, CMA's central requirement isn't fulfilled.

Also, the whole point of doing OFDM is that you don't need an overall equalizer. So don't use one. OFDM in itself is designed to transform your equalizer-necessitating wideband channel into narrowband channels that are "flat" enough. So, whenever you apply an equalizer to an OFDM system: you're either doing something wrong, or the OFDM system was designed incorrectly (i.e. not for the channel that you're having).

(That is not to say that there's absolutely no room for time-domain equalization in OFDM systems, but it's going to be extremely specific corner cases.)


Why Is Bi Quadratic Interpolation for Image Resampling / Interpolation Rarely Done? Related question: What are the practically relevant differences between various image resampling methods?

Bilinear and bicubic interpolation for image resampling seem to be fairly common, but biquadratic is in my experience rarely heard of. To be sure, it's available in some programs and libraries, but generally it doesn't seem to be popular. This is further evidenced by there being Wikipedia articles for bilinear and bicubic interpolation, but none for biquadratic. Why is this the case?

You want to have an interpolation with a equal number of weights on both sides of the points you want to interpolate in between. So you choose either one or two weights on each side, resulting in an interpolation of two (linear) or four (cubic) points.

An quadratic interpolation would need three points, which would only make sense at the border of a grid, where you can choose only one point on one side and two on the other.


Adding noise to real MRI data I designed an algorithm for MRI data denoising which has good properties under heavy Rician noise (sigma is greater than 80). The method was tested on some phantoms and Rician noise generators. Now I want to test the algorithm on real MRI data. In the MRI data, which I have, the noise level is low. Is it methodologically correct adding noise to a real MR image to increase its level based on Rician noise generator?

Using real data is a way to validate the assumption that your algorithm will behave with real data in a similar way to the behavior with generated data.
In a general sense.
When we want to test an algorithm we can create tests, that are other algorithms that provide inputs and evaluate outputs of the algorithm we are design (your denoising algorithm).
Sometimes acquiring real data is difficult, but we can, making some assumptions generate data with similar distribution from data generator (the procedure by which you apply noise using your Rician noise generators), this makes easier to test your algorithm in much more data, and get more confidence about the algorithm.
When you write your data generator you are making one assumption about the types of distortion the images may have. Testing with real data you may find effects that are not documented, not studied, not expected. This is why in the real world the performance of algorithm may be degraded with measured data compared to generated data. Testing with real data will give you one idea of what is the impact of the neglected effects.
Different neglected effects may occur. The response may be non-linear and the  noise do not follow Rice distribution. The parameters of the noise distribution may depend on the position differently from what you assumed (you may have assumed that the distribution of the noise is independent of the position). Your scanner(s) may have different responses due to different temperature. Different scanners (same of different models/manufacturer) may produce data with different characteristics.
When you add noise to your images you are testing how the algorithm performs with the same real images and different Rician additional noise. There is nothing wrong in doing so, but you must know that worse quality images may have more of the distortions you are not considering.


Why is there a mismatch between normalized cross-correlation and perception for basic waveforms? To my ears, sawtooth and square waves are more similar to each other than either one is to a triangle wave. However, when I compare their magnitude spectra against each other using normalized cross-correlation, the square and triangle waves are the most similar to each other, followed by the sawtooth and triangle waves. Does anyone know what the source of this apparent mismatch is?

Sawtooth:

Even and odd harmonics, falling off at 1/n

Square:

Odd harmonics only, falling off at 1/n

Triangle:

Odd harmonics only, falling off at 1/n2

So while triangle and square have the same set of frequencies present, the spectral envelope of sawtooth and square is more similar, with much stronger high frequency harmonics.  Triangle waves are simpler/closer to a sine wave.

(Also, if you're generating them digitally, it's uncommon for them to be bandlimited, which produces worse aliasing artifacts with the square and sawtooth than with the triangle.  This would be most audible with frequency sweeps or high fundamental frequencies.)


Difference Between Convolution and Multiplication of Frequency Response and Frequency Spectrum Suppose I have  impulse response like [1/25,2/25,3/25,4/25,5/25,4/25,3/25,2/25,1/25].
I did convolution with 600 samples of test signal (it seems that I did some filtering).
Then I calculated the same thing in frequency domain:

Found frequency spectrum of test signal and frequency response.
Multiplied them and via IDFT convert back to time domain.

Then I plotted difference between these two signals and I got that difference between
these to signals is about |0.4|.

Is it ok? It seems a little bit huge for me -- the difference is about one-tenth.

function generateConvolution2(x, impResp)
{
    var xVal = x();
    var irVal = impResp();

    var convolution = [];
    var data = [];

    for (var i = 0; i < xVal.length + irVal.length - 1; i++)
    {
        convolution[i] = 0;
        for (var j = 0; j < irVal.length; j++)
        {
            if (i - j < 0) convolution[i] += 0;
            else convolution[i] += irVal[j] * xVal[i - j];
        }

        data.push(convolution[i]);
    }
    return data;
}

//calculate dft magnitude and spectrum of a function
function magnitudeAndPhase(length, dftLength, func)
{
    var data = [];
    for (var j = 0; j <= dftLength; j++)
    {
        var f = 0;
        var re = 0;
        var im = 0;

        for (var i = 0; i < length; i++)
        {
            f = func(i);
            re += f * Math.cos(2 * Math.PI * j * i / length);
            im -= f * Math.sin(2 * Math.PI * j * i / length);
        }

        //calculate magnitude
        var mag = Math.sqrt(re * re + im * im);

        //calculate phase
        var phase;
        if (re === 0)
        {
            re = 1e-20;
        }
        phase = Math.atan(im / re);
        if (re < 0 && im < 0) phase -= Math.PI;
        if (re < 0 && im >= 0) phase += Math.PI;
        data.push({mag: mag, phase: phase});
    }
    return data;
}

function multDFt(length, dftLength, func1, func2)
{
    var data1 = magnitudeAndPhase(length, dftLength, func1);
    var data2 = magnitudeAndPhase(length, dftLength, func2);
    var result = [];
    //1)multiply magnitudes
    //2)add phase
    for (var i = 0; i <= dftLength; i++)
    {
        result.push({mag: data1[i].mag * data2[i].mag, phase: data1[i].phase + data2[i].phase });
    }

    return result;
}

//using magnitude and phase calculate IDFT
function IDFT2(mainLength, length, dftLength, testSignal, imprResp)
{
    var rea = [];
    var ima = [];

    var dftMagResult = multDFt(length, dftLength, testSignal, imprResp);
    //create re and im part of coefficients
    for (var i = 0; i <= dftLength; i++)
    {
        rea.push(dftMagResult[i].mag * Math.cos(dftMagResult[i].phase));
        ima.push(dftMagResult[i].mag * Math.sin(dftMagResult[i].phase));
    }

    //prepare re and im coefficients
    rea[0] /= length;
    rea[dftLength] /= length;
    for (var j = 1; j < dftLength; j++)
    {
        rea[j] /= (length / 2  );
        ima[j] /= -1 * (length / 2  );
    }
    var vala = [];
    for (var i = 0; i < mainLength; i++)
    {
        vala.push(0);
    }
    //calculate idft
    for (var j = 0; j <= dftLength; j++)
    {
        for (var i = 0; i < length; i++)
        {
            vala[i] += rea[j] * Math.cos(2 * Math.PI * j * i / length);
            vala[i] += ima[j] * Math.sin(2 * Math.PI * j * i / length);
        }
    }

    return vala;
}

UPDATE: Added my js code. Don't bother with my js style (I'm beginner) and some inconsistent naming. Every func. return array because it is needed for plotting libs.

Thanks in advance.

Few notes:

Pay attention to the "Factors" in the DFT. Some use $ \frac{1}{2 \pi} $ some use others, pay attention you normalize accordingly.
Multiplication of the DFT is equivalent of Circular Convolution.
You should pad the signals accordingly to get the "Classic Convolution".


Why are Gold codes and Kasami codes used instead of pure m-sequences? I'm trying to understand why Gold codes and Kasami codes are used instead of pure m-sequences, in direct-sequence spread-spectrum (DSSS) communication systems, to prevent interference between multiple transmitters.

If I understand correctly, Gold codes are defined as the XOR between two m-sequences with different polynomials of the same degree (e.g. one LFSRs based on $p_1(x)$, and another LFSR based on $p_2(x)$ both of degree 20; the LFSR outputs are then XOR'd together), and within one system, multiple transmitters use the same pairs of polynomials but with a different time shift. (Kasami is also an XOR so I guess it's similar.)

How is this better than, say, having a system where all transmitters use m-sequences but each of them uses a different polynomial? (In both cases, the autocorrelation values are higher than the cross-correlation values)

There are very few m-sequences of any given length with good cross-correlation properties.  Their autocorrelation properties are excellent, but the cross-correlation properties are variable.  For example, there are 18 m-sequences of period 127 but to have good cross-correlation properties, one must choose a set of no more than $\require{cancel}\cancel{\text{3}}$ 6 m-sequences out of the 18 available. As a consequence, when several sequences with good cross-correlation properties are needed, one needs to use Gold sequence sets  or Kasami sequence sets which have good cross-correlation properties and good (though not excellent) autocorrelation properties.

For more than you probably want to know, I refer you to the paper D.V. Sarwate and M.B. Pursley, "Cross-correlation properties of pseudorandom and related sequences," Proc. IEEE, vol.68, pp.593-619, May 1980. It includes a detailed discussion of the Gold sequences as well as the small sets and large sets of Kasami sequences.


Threshold for peak detection of noisy signal in frequency domain I have a noisy signal (Gaussian Noise) with a known SNR and known Noise variance.
This signal contains a number of Frequency components.
My ultimate objective is to detect the components of the signal.
What is the Threshold formula should I use to detect the peaks of the signal in frequency domain given a certain Pfa (probability of false alarm)?

The following plot is an example of my objective. I want a Threshold that separates my useful components from the noise components (Let's say between 0.6 and 1)

This is a pretty common problem in radar signal processing. Typically, we’ll use what’s called a CFAR detector (constant false alarm rate). This essentially boils down to a non-linear smoothing filter where to attempt to estimate the local average power of the signal (excluding the current cell under test (CUT)). Afterwards, a “bias factor” is applied; this is typically determined a priori. Anything above the CFAR level + bias factor is labeled as a detection, and everything else is ignored. I’m sure there is a wealth of publicly available information on the topic, and I’d highly recommend doing some googling, or even looking around more on this website.

EDIT:
Neyman-Pearson is indeed the standard technique.  There are loads of papers about this topic on IEEE Xplore, I highly recommend you check them out.  That "bias factor" I talked about earlier is the thing you'll want to calculate based on your PFA or PD curve, which is a function of SNR.


Understanding zero-phase filters I'm reading the book "Digital Image Processing" by Gonzalez and Woods and I'm wondering how their definition of zero-phase-shift filter is equivalent to the one given here.

"Digital Image Processing" gives the following definition:

The DFT of $f(u, v)$ in Polar coordinates is:

$$F(u, v) = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1} e^{-j 2\pi (\frac{ux}{M} + \frac{vy}{N})} = |F(u, v)|e^{j\phi(u, v)}$$

where $\phi(u, v) = \arctan\left(\frac{I(u, v)}{R(u, v)}\right)$ and $I, R$ is the imaginary and real part of $F(u,v)$. Then filtering means multiplying the filter $H(u, v)$ with $F(u, v)$ and then taking the inverse DFT i.e. $F^{-1}\left(H(u, v)R(u, v) + j H(u, v)I(u, v)\right)$.

"The phase angle is not altered by filtering in the manner just described because $H(u, v)$ cancels out when the ratio of the imaginary and real parts is formed. Filters that affect the real and imaginary parts equally, and thus have no effect on the phase, are appropriately called zero-phase-shift filters.".

If I understand it correctly, they are saying $\phi(u, v) = \arctan\left(\frac{\text{Im}(H(u,v))I(u, v)}{\text{Re}(H(u,v))R(u, v)}\right) = \arctan\left(\frac{I(u, v)}{R(u, v)}\right)$.

But the other definition looks different. They focus on the imaginary part of $H(u, v)$ being zero. This makes more sense to me, because then $\phi(u, v) = 0$. Right? But how does this affect "the real and imaginary parts equally"?

Filters that affect the real and imaginary parts equally, and thus
  have no effect on the phase, are appropriately called zero-phase-shift
  filters."

The filter itself is zero-phase, but the input (and therefore the output) waveform can be complex. In this case the real and imaginary components of the input would be effected equally since the filter is not adding any additional phase shift.


The Technology Behind Animoji by Apple Last year Apple releases their new phone Apple X that has an exciting feature called Animoji. I am so fascinated by this technology that I really want to make my own version of it. But the problem is I am really new in image processing stuff and I've no idea about the technology behind Animoji. So what is the technology behind the Applle Animoji and what sort of knowledge should I have if I want to build my own version of Animoji from scratch?

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


What is the significance of the DFT in SC-FDMA systems? Long term evolution (LTE) systems use single-carrier frequency division multiplex (SC-FDMA) in the uplink. In such a system, there is a DFT block and a sub carrier mapping block before the IFFT block. What is the significance of DFT and sub carrier mapping? And how do these blocks make the system single carrier?

In the LTE uplink, user equipments (UE) of different users share the same frequency band when sending data to the base station (BS). This frequency band is divided into smaller resource units which must be used by only one UE at a time. This technique is called frequency-division multiple access (FDMA).

Let's first assume that there are $M$ UEs and that every UE uses conventional orthogonal frequency division multiple access (OFDMA). That is, every UE has an IDFT block of length $N$. Different UEs must not use the same subcarriers at a time so their transmit data is distributed onto orthogonal sets of subcarriers.

  Example: Let $N=30$ and $M=3$. To implement FDMA we could assign
  subcarriers 1 to 10 to UE1, 11 to 20 to UE2 and 21 to 30 to UE3.

The drawback of this technique is that the resulting transmit signals have a relatively high peak-to-average power ratio (PAPR).

The prinicple of SC-FDMA is to generate (in every UE) a single-carrier signal and to shift this signal to the frequency resource that is assigned to the respective UE. The necessary frequency shift is implemented by the $N$-IDFT. The three principal blocks are:

DFT: Take $M$ QAM symbols $y(n)$ (which can be interpreted as $M$ consecutive symbols of a single-carrier signal) and calculate the frequency-domain representation $Y(k)$ of this signal by an $M$-DFT.
Subcarrier mapper: "Place" $Y(k)$ in the freq. domain at the freq. resource that is assigned to the UE. This is done by assigning $Y(k)$ to the  $M$ IDFT inputs that represent the wanted freq. resource and by setting the remaining $N-M$ IDFT inputs to zero.
IDFT: Calculate the time domain signal.

  Example: Again, $N=30$, $M=3$. In every UE, calculate the DFT of $N/M=10$ QAM transmit symbols. In UE1, assign the result to IDFT inputs 1 to 10. In UE2, assign the result to IDFT inputs 11 to 20. In UE3, assign the result to IDFT inputs 21 to 30. The remaining IDFT inputs are set to zero, respectively.

The DFT makes SC-FDMA a single-carrier signal because the input to the IDFT is a "spectrum" already not QAM symbols. Therefore, the result of the IDFT is a frequency-shifted single-carrier signal.

I'd like to point out, that the whole purpose of SC-FDMA is to reduce the PAPR. Especially, SC-FDMA is not necessary to implement FDMA. This can also be done by OFDMA.

Also, there are different subcarrier mapping schemes. In my answer I have assumed the so-called localized FDMA (LFDMA) which assigns $Y(k)$ to neighbouring IDFT inputs. Another technique is the interleaved FDMA (IFDMA) which assigns $Y(k)$ to IDFT inputs with distance $M$ in an interleaved way. It has an even lower PAPR. The PAPR of OFDMA, LFDMA and IFDMA are compared in [1].

[1] Myung, H., Lim, J., & Goodman, D. (2006). Peak-To-Average Power Ratio of Single Carrier FDMA Signals with Pulse Shaping. 2006 IEEE 17th International Symposium on Personal, Indoor and Mobile Radio Communications


This will give you a plot of the autocovariance up to lag 100 samples:

plot(autocov(a,100))

There you can clearly get the period of your signal.

Another approach is to explicitly get the time of each pulse:

pulses = a > 0.1;
leading_edges = diff(pulses) > 0;
times = find(leading_edges > 0);
periods = diff(times)

ans =
    55   56   56   56   56

Explanation: pulses is 1 wherever there is a pulse. leading_edges gets a 1 only on the first sample of a pulse. times gives the sample number of each leading edge.

With this result you can perform your own statistics or estimation. This second method is computationally more efficient than using transforms. Computing time is linear with number of samples.


Follow up question regarding: "Complex sampling can break Nyquist?" I'm having some trouble understanding the sample rate limitations when considering a complex baseband signal.

I understand (based on the linked SE questions below), the either (1) physically sampling a signal with frequencies between $0$ and $B$ using 2 ADCs at sample rate $B$, one 90 degrees out of phase, or (2) sampling with one at $2B$, can satisfy the Nyquist criteria.

I'm working with a system where we have a real digital signal of bandwidth $B$, sampled at $2B$. It has then been converted to complex baseband representation by calculating the I/Q components from the existing samples. After that, the center frequency (also carrier in this example) is moved to zero and the signal now occupies $-B/2$ to $B/2$.

I suspect that even though the complex representation is twice the data, I still cannot down sample without breaking the Nyquist criteria because my I/Q components are dependent. For example, each real sample in the original signal is used to compute an I and Q component for that time.

Is this actually the case? Or can I somehow down sample by 2 and still fully represent the signal in complex form with a sample rate of $B$? To clarify, I believe that I could if the I/Q components were physically sampled independently.

I've done a few simple tests of this in python using a chirped signal, I could share some of the code if it will help. So far, if down sampling does violate Nyquist, I'm still able to reconstruct the original signal nearly as well as in the case with no down-sampling.

This question: Complex sampling can break Nyquist?

This answer on another question: Answer

To represent a complex baseband signal in the range $[-B/2,B/2]$ you need either $2B$ real samples or $B$ complex samples per second.

To take real samples, you write the signal as $I(t)+jQ(t)$, with $I(t)$ and $Q(t)$ real, and sample each one at $B$ real samples per second, for a total of $2B$ real samples per second. It should be clear that you can't downsample without aliasing.

You can also take complex samples at $B$ samples per second (which is just a different way of thinking about the real sampling described above). This allows you to represent the entire band $[-B/2,B/2]$. Downsampling in any way will also result in aliasing.

To be more specific regarding your application: You have a bandpass signal with bandwidth $B$ and take real samples at rate $2B$. When you calculate the complex envelope, the I and Q components can be sampled at a rate of $B$ real samples per second each. The complex envelope can also be sampled at a rate of $B$ complex samples per second. In any case, you end up with a total of $2B$ real samples per second.


How to get the autocorrelation of this finite discrete function? This is the function:

$$
x(n) =
\begin{cases}
1, & \text{-N  $\leq$ n $\leq$ N ;}\\
0, & \text{otherwise.}
\end{cases}
$$

And I tried to solve it, using sigma notation or summation:

\begin{equation}
\ r_{xx}(l) = \displaystyle\sum_{n=-\infty}^{\infty} x(n)x(n-l) \end{equation}

I get 2N-1 but that is wrong. Because when I do the autocorrelation in MATLAB I get something like:
For exmple if N = 2, the discrete function rxx(l) would be like:

Which is something like rxx(l) = 2N-1+|l|
But I don't know how to get there... Can you help me?

Calculating the autocorrelation of this function is relatively easy. You just need to follow the formula you have written:

\begin{equation}
\ r_{xx}(l) = \displaystyle\sum_{n=-\infty}^{\infty} x(n)x(n-l) \end{equation}

It gets much easier because your function is $\neq0$ for only a few values ($2N+1$ exactly).

Because of this the previous sum can be replaced with this much simpler one:

\begin{equation}
\ r_{xx}(l) = \displaystyle\sum_{n=-N}^{N} x(n-l) \end{equation}

Which is basically just counting how many times $x(n-l)$ is not zero on that finite interval. Because of this it appears easily that

\begin{equation}
r_{xx}(l)=max(2N-1-|l|,0)  \end{equation}

You almost had it in your question but the sign was wrong.


The Meaning of the Terms Isotropic and Anisotropic in the Total Variation Framework The isotropic TV is defined as the estimation of 2-norm of gradients $\sqrt{(y_{i+1,j}-y_{i,j})^2+(y_{i,j+1}-y_{i,j})^2}$, while the anisotropic TV is defined as the estimation of 1-norm of gradients $|y_{i+1,j}-y_{i,j}|+|y_{i,j+1}-y_{i,j}|$.

Now I am wondering why the second one will be called anisotropic. Since from my perspective, the isotropic TV bounded the 2-norm of the gradient, and thus it can bound the norm of gradients for every directions. But for the second one, since the  absolute values are bounded, it means that the 2-norm of every directions can also be bounded. Isn't it also an isotropic way to bound the gradient?

In the Total Variation framework we define 2 flavors:

$$ \text{Isotropic TV} \; {TV}_{ {L}_{2} } \left( X \right) = \sum_{ij} \sqrt{ { \left( {D}_{h} X \right) }_{ij}^{2} + { \left( {D}_{v} X \right) }_{ij}^{2} } $$

$$ \text{Anisotropic TV} \; {TV}_{ {L}_{1} } \left( X \right) = \sum_{ij} \sqrt{ { \left( {D}_{h} X \right) }_{ij}^{2} } + \sqrt{{ \left( {D}_{v} X \right) }_{ij}^{2} } $$

Where $ {D}_{h} $ and $ {D}_{v} $ are the horizontal / vertical finite differences operators.

Quoting Isotropic and Anisotropic Total Variation Regularization in Electrical Impedance Tomography:

  It is worth noticing that referring to functional as ”anisotropic” is slightly misleading, because the weights corresponding to all coordinate directions are equal.

The term Anisotropic is used because this flavor separate Vertical and Horizontal components as done in Anisotropic Diffusion. Hence it doesn't have rotational invariant term as the Isotropic Flavor. As you know, sum of quadratic terms is connected to Ellipse and with equal weights it is a Sphere / Circle which is Isotropic.

Great reading about the subject is Laurent Condat - Discrete Total Variation - New Definition and Minimization (DOI Address). It analyze the properties of different TV definitions regarding their spatial smoothing.


Difference Between Contrast Enhancement and Texture Enhancement in Images I have a very basic question. I understand that contrast enhancement is exaggerating the difference between the intensity of the adjacent structures and texture in an image is the repeating geometric patterns in an image. But what is exactly texture enhancement? Does enhancing contrast in an image automatically enhance texture also? Can someone explain this with an example. By looking at an image, how can one identify if the image needs contrast enhancement or texture enhancement? Or how does one identify, if an algorithm has done contrast enhancement or texture enhancement?

Usually Contrast Enhancement is Element / Pixel Wise Operation in Digital Image Processing.

It means it adjusts the value of a pixel based only in its value.
In order to increase contrast we usually push bright values to be brighter and dark values to be darker (S Shaped Contrast).

When people try to enhance Texture they usually refer to the term Local Contrast.
Local contrast means we want to enhance contrast but using a Neighborhood Operation.
In the sense of contrast, we will push a pixel to be darker even if its value is bright but it is darker than its neighbors. Namely we have the locality sense of the operation.
For instance High Pass will increase local contrast.


Sum of correlated exponential random variables How do we find PDF of sum of correlated exponential random variables. I know for independent random variables. But how to find it for correlated exponential random variables.

The pdf $f_Z(z)$ of the sum $Z=X+Y$ of any two jointly continuous random variables $X$ and $Y$ with joint pdf $f_{X,Y}(x,y)$ is as follows:
$$\text{For all } z, -\infty < z < \infty, ~~ f_Z(z) = \int_{-\infty}^\infty f_{X,Y}(x,z-x) \, \mathrm dx.\tag{1}$$

For the special case when $X$ and $Y$ are nonnegative random variables (including as a special case, exponential random variables) and so take on nonnegative values only, $f_{X,Y}(x,y)$ has value $0$ if at least one of $x$ and $y$ is smaller than $0$. Hence, in this case, the integrand $f_{X,Y}(x,z-x)$ in $(1)$ has value $0$ if $x < 0$ or if $z < x$. Consequently, if $z$ is a negative number, then the integrand in $(1)$ is always $0$ regardless of the value of $x$ and therefore so is the integral. All of which is just a long-winded way of saying that $f_Z(z)$ has value $0$ when $z<0$, that is, $Z$ takes on nonegative values only, which any idiot could have deduced from the fact that $Z=X+Y$ and both $X$ and $Y$ are nonnegative. But the approach is useful even for $z>0$ since now we have that the integrand in $(1)$ is zero when $x<0$ or when $x >z$ and so for nonnegative $X$ and $Y$, we can simplify $(1)$ to
$$f_Z(z) = \begin{cases}\displaystyle\int_0^z f_{X,Y}(x,z-x) \, \mathrm dx, & z \geq 0,\\\quad\\
0, & z < 0\end{cases} \tag{2}$$

  No further simplification of $(2)$ is possible in general.

For the special case when $X$ and $Y$ are independent random variables, $f_{X,Y}(x,y)$ factors into $f_X(x)f_Y(y)$ and so $(1)$ becomes the familiar convolution integral and $(2)$ the somewhat-less-familiar convolution integral for causal signals. But no such simplification is possible for nonindependent random variables $X$ and $Y$; we need the joint pdf to calculate $f_Z(z)$ and just knowing that $X$ and $Y$ are correlated random variables (whether exponential or Gaussian or whatever) is not enough.


SC-FDMA disadvantages I understand the advantages of SC-FDMA over OFDM and I also understand the disadvantage of additional complexity (FFT block), are there another disadvantages?

Also, if $N$ has to be bigger than $M$ ($N>M$) then I don't understand why I don't lose rate as a result of padding with zeros instead of transmitting another data subcarriers.

Some more disadvantages:

Channel estimation using pilots is harder, because you do not have orthogonal data on each frequency bin, compared to OFDM. That's why e.g. LTE uplink resorts to a midamble for channel estimation in uplink
Application of non-linear MIMO detection algorithms is harder (i.e. more complex), because you essentially have coupled MxM elements from each transmit antenna, compared to 1x1 elements when using OFDM.
In high SNR regions I think SC-FDMA does not perform as good as OFDM when considering a coded transmission with linear detectors.

Regarding the data loss: Even if $M<N$, there is no real loss, because the used bandwidth of your signal is only M frequency bins (ignoring any guard carriers). So, any other system can use the remaining N-M frequency bins.


System classification: unit-time delay I'm reading a book on linear systems and I can't understand why the unit-time delay is a distributed system. This is the example given in the book:

I understand that the initial state of the system is the input, that is:

$u(t),\,t_0-1\le t<t_0$

However, for me, a distributed system is something like a transmission line, where the parameters can't be modeled accurately by "discrete" and finite elements.

I can't grasp why this system is classified as distributed just because there are an infinite number of points describing the initial state.

You are right that a distributed system could be "something like a transmission line". Note that the system

$$y(t)=x(t-T)\tag{1}$$

is a simple model of a transmission line, where just a frequency-independent delay $T$ is taken into account, and the attenuation is neglected.

Note that lumped electrical systems, described by resistors, capacitors and inductors, result in ordinary differential equations, i.e., the output $y(t)$ depends on the time derivatives of the input and output signals. If the system is also time-invariant, i.e., the component properties do not change over time, you can define a transfer function, and with lumped elements that transfer function is always rational. The transfer function describing the system given in $(1)$ is

$$H(s)=e^{-sT}\tag{2}$$

which is clearly non-rational.


How to remove reflections from images? I've been playing around with making images better and removing all the unwanted disturbances from it. I found out a major concern would be removing the unwanted reflections that might come.

I've gone through a paper that uses ICA to get rid of these reflections. What other methods would work?

Say I have an image:

Another image:

Converted it into CMYK

How should i remove this semi black portion?

I am not a specialist in image processing but I highly doubt this is feasible with your requirements (single image, no prior information). The challenges is twofold:

Detect which area of the image is a reflection. This looks like a quite complex scene understanding problem. In your example, I know that the right side is a reflection because I can infer from the image that a person is depicted here, and I have prior knowledge about what the back of a person looks like.
Reconstruct the background in the identified region. While it would be possible to just inpaint the area detected as a reflection with textures from the surrounding, the result would hide many details! Whatever I've seen in the field of single input image separation was relying on very simple image models (stationary textures), not real world images. The more traditional image separation setups (multiple channels, recto/verso of translucent pages, multiple polarization) do not apply to your case.


Obtaining phase and magnitude of 2 spectral components after ADC I am currently implementing a digital filter to obtain the output of a PID controller. The control loop is carried out by applying a 1kHz sine wave at the input of the system and reading the output signal. The signal is distorted by the system and the output signal contains a 2nd harmonic component of the 1kHz input sine wave. I am currently applying FIR filters to obtain the phase and magnitude of the fundamental and harmonic tones but in order to reduce the noise I have to use very high order FIRs (around 1000 taps) and it is a very slow process. I would like to know if a matched filter would be a better approach, or if there are any other suggestions. Thanks!

To obtain the phase and magnitude of a known frequency, I recommend multiplying each sample by $exp^{-j\omega t}$ where $\omega$ is $2\pi 1000$ Hz for the case of the fundamental (and $2\pi n 1000$ for each harmonic n).  The output can then be "accumulated and dumped" over a number of samples depending on how much noise is present, traded with possible offset errors in the actual frequency, and how stationary the process is. For every N samples accumulated, the noise if white will go down relative to the accumulated output by $\sqrt{N}$.  If this measurement is being used in the control loop itself then there is the additional complexity of ensuring that the accumulation time is significantly shorter than your loop response time or this process will become part of your loop filter, which you probably don't want!  (In this case I would likely let the loop filter take care of all the filtering.)

This "accumulate and dump" is equivalent to a moving average filter followed by a decimator, so a 1 ms accumulation time would be very attractive as the nulls would all be spaced by 1 KHz in that case (conveniently nulling out all the harmonics for all cases)!   Further this "correlator" which is what this is would have an amplitude vs frequency offset sensitivity given by the Sinc frequency response of the moving average filter.

Then with the accumulate and dump output (I and Q using the complex method I described) phase, frequency and magnitude estimates can be easily obtained with minimum processing which can additionally be further averaged to increase measurement SNR (if the system is sufficiently stationary), or include in a correction loop which I suspect the OP may be trying to do.

Here is a common magnitude algorithms for efficiently estimating magnitude from I and Q samples (if one doesn't want to or can't calculate $\sqrt{I^2+Q^2}$ directly): http://dspguru.com/dsp/tricks/magnitude-estimator/.

The phase information is the phase of the I,Q sample since we have translated the signal of interest to baseband with a complex rotation (so in that case we need to compute $ATAN2(Q,I)$ which can be done with a look-up table or the Cordic Rotator if an iterative approach is allowed.

NOTE: If the intention was to implement a phase lock loop, then this becomes even simpler in that a multiplier can be used as the phase detector (as the output of a multiplier would be proportional to the cosine of the phase between the two input signals), knowing that the loop will drive the phase detector to it's linear region at the 90° crossing, and in this case a real 1KHz signal can be multiplied by the test signal.

For efficient frequency estimation please see the figure below that shows how a frequency estimate can be extracted from this simple operation on two consecutive complex samples:

$$Freq \varpropto I[n]Q[n-1]-I[n-1]Q[n]$$

This comes about because the complex conjugate multiplication of any two vectors will be a vector with a phase that is equal to the phase between the two vectors (so for small angles and normalized vectors, the Q of that vector is the phase!).    This is known as a cross-product discriminator for estimating phase (and in the case when the two vectors are the same signal delayed in time, this would be change in phase over a change in time which is frequency).  Any such frequency measured would be the frequency offset from the 1KHz, or whichever harmonic was selected.

This plot below shows how the phase is derived (and here the vectors $V_1$ and $V_2$ do not relate to the same vectors in the previous plot).   The dashed line shows the conjugate of $V_1$. When two vectors are multiplied, their magnitudes multiply and their angles add, resulting in the product vector drawn in red that I will call $V_3$. If $V_1$ and $V_2$ were normalized to be 1, and if the angle between $V_1$ and $V_2$ was small, then the imaginary portion of $V_3$ (Q) would be equal to the angle with no arctan calculation necessary!


Confusion regarding pdf of circularly symmetric complex gaussian rv Considering a random variable $x$ that takes in values from a complex domain. Its real and imaginary components are totally uncorrelated. I am following this link and
also studying this document.

In the wikipedia link, the pdf for a single observation is given and $k =2$ for bivariate gaussian as my assumption is that the real and imaginary are totally uncorrelated and are gaussian respectively. There is a 2 in the denominator in the pdf but the log likelihood for the complex case does not have any 2 in the denominator. The log-likelihood for $k=2$ would be $$-N \ln| \Sigma| - (x_n) {|\Sigma|}^{-1}(x_n)^H - k N \ln \pi$$
where $\mu=0$ as I am assuming zero mean r.v $x$

Confusion 1:

I was thinking that the variance of $x$ was
$$\Sigma = \begin{bmatrix}\sigma_1^2 & 0\\ 0& \sigma_2^2\end{bmatrix}\,\text.$$
So for $N$ samples (observation), the joint pdf would turn out to be
$$P_x(x_1,x_2,...,x_N) = \prod_{n=1}^N\frac{1}{\pi \sigma^2_x} \exp \bigg(\frac{-{({x_n})}^H ({x_n})}{\sigma^2_x} \bigg)$$
where $\sigma^2_x = [\sigma_1^2,\sigma_2^2]$. Would there be a 2 in the denominator of the pdf and what is the correct pdf expression?
Confusion 2:

In the document, the expression for the pdf for complex case looks different from the wikipedia link. Are they the same but maybe I am missing some link between them? Which pdf should I use?
Confusion 3:

When simulating the r.v in Matlab with zero mean and variance 1, I am halving the variance because the real and imaginary part's variance should add up to 1. Then, would the pdf also contain half of the variance as $\sigma^2_x/2$?

Let me try to establish the relation between the univariate PDF for a real Gaussian and the univariate PDF for a complex proper (i.e. circular symmetric) Gaussian.

You know that
$p_x(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp(\frac{-x^2}{2\sigma^2})$
is the PDF of a real-valued Gaussian with variance $\sigma^2$. We write $x\sim\mathcal{N}(0,\sigma^2)$ to denote that $x$ is a random variable that follows a real-valued Gaussian with zero mean and variance $\sigma^2$.

Now, let $z=a+jb$ be a circular symmetric random variable with real part $a$ and imaginary part $b$. In a circularly symmetric Gaussian random variable, the real and imaginary part are i.i.d., i.e. $a\sim\mathcal{N}(0,\sigma^2)$ and $b\sim\mathcal{N}(0,\sigma^2)$. Since $a,b$ are independent, their joint PDF is the product of their PDFs.

$$p_z(z=a+jb)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(\frac{-a^2}{2\sigma^2}\right)\cdot\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(\frac{-b^2}{2\sigma^2}\right)$$

which gives the PDF of the complex variable $z$ at the value $a+jb$. We can reformulate this to
$$\begin{align}p_z(z=a+jb)&=\frac{1}{2\pi\sigma^2}\exp\left(\frac{-(a^2+b^2)}{2\sigma^2}\right)\\&=\frac{1}{2\pi\sigma^2}\exp\left(\frac{-z^*z}{2\sigma^2}\right)\\&=\frac{1}{2\pi\sigma^2}\exp\left(\frac{-\|z\|^2}{2\sigma^2}\right)\end{align}$$

We also write for this case that $z$ follows a circular Gaussian distribution and denote this by $z\sim\mathcal{CN}(0,2\sigma^2)$. Why is it suddenly $2\sigma^2$ (i.e. the double variance compared to the uni-variate, real case? Well, because $E[\|z\|^2]=2\sigma^2$ since it consists of a real and an imaginary part (which are independent and hence their variance adds up together).

So, to sum up:

if $x\sim\mathcal{N}(0,\sigma^2)$, then $p_x(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{x^2}{2\sigma^2}\right)$
if $z\sim\mathcal{CN}(0,\sigma^2)$, then $p_z(z)=\frac{1}{\pi\sigma^2}\exp\left(-\frac{z^*z}{\sigma^2}\right)$

Does this answer your confusions?


Tool for Harmonic Analysis & Synthesis I've been using a program for identifying harmonics for a wave after doing a fourier analysis.

http://www.phon.ucl.ac.uk/resource/sfs/esynth.htm

I would like to be able to tweak some of the harmonics involed in a wave and also to add harmonics.

The program mentioned above permits this kind of operations but it isn't precise.

Does anyone know of a good tool for harmonic analysis and synthesis like the one mentioned above?

ps:It's important that the program exports the harmonics to a file.

thanks.

Try NCH Software

From Wikipedia:

  NCH Software is a software development company founded in 1993 in
  Canberra, Australia. The Colorado office was started in April 2008 due
  to the large U.S. customer base. NCH Software primarily sells to
  individuals via their website.

And further

  Browser Toolbar

  During the installation process NCH Software also by default (but you
  can uncheck the option during install) installs Google Chrome. This
  can be unselected during install or removed after install. In older
  versions of their software a browser toolbar powered by Conduit,[11]
  was installed as well as making changes to the default homepage and
  search engine.


Machine learning for denoising MRI images I'm currently a sophomore in college and pretty new to the field of research. I'm currently working on an existing algorithm for MRI de-noising and the results are nothing that great. I can't see how I can improve the results.

I've tried a few different tweaks and played with the code and parameters for sometime now. I also saw some papers on using neural networks for image de-noising as well.

I'm thinking of trying the same but I don't know if it makes sense to use neural networks on an already existing algorithm. Won't the neural network overpower my algorithm and all the results be based on it only? So if I use a pre-existing MatLab filter and use neural networks in it, will it differ much if I first use my algorithm and then use neural networks in it?

I have little experience with machine learning and hope that I made sense in my question. Any kind of help will be appreciated. Thanks.

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

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

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


Frequency Support of Window Functions All window functions, which go from zero amplitude, up, then back to zero amplitude, are all base-band signals. Why is this? Are all functions that have this behavior always band-band aka low-pass?

A window function $w(t)$ satisfies $w(t)\ge 0$. The value of a window's frequency response $W(f)$ at DC ($f=0$) equals its integral
$$W(0)=\int_{-\infty}^{\infty}w(t)dt>0\tag{1}$$
which is clearly greater than zero because $w(t)\ge 0$. For all other frequencies we obtain the following bound:
$$|W(f)|=\left|\int_{-\infty}^{\infty}w(t)e^{-j2\pi ft}dt\right|\le\int_{-\infty}^{\infty}\left|w(t)e^{-j2\pi ft}\right|dt=\int_{-\infty}^{\infty}w(t)dt=W(0)\tag{2}$$
Eq. $(2)$ means that $|W(f)|$ cannot be greater than its DC value. Actually, it can be shown that $(2)$ is a strict inequality for any $f\neq 0$: $|W(f)|<W(0)$ for any $f>0$. This shows that $W(f)$ has a lowpass characteristic.
To give a more intuitive explanation, a rectangular window's frequency response is a sinc function, which has a lowpass characteristic. All other windows are just smoother versions of a rectangular window, which means that their spectra decay faster than a sinc. I.e., non-rectangular windows have an even more pronounced low pass characteristic than a rectangular window.


Books/resources for implementing various mathematical functions in fixed point arithmetic for DSP purposes I'm looking for books or resources that cover the following in detail:

implementing mathematical functions (e.g., logarithm, exponential, sine, cosine, inverse) in fixed point arithmetic for DSP purposes.
techniques like using lookup tables, Taylor series, etc.

I am fairly familiar with C programming and am more interested in the algorithms on how to go about implementing various mathematical functions in an efficient manner.


Calculating the phase shift between two signals based on samples I am comparing two signals in MATLAB Simulink for finding the phase between them. To do this I am inspired by using the code found here.

I have two vectors of the same size which are a collection of samples of the two signals (sampling is more than fast enough). They are sine-signals with mostly the same frequency. To find the phase between them the code below is used

signal_one; %A vector with values of signal 1
signal_two; %A vector with values of signal 2

dot_product = dot(signal_one,signal_two);
norm_product = (norm(signal_one)*norm(signal_two);
phase_shift_in_radians = acos(dot_product/norm_product)

This code works satisfactory, and I am able to get the phaseshift from it, I just don't understand why. Can anyone show me why this gives the phase shift? A mathematical demonstration or a reference to an equation would be great.

It's an approximation as also admitted by your statement "satisfactory performance". Here is why.

Let $a[n]=\cos(w_0 n)$ and $b[n]=\cos(w_0 n + \theta)$ be two sequences of size $N$, represented by $1 \times N$ row matrices $\bar{a}$ and $\bar{b}$ in a program (such as Octave or Matlab)

Treating those two sequences $a[n]$ and $b[n]$ which are represented by two row matrices $\bar{a}$ and $\bar{b}$ as components of vectors $\vec{a}$ and $\vec{b}$, we can take advantage of the dot product stated for vectors and (approximately) compute the phase angle $\theta$ in between those two sinusoidal sequences $a[n]$ and $b[n]$

Geometric evaluation of the dot product between two vectors of same size is:
$$ \vec{a} \cdot \vec{b} = |\vec{a}||\vec{b}| \cos(\phi)$$

And the angle between those two vector therefore is:
$$\cos(\phi) = \frac{\vec{a}\cdot \vec{b}}{|\vec{a}||\vec{b}|}$$
Where $\phi$ is the angle between those vectors,  which will be equal to $\theta$ as shown:

The dot product between the vectors $\vec{a}$ and  $\vec{b}$ can also be algebrically computed from their MAC sum of components which are contained in the elements of the matrices  $\bar{a}$ and $\bar{b}$:
$$ \vec{a} \cdot \vec{b} = \sum_{i=1}^{N} \bar{a}(i)\bar{b}(i) $$

The product $\bar{a}(i)\bar{b}(i)$ can be shown to be equal to the following (using trigonometry):
$$\bar{a}(i)\bar{b}(i) = \cos(w_0 i) \cos(w_0 i + \theta)= 0.5\cos(2 w_0 i + \theta) + 0.5 \cos(\theta)$$

And the dot-product summation becomes:
$$ \vec{a} \cdot \vec{a} = \sum_{i=1}^{N} a(i)b(i)  =  0.5 N \cos(\theta) + \sum_{i=1}^{N} 0.5 \cos(2 w_0 i + \theta)$$

The absolute values $|\vec{a}|$ and $|\vec{b}|$ of the vectors $\vec{a}$ and $\vec{b}$ can also be computed from their norms computed over the matrices $\bar{a}$ and $\bar{b}$ as follows:

$$ |\vec{a}| = ||\bar{a}||=  ( \sum_{i=1}^{N} a(i)^2 )^{0.5}$$
$$ |\vec{b}| = ||\bar{b}||=  ( \sum_{i=1}^{N} b(i)^2 )^{0.5}$$

Again expanding the squared terms result in:
$$ |\vec{a}| = ||\bar{a}||=  ( 0.5 \sum_{i=1}^{N} 1 + \cos(2 w_0 i)  )^{0.5}$$
$$ |\vec{b}| = ||\bar{b}||=  ( 0.5 \sum_{i=1}^{N} 1 + \cos(2 w_0 i + \theta) )^{0.5}$$

Now, it can easily be shown that the following summations are small compared to large $N$:
$$ N+\sum_{i=1}^{N} \cos(2 w_0 i) \approx N , ~~~ N+\sum_{i=1}^{N} \cos(2 w_0 i + \theta) \approx N  $$
for large $N$ and suitable $w_0$:

Simplify the above equations to get the angle:
$$\cos(\phi) = \frac{\vec{a}\cdot \vec{b}}{|\vec{a}||\vec{b}|}$$
$$\cos(\phi) \approx \frac{ 0.5 N \cos(\theta)}{|\vec{a}||\vec{b}|}$$
$$\cos(\phi) \approx \frac{ 0.5 N \cos(\theta)}{0.5 N}$$
$$\cos(\phi) \approx \cos(\theta)$$
$$\phi \approx \theta$$

The approximation depends on the length $N$.


filter bank of gabor filters with different frequencies and orientations I was reading a paper which uses a filter bank of gabor filters with different frequencies and orientations. However, I often see the construction of a filter bank of gabor filters with different frequencies and scales.

I am not sure on the construction of it. How do I construct it, any generic example will be ok. Thanks!

Look on wikipedia. https://en.wikipedia.org/wiki/Gabor_filter The $\theta$ controls the orientation.


why is my FFT magnitude of real time recording so volatile? My goal is to do spectrum analysis of a gear motor. I used Python with the Pyaudio package and made a crude spectrum analyzer that displays the FFT as the sound is being recorded. I recorded 10 seconds of sound just for test. Sample rate is 44100, each frame is 1024 samples. For each frame I plotted the FFT, and made a dynamic plot thru the 10 seconds.

I also downloaded a FFT analyzer to my smartphone and compared the results as I turned on a motor. The smartphone app is able to display relatively stable spectrum. In my plot, the magnitudes are very high upon turning on the motor, and only "settled" after a couple of seconds.

I wonder why that is. Is it due to my sound card or laptop mic? I tested using another laptop and got similar result. Is it because I recorded in Mono mode? The spectrum changes too much as I move the motor around, whereas on the smartphone the spectrum is more stable.

Or is it because I didn't use a Window function? Here is my code:

import numpy as np
import os
import time
import pyaudio
import matplotlib.pyplot as plt
from scipy.fftpack import fft

CHUNK = 1024
FORMAT = pyaudio.paInt16
CHANNELS = 1
RATE = 44100
RECORD_SECONDS = 10
WAVE_OUTPUT_FILENAME = "test.wav"

p = pyaudio.PyAudio()

stream = p.open(format=FORMAT,
                channels=CHANNELS,
                rate=RATE,
                input=True,
                frames_per_buffer=CHUNK)

input('Press key to start recording:\n')
print('1 sec delay started') #delay so keystroke doesn't get recorded
time.sleep(1)
print("* recording")

frames = []
n = 1024
k=np.arange(n)
T = n/RATE
frq = k/T
frq = frq[range(int(n/2))] # one side frequency range
for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
    data = stream.read(CHUNK)
    frames.append(data)
    decoded = np.fromstring(data, dtype=np.int16) #grab the data in stream
    fft_decode=fft(decoded)/(len(decoded)/2) #normalized FFT
    mags=np.absolute(fft(decoded)) #
    plt.ylim(top=55000)
    plt.xlabel('Freq (Hz)')
    plt.ylabel('|Y(freq)|')
    plt.plot(frq, mags[range(int(n/2))],'b')
    plt.pause(.01)
    plt.gcf().clear()

print("* done recording")

plt.close()

How can I improve my analyzer?
 Is there anything wrong with my code or the way I approached this?

Edit: Nevermind. The reason turns out to be the boost function of my microphone, and disabling it solved my problem.

Are you appending new data to your FFT input vector?  If so, then every FFT result will have a different length; and it's likely the early ones are too short for the frequency resolution you desire.  The smartphone app may also be averaging multiple successive FFT results to reduce noise.


Sampling bandpass signal with spectral replications I want to clear my thinking on how to approach problems of the type characterized below.

I guess what we want to do is select a sample rate of a sufficient size such that this signal can be fully captured  therein but also small enough so as not to waste all the space from 0 to, 156? (77*2).

I've seen in some examples that if we select a sampling rate that is some integer multiple of this rate then we can use one that is much lower than the nyquist frequency, taking advantage of the property of replication.

So maybe for part (a), this would just have the same thing repeated over and over every 56 MHz? and How about specral inversion, would every other sample be inverted?

For part (b) maybe it's 56/2 = 28? or maybe it's 77/2?

For part (c) are we meant to think about the

fs' < 2fc-B/m

or maybe the

2fc-B/m > fs > 2fc+B/m+1

I guess m in this situation is 56? And B is 14? But how about the transition of B, it doesn't drop down straight but it flares out at the bottom, how to think about that? And I guess fc = 70?

EDIT: (reaction to Matt L.'s answer)

Is the image below what you had in mind for part (a)?

But why don't we take the bandwidth into account? Furthermore, in the image it seems that the 14MHz indicated there as the bandwidth doesnt completely describe the size of B since it seems to flare out at the ends, how about that?

So can we generalize that equation to be that for any time we are performing that undersampling we can get the location of all the replications by:

'Highest frequency component for which we must account' - ('number of K term' * 'sampling rate')

and sampling rate must be > 2B

Why did you only calculate for two terms? could be go on? can we do k = 3, 4, 5, 6... ?

I guess you mean the answer to part (b) is 14 because it's the smallest positive spectral component right?

& for (c) you mean f1 = 14 MHz = 70MHz - fs

For (a) you just need to know that sampling will create replicas of the original spectrum shifted by multiples of the sampling frequency. But don't forget the negative frequencies. Since $s(t)$ is real-valued, you have $|S(f)|=|S(-f)|$, and you also have to shift the contribution at negative frequencies (which is a mirror image of the part of the spectrum shown in the figure). Considering the frequency range [-70MHz,70MHz], replicas of the positive part of the spectrum (as shown in the figure) occur at $70MHz-k\cdot56MHz$, $k=1,2$:

$$f_1=14MHz,\quad f_2=-42MHz$$

Replicas of the negative part of the original spectrum in the range [-70MHz,70MHz] occur at $-70MHz+kf_s$, $k=1,2$:

$$f_3=-14MHz,\quad f_4=42MHz$$

So the answer to (b) is $f_1=14MHz=70Mhz-f_s$, with the last equality being the answer to (c).

Check this link for more information on bandpass sampling.


Steps for recognizing hand written signs What are the steps to built a system which recognizes hand written signs in an image and returns the most similar sign in Python? I know that depends on the image and many variables, but isn't there a sequence of steps, a main path to follow for recognizing images with at most 2 or 3 signs per image?

I thought the steps should be:

Take an image and convert to binary
Apply find_contour to binary image
Extract scale- and rotation-invariant properties
Create a feature vector with scale/rotation-invariant properties (if there is much more than 1 sign, like for example a star and a square not overlapped, is that a problem?)
Repeat steps 1 to 4 for every image
Put them into a k-NN or SVM

These next steps, I do not know how to implement:

I choose an image
System spits out the result (most similar image)

That's what you would do with images, not so much with handwriting, which needs more preprocessing and robust features because people don't write exactly the same way twice. See for example:

Improving Offline Handwritten Text Recognition with Hybrid HMM/ANN Models
Adaptive Membership Functions for Handwritten Character Recognition by Voronoi-Based Image Zoning
A visual approach to sketched symbol recognition
Feature extraction and classifier combination for image-based sketch recognition
Sketched symbol recognition with auto-completion
HBF49 feature set: A first unified baseline for online symbol recognition
Neural network-based symbol recognition using a few labeled samples


A Music Recommender System by Using Basis Functions and Inter Correlations So my university project is about music recommender system.
My teacher not saying too much. But he only said it will use basis functions and convolution technique.

I want some ideas about this project about how to create this.
I have programmed in MATLAB and know basics of DSP like DFT (FFT), Correlation, Convolution etc.

Any idea from experts?

I would like to inform that this music recommendation system will be content based and it will use correlation technique to find music similarity with basis signal.

Many Recommender Systems are based on Matrix Factorization.
There is a technique called Convolutional Matrix Factorization which is might what your teacher implied.

You may have a look here Tutorial - Large Scale Matrix Factorization (IEEE Big Data 2016 Tutorials) for Convolutional Matrix Factorization.


Implementation of Overlap Save Filtering with Decimation I am trying to implement Matlab code which performs FIR filtering with complex bandpass filter using Overlap Save method and perform decimation as well by using different length forward and inverse FFT as mentioned in this article:
https://www.iro.umontreal.ca/~mignotte/IFT3205/Documents/TipsAndTricks/MultibandFilterbank.pdf
Below is the code snip I have written:
clear all;close all;clc
fs=1e6;
sim_time=20e-3;
Nfft=256;

t=0:1/fs:sim_time-1/fs;
freq = fs/Nfft*(0:Nfft-1);

%create composite signal
fsig1= 50e3;
fsig2= 100e3;
fsig3= 300e3;
fsig4= 872e3;

sig1=2*cos(2*pi*(fsig1-5e3)*t)+4*cos(2*pi*fsig1*t)+8*cos(2*pi*(fsig1+5e3)*t) + j*(2*sin(2*pi*(fsig1-5e3)*t)+4*sin(2*pi*fsig1*t)+8*sin(2*pi*(fsig1+5e3)*t));
sig2=8*cos(2*pi*(fsig2-5e3)*t)+4*cos(2*pi*fsig2*t)+2*cos(2*pi*(fsig2+5e3)*t) + j*(8*sin(2*pi*(fsig2-5e3)*t)+4*sin(2*pi*fsig2*t)+2*sin(2*pi*(fsig2+5e3)*t));
sig3=2*cos(2*pi*(fsig3-5e3)*t)+8*cos(2*pi*fsig3*t)+4*cos(2*pi*(fsig3+5e3)*t) + j*(2*sin(2*pi*(fsig3-5e3)*t)+8*sin(2*pi*fsig3*t)+4*sin(2*pi*(fsig3+5e3)*t));
sig4=10*cos(2*pi*fsig4*t)+j*10*sin(2*pi*fsig4*t);

sig=sig1+sig2+sig3+sig4;

% Create Filter
bpFilt1 = designfilt('lowpassfir', 'FilterOrder', 40, ...
             'PassbandFrequency', 100e3, 'StopbandFrequency', 150e3,...
             'SampleRate', fs);

%filter impulse response
h1= bpFilt1.Coefficients;
%create complex band pass by multiplying with complex exponential
fc=800e3;
cmp_exp=exp(j*2*pi*fc*t);
h1_c=h1.*cmp_exp(1:length(h1));
%take FFT of filter impulse response
H1=fft(h1,Nfft)/Nfft;% low pass prototype
H1c=fft(h1_c,Nfft)/Nfft; %complex bandpass filter

%implement overlap save
P=41;%filter length

for i=1:10 % This code snippet compute the OSA filter with decimation by 4
    if(i==1)
        sig_win=sig(1:Nfft);
        index_fft_in=Nfft;
    else
        %sig_win=sig((i-1)*256-(P-1):(i-1)*256-(P-1)+Nfft-1); %overlapped input
        sig_win=sig(index_fft_in-(P-1):index_fft_in-(P-1)+Nfft-1); %overlapped input
        index_fft_in=index_fft_in-(P-1)+Nfft-1;
    end

    X_sig=fft(sig_win,Nfft)/Nfft;

    Y_sig=H1c.*X_sig;%multiply with complex bandpass filter
    y_sigd(i,:)=ifft(Y_sig((3*Nfft/4)+1:4*Nfft/4),Nfft/4)*(Nfft/2)^2;%extract complete for now to test filtering

    if(i==1)
        y_osa=y_sigd(i,((P-1)/4)+1:end);%discard the first (P-1) points
        index=length(y_osa);
    else
        y_osa(index:index+(Nfft/4)-(P-1)/4-1)=y_sigd(i,(P-1)/4+1:end);
        index=length(y_osa);
    end

end

% figure;plot(freq/1e3,abs(H1),'-*r')
% hold on;plot(freq/1e3,abs(H1c),'-ob')%plot filter responses

y1 = filter(h1_c,1,sig(1:10*Nfft)); % apply filter 1 %reference for comparison

%compare decimation in FD vs TD
y1_d = downsample(y1,4);

figure;
plot(real(y_osa),'-ob')
hold on;plot(real(y1_d),'-*g')

I am also doing the same filtering in time domain along-with downsampling to ensure that the time domain and frequency domain implementations are equivalent however there is some glitch in the output when compared as shown in image below:

Can anyone point out the mistake in the code or the method and suggest solution?
Edit: The first block of the result appears to match the reference time-domain filtered signal but some issue occurs in the subsequent blocks.
**Also, the result is fine when the same length of IFFT is performed (without any decimation) and the result matches exactly that of the time-domain reference signal. The issue is only when smaller IFFT size is used for required bins.

Attempt 2
I am not quite sure what is the interpretation of your filtering, in special, the IFFT part, you are taking a slice of the spectrum shifted by 3*Nfft/4, I think it should be (circularly) shifted by 7*Nfft/8.
  % your method
  xa = ifft(X((3*Nfft/4)+1:4*Nfft/4),Nfft/4) / 4;
  % what I would do
  xb = ifft(fftshift(fftshift(X)(3*Nfft/8+1:5*Nfft/8))) / 4;

Notice that you factor is (Nfft/2)^2 because you multiplied the spectra by Nfft.
To make sure I am not missing any constant I plotted the histogram of the correlations between the fft(trimmed(signal)) and decimated(ifft(signal)) for your method and the modified formula (correlation of 1 is exact up to an affine transformation).

The complete code I used to check is this
cora = []
corb = [];
for ii = 1:10000
  x = zeros(Nfft, 1);
  x(1:4:end) = randn(Nfft/4, 2)*[1;1j];
  x = cumsum(x);
  X = fft(x);
  % your method
  xa = ifft(X((3*Nfft/4)+1:4*Nfft/4),Nfft/4) / 4;
  % what I would do
  xb = ifft(fftshift(fftshift(X)(3*Nfft/8+1:5*Nfft/8))) / 4;

  % make sure this is not for a constant factor
  cora(ii) = abs(corr(xa, x(1:4:end))); % 0.66
  corb(ii) = abs(corr(xb, x(1:4:end))); % 0.98
end
figure
subplot(211),hist(cora, 25),xlim([0, 1]);
title('correlation histogram for your method'), xlabel('correlation')
subplot(212),hist(corb),xlim([0,1])
title('correlation histogram for the proposed method'), xlabel('correlation')

Working your example.
A few modifications that makes your frequency and time domain filtering consistent
clear all;close all;clc
fs=1e6;
sim_time=20e-3;
Nfft=256;

t=0:1/fs:sim_time-1/fs;
freq = fs/Nfft*(0:Nfft-1);

%create composite signal
fsig1= 50e3;
fsig2= 100e3;
fsig3= 300e3;
fsig4= 872e3;

sig1=2*cos(2*pi*(fsig1-5e3)*t)+4*cos(2*pi*fsig1*t)+8*cos(2*pi*(fsig1+5e3)*t) + j*(2*sin(2*pi*(fsig1-5e3)*t)+4*sin(2*pi*fsig1*t)+8*sin(2*pi*(fsig1+5e3)*t));
sig2=8*cos(2*pi*(fsig2-5e3)*t)+4*cos(2*pi*fsig2*t)+2*cos(2*pi*(fsig2+5e3)*t) + j*(8*sin(2*pi*(fsig2-5e3)*t)+4*sin(2*pi*fsig2*t)+2*sin(2*pi*(fsig2+5e3)*t));
sig3=2*cos(2*pi*(fsig3-5e3)*t)+8*cos(2*pi*fsig3*t)+4*cos(2*pi*(fsig3+5e3)*t) + j*(2*sin(2*pi*(fsig3-5e3)*t)+8*sin(2*pi*fsig3*t)+4*sin(2*pi*(fsig3+5e3)*t));
sig4=10*cos(2*pi*fsig4*t)+j*10*sin(2*pi*fsig4*t);

sig=sig1+sig2+sig3+sig4;
sig=randn(size(sig)); # white noise for testing
P=41;%filter length
h1 = exp(-((1:P)/10).^2); % (Octave) without designfilt

%create complex band pass by multiplying with complex exponential
fc=800e3;
cmp_exp=exp(j*2*pi*fc*t);
h1_c=h1.*cmp_exp(1:length(h1));
%take FFT of filter impulse response
H1=fft(h1,Nfft);% low pass prototype
H1c=fft(h1_c,Nfft); %complex bandpass filter
%implement overlap save
wLen = Nfft - (P - 1);
BLOCKS = 10;
y_osa = zeros(1, wLen * BLOCKS);
for i=1:BLOCKS % This code snippet compute the OSA filter with decimation by 4
    sig_win = sig(wLen*(i-1)+(1:Nfft));
    X_sig=fft(sig_win,Nfft);
    Y_sig=H1c.*X_sig;%multiply with complex bandpass filter
    y_sigd = ifft(Y_sig);
    y_osa(wLen*(i-1)+1:wLen*i)=y_sigd(P:Nfft);
end

y1 = conv(h1_c, sig(1:(BLOCKS+1)*wLen));
y1 = filter(h1_c, 1, sig);
y1 = y1(P:wLen*BLOCKS+P-1);

figure
hold on
plot(real(y1))
plot(real(y_osa))
```


Classifying and filtering out noise from zero crossings Note: I know absolutely nothing about zero crossings nor signal processing in general. However, I am doing programming exercises that require them. I've searched around but I cannot understand most of the information I'm getting. Hopefully I get clarification here. Feel free to close if it is off topic.

Here are two exercises from a computer science book (here https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-24.html#%_sec_3.5.3). I'll only include the relevant details.

  Exercise 3.74: Alyssa P. Hacker is designing a system to process
  signals coming from physical sensors. One important feature she wishes
  to produce is a signal that describes the zero crossings of the input
  signal. That is, the resulting signal should be 1 whenever the input
  signal changes from negative to positive, -1 whenever the input signal
  changes from positive to negative, and 0 otherwise. (Assume that the
  sign of a 0 input is positive.) For example, a typical input signal
  with its associated zero-crossing signal would be   … 1 2 1.5 1 0.5
  -0.1 -2 -3 -2 -0.5 0.2 3 4 … … 0 0  0  0  0   -1   0  0  0   0   1  0 0 …  ...

The associated second exercise is:

  Exercise 3.75.  Unfortunately, Alyssa's zero-crossing detector in
  exercise 3.74 proves to be insufficient, because the noisy signal from
  the sensor leads to spurious zero crossings. Lem E. Tweakit, a
  hardware specialist, suggests that Alyssa smooth the signal to filter
  out the noise before extracting the zero crossings. Alyssa takes his
  advice and decides to extract the zero crossings from the signal
  constructed by averaging each value of the sense data with the
  previous value. She explains the problem to her assistant, Louis
  Reasoner, who attempts to implement the idea, altering Alyssa's
  program as follows:

(define (make-zero-crossings input-stream last-value)
   (let ((avpt (/ (+ (stream-car input-stream) last-value) 2)))
     (cons-stream (sign-change-detector avpt last-value)
        (make-zero-crossings (stream-cdr input-stream)
                               avpt))))

  This does not correctly implement Alyssa's plan. Find the bug that
  Louis has installed and fix it without changing the structure of the
  program. (Hint: You will need to increase the number of arguments to
  make-zero-crossings.)

note: sense-data is the input signal

My questions are:

What is Exercise 3.75 referring to? Is it the large amounts of 0s?
Suppose we filter out the noise from the signals in "… 0 0  0  0  0   -1   0  0  0   0   1  0 0 …", should the resulting filtered signal be "… 0 -1   0  1  0 …" or "… -1   1  …"?
In this exercise, smoothing is done by averaging each value of the input signal with the previous value, for the life of me, I don't understand what that should do and why we do it?

Without knowing the context, here is what probably is meant:

3.75 refers to a more real-world sensor. The sensor described in 3.74 apparently is assumed to be perfect, ie whenever there is a negative/positive sensor data, the true signal does have the associated property. Ideal sensors, however, do not exist and 3.75 tells you, that your sensor does have imperfections that lead to spurious false sign alterations. The sensor shall be assumed to produce randomly additional zero-crossings that might look like that:

Instead of: 1 0  1 2 3  1 -5 -2 -2 -1 5 3
You get:    1 3 -2 3 2 -2 -4  3 -2 -1 2 3

Your example data suggested that a positive/negative periods are typically longer than just one measurement. You are supposed to "harden" your filter a bit to still provide the valid zero-crossings. This filter will never be perfect, but it should be better than the naive filter you probably implemented for 3.74.

In 3.74 you probably implemented a simple if-check: If the sign alternates between measurement $n-1$ and $n$, spit out the corresponding 1 or -1. A hardened filter could be implemented by looking back $n-k$ points with $k=1,2,3,...$ chosen appropriately and spit out a 1 or -1 if all of the $k$ points were of the same sign. This, however, would introduce a delay into your filter, and it would miss zero-crossings that are shorter than $k$ points. So, key would be to know the shortest duration between two zero crossings and choose $k$ smaller than that, and it would be important to know whether or not the delay is acceptable.

And no, noise does not refer to the zeros in your signal ;)


Detecting drum bpm in a noisy .wav file I am looking for algorithm(s) to solve the following problem:
Given a noisy .wav sound capture (some wind + friction noise on the microphone), how to detect the BPM of a soft drum beat?

I have attempted googling the subject, but the results are quite poor, due to high amount of mp3 related software for both analysis and fingerprint id generation. None of them supply information about how to actually do it.

I am aware of algorithms to remove the noise, but that still leaves me with the problem of detecting BPM. And depending on how the BPM problem is solved, it's possible that I don't even need to denoise (since drum tends to be in the lower frequencies and noise in higher, a simple low-pass might be sufficient pre-processing).

One method that works if there's a relatively strong drum beat is to take the magnitude of the STFT of the waveform, and then auto-correlate it in only the time dimension.  The peak of the auto-correlation function will be the beat, or a submultiple of it.

This is equivalent to breaking up the signal into a lot of different frequency bands, finding the amplitude envelope of each, autocorrelating each envelope, and then summing them.  The noise and other parts of music are averaged out by the cross-correlation operation.

This is because drum beats produce short-lived sound at many frequencies (vertical lines), while other parts of music are long-lived at only a few frequencies (horizontal lines), and noise is long-lived but random at all frequencies.  You can see the beat repetition if you look at an STFT:

I came up with this for a school project to find a single BPM value for entire music files, but it could be adapted to a stream of audio with changing BPM, too.  You'd need to process chunks that are at least twice as long as the period of the BPM you're looking for.


Conceptual Question from Signal Processing - Impulse Response and AR Coefficients In continuation to the previous question Conceptual questions from signal processing I have a doubt which is: Consider an Autoregressive model (AR(2)):
$$
y(t) = ay(t-1) + by(t-2)
$$

and a FIR (Moving Average, MA(2)) model

$$
x(t) = a\epsilon(t-1) + b\epsilon(t-2).
$$

According to the reply in the prev question, in time domain

$$
y[n] = h[n]\star x[n]
$$
$h$ is the impulse response.

Is there any relation between impulse response and the coefficients of AR and MA model?
What is the intuition of the coefficients and how do we get them?

Impulse Response is basically the FIR coefficients of the system.
Namely, a system $ H $ with an impulse response given by $ f [n] $ and a Filter $ F $ with an FIR representation of $ {f[0], f[1], \cdots, f[n]} $ are equivalent.

Now, systems with Feedback are equivalent of both FIR and IIR (AR) filters.
But given infinite length of FIR model any LTI system can be represented by FIR coefficients.

The relation could be easily displayed by the Laplace Transform of the system.


Determine whether a signal is increasing or decreasing over time I have a signal (x-axis: Time, y-axis: Data) which has more than 10000 samples per second.
The signal itself represents some physical data (e.g. speed, acceleration, etc.) and has some noise on it.

Details about the considered signal:

The signal increases/decreases over the considered capture interval with a small gradient
Let's consider the case where the signal decreases:

Due to the arbitrary noise (which has a very small amplitude) it regularly happens that one (or more consecutive) samples are not smaller than their predecessor.
However due to the high sampling rate those deviations are just an issue for the evaluation routine (Python) trying to decide whether the signals is continously decreasing in the considered time interval.
For a human being, with optical inspection of the plot it is very clear that the signal is decreasing

What is the most robust way to determine, whether the signal is increasing or decreasing over time?

I thought about using a moving average filter but I am curious to know whether there is a more proper way to reach my goal.

Any answer is highly appreciated.

You may use a bandlimited differentiator (i.e. the cascading of a low-pass filter to remove noise, and a differentiator to get the value of the slope).


Looking for pratical quantitative comparison metrics for scaled, delayed and warped Signals [Nota: I have validated the only answer so far, yet more explicit versions are still welcome]
The following question is detailed in 1D, with time as the ordinal variable. Similar questions could apply in other dimensions.
In several signal processing techniques, such as blind source separation (BSS), filter banks, or deconvolution, one may wish to estimate a signal $x(t)$ and only recovers $s.x(t+d)$, a scaled and delayed estimate. Rotations and shears can be added in higher dimensions, and many other more.  $s$ is a scale factor, $d$ a delay. One might even stumble with warped data ($x_{s,d,w} = s.x(t/w+d)$), as  in super-resolution for instance.
In theory, one can estimate continuously $s$ and $d$ with local correlation or Fourier transforms (How to match 2 signals that have same information, though shifted and scaled). The warping $w$ might be estimated with the scale-transform or wavelet representations. I have read several BSS papers and books, asked people, been in conferences, and was not able to find a standard, or at least a usable metric.

In image (it works on signals as well), the Structural Similarity index somehow compensates offset and variance.

Are there practical error metrics to compare the original $x(t)$ with the transformed $x_{s,d,w}(t)$ in the context of sampled signals and noise conditions? Indeed, the discretization induced by sampling complicates the comparison task (imagine for instance a $1$-sample spike on the sampling grid that would be delayed by a non-integer time), as well as the noise.
Should one resort to asymmetric quantities such as divergences?
Can other signal properties help (bandpass, sparse, positive, etc.)?

Forgetting about the warping, I have tried to minimize a standard $\ell_p$ norm, with $s$, $d$, and $w$ as parameters, and to smooth both signals. I am not satisfied with the complexity and the outcomes, and this is a bit tedious.

I'm answering the question the way I understood it - How can one find a similarity measure which isn't sensitive to scaling and shifting.

An approach could be borrowed from the Computer Vision world by comparing Shift and Scale Invariant features between the two signals.

I'm not sure it will work for measuring the quality of recovering signals but it certainly will tell two signals are similar given the hierarchy between their features is similar and the features themselves are similar.


Scipy Butterworth highpass filter results in distortion and bloated output Trying to use Scipy Butterworth filter to filter frequencies below 54 Hz from audio files. Audio files are stored as WAVs. Sample rate is 44100.
When I try using this filter the result is distorted (very). Also, I assume with fixed time and sample rate that input and output files should be exactly the same size. Input file is 441 Kib, but output is 1.8 Mib. But the duration is unchanged.
Problem is with the filter or with application of the filter. I'm opening and saving WAV files with scipy.io.wavfile and have tested without applying the filter -- what I get out is what I put in, as expected.
Here's my code:
import scipy.io.wavfile
from scipy import signal

# Scipy reads audio data from WAV files
rate, data = scipy.io.wavfile.read('./test.wav')
# rate is the sample rate, data is the data
# NOTE: Sample rate of my input is 44100
assert rate == 44100

cutoff = 54  # want this to be 54 Hz
nyquist = 0.5 * rate
normal_cutoff = cutoff / nyquist
order = 5

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html
# Create a filter and return coefficients
b, a = signal.butter(order, normal_cutoff, btype='highpass', analog=True)

# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html
# Apply filter to audio signal
filtered = signal.filtfilt(b, a, data)
scipy.io.wavfile.write('filtered.wav', rate, filtered)

Have also tried filtered = signal.lfilter(b, a, data) for application of filter, and have tried analog=False. Behavior is similar in both cases.
What am I doing wrong? Does this have to do with the sample rate? Bit depth? I'm stuck.

For digital filters, Wn are in the same units as fs. By default, fs is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (Wn is thus in half-cycles / sample.)
For analog filters, Wn is an angular frequency (e.g. rad/s).
As yours is an analog filter, your Wn parameter is 2π·54Hz
from scipy import constants
...
Wn = 2*constants.pi*54
b, a = signal.butter(order, Wn, btype='highpass', analog=True)


DCT step in MFCC process During MFCC process FFT is done to get the power spectrum but last step is DCT in order to get back to time domain. Cannot DCT be used in the first case to convert the signal to frequency domain?

DCT is not used only for going back to time domain, it is also used due to it's "de-correlating properties". By saying that I mean compression of most of the information in smaller number of coefficients. Although in theory we could use the DCT to calculate the spectrum, but I never saw anyone doing that. Usually thing is simple:

want to analyse the spectrum? Use DFT
want to compress and de-correlate data? Use some flavour of DCT

Also one very important thing. In case of DCT, the equivalence property of convolution and multiplication in frequency domain doesn't hold:
$$x[n]\star y[n] \nLeftrightarrow X_{DCT}[k]\cdot Y_{DCT}[k] $$
which is the fundamental property (homomorphic) used when calculating the MFCC's. It allows to perform the Cepstral Mean Normalization and remove the "convolutional distortions", such as vocal tract. For more info please refer to my previous answer.


Is resizing of frequency spectrum a valid method of resampling? One method of resampling is to perform a Fourier Transform on a signal, resize the resulting frequency spectrum, and then return to the resampled version of the signal with an inverse fourier transform.
Is this a method used in practice?
Are there methods that perform better down/upsampling (it terms of interpolation error)?
Are there faster methods?

Interpolation by zero padding and computing the IDFT is used in practice but is not optimum and will have greater error over other approaches with the same or less resources, for the same reason that filtering by zeroing frequency bins in a DFT is not recommended. This point is detailed at this post Why is it a bad idea to filter by zeroing out FFT bins?.
Methods that perform better in terms of efficiency and performance include polyphase resampling structures using filter coefficients based on least squares algorithm, and specifically for multiband fiters that concentrate rejection at the image locations associated with the spectral distortion.
For further details on such approaches see:
How to implement Polyphase filter?
and
Computational Complexity of Polyphase Resampling
When a Sinx/x droop can be tolerated (or easily compensated for, see how to make CIC compensation filter) cascade-integrator-comb (CIC) interpolators are highly attractive for their simplicity and efficiency (there are many posts here that detail the beautiful CIC (otherwise called the Hogenhauer filter after the original author) further, see:
CIC Cascaded Integrator-Comb spectrum
Puzzled by CIC filters
Interpolating/Decimating CIC Filter Group Delay


Clarification on eigenfunction property from Oppenheim's Discrete-time signal processing, 3rd ed Could anybody explain what is meant  by the following statement on page 50 of Oppenheim and Schafer's book Discrete-Time Signal Processing (third edition)?

The eigenfunction property of complex exponentials depends on stability
of the system, since at finite $n$, the transient response must have become zero, so that we only see the steady-state response $H(e^{j\omega})e^{j\omega n}$ for all finite $n$.

This statement can only be understood given its context. Right before the sentence you quoted, we have this statement:

It is no surprise that the condition for existence of the frequency
response is the same as the condition for dominance of the steady-state solution.
Indeed, a complex exponential that exists for all $n$ can be thought of as one that is applied at $n = -\infty$.

Consequently, starting out at $n=-\infty$, if the system is stable, all transients must have died out at "finite $n$" because at finite $n$ the moment that the exponential was applied is infinitely long ago.
We also depend on stability because otherwise the eigenvalue $H(e^{j\omega})$ wouldn't exist, at least not in the classical sense.


How to find SNR of received DSB:AM signal with white noise after demodulation & LPF'ed? This is a homework problem so please do not give a full solution. I intend to keep this at discussion-level.

Problem Statement

I am given the Power Spectral Density (PSD) of a DSB-Amplitude modulated signal and am told that it is corrupted with additive noise with PSD = N/2 "within the passband of the signal". That signal (with noise) is demodulated and passed through a Low-Pass Filter (LPF) whose frequency spectrum is unity between -W and W. I am tasked to find the SNR at the output of the LPF.

Attempt

SNR is defined to be the ratio of (1) the power of the message component of the demodulated received signal to (2) the power of the noise component of the same demodulated received signal.

We just need to take the ratio: (1)/(2).

Assuming the answer to question 2 (below) is yes. The power of the white noise component between -W and W is $2W*N/2 = WN$.

Questions

I'm given the (plot of the) PSD of the modulated signal. Is this (1)?
How does LPF affect the PSD of a signal? (Suppose we're given the plot of both the LPF (unity gain) passband and the PSD of the signal. Do we just point-wise multiply the PSD by the LPF function?).

Since you don't want a full answer, I'll just provide a couple of hints. Your calculation for the noise power is correct. Regarding your questions:

The power of a signal is the integral of its PSD.
You're looking for the Wiener-Khinchin theorem, and especially its application to LTI systems.


Question about convolution in time and frequency domain My question is about the number of data points when doing a convolution (or correlation) in the the time or frequency domain. let's say we have two signals, $a(t)$ and $b(t)$ each of length $16$. When we do the convolution as follows:

$$c(t) =a(t)\star b(t)$$

$c(t)$ consists of $31$ points, right? Now if we transform these into $A(f) = \mathrm{FFT}(a(t))$ and $B(f) = \mathrm{FFT}(b(t))$, $A(f)$ and $B(f)$ each consist of $16$ points still, were the first half is a mirror of the second half. The convolution is just a multiplication in the frequency domain:

$$C(f) = A(f)\cdot B(f)$$

I believe $C(f)$ should be $16$ points. I'm fine with that much, but if we transform $C(f)$ back to the time domain, doesn't  $c(t)$  end up having $16$ points?

I'm sure I'm just not thinking correctly about this. Thanks for any clarification!

Normally, the variables $t$ and $f$ are used for continuous time and frequency. But from your question I understand that you're talking about discrete time and frequency, and their relation via the Discrete Fourier Transform (DFT).

If you have discrete-time sequences $a[n]$ and $b[n]$ and their DFTs $A[k]$ and $B[k]$, then the DFT of the linear convolution $c[n]=a[n]\star b[n]$ is NOT equal to $A[k]\cdot B[k]$. The product of the DFTs corresponds to circular (or cyclic) convolution in the time domain.

The linear convolution of two $16$-point sequences has indeed $2\cdot 16-1=31$ points, whereas the circular convolution of two $16$-point sequences also has $16$ points.