Detection of low frequency signal I have a strong primary signal with frequency range between 100–200 kHz. This signal has very low frequency and lower SNR secondary signal phase modulated to it.
The modulated signal is around 10-30 Hz.

The fundamental signal frequency changes over time, fast enough that I cannot take enough samples to observe the secondary frequency. For example signal goes from 200 to 170 kHz in short enough time that if I take a simple FFT, the peak will be spread from 200 to 170 kHz and secondary signal will be buried inside.

The good news is that the second signal frequency doesn't change with fundamental frequency change.

I like to detect this secondary, modulated lower SNR signal with 1 Hz accuracy.

What should be my sampling frequency for a given maximum 200 kHz primary signal?
What should be the DSP approach I should take to determine the exact frequency of the modulated signal?

The hard one first

  What should be my sampling frequency for a given maximum 200 kHz primary signal?

Assuming your sampling is real-valued, $f_s > 2(200 + 30)\text{ kHz}$, as per the Shannon-Nyquist theorem, would sound sufficient – but the problem is that it isn't, for reasons of ambiguity, see below.

  What should be the DSP approach I should take?

First, model your resulting signal; this isn't really hard to do here. People mean different things when they say "modulation", but I'm going to go ahead and assume what seems most likely: primary signal being modulated by the secondary means that your resulting signal contains the sum frequency.

Hence, spectral estimation is really what you're after.

Furthermore, you need an estimator that is more efficient than the DFT, but doesn't have to be able to estimate a lot of spectral components – in fact, you're only interested in a single one, namely the sum of secondary and primary frequency.

I'd recommend using ESPRIT, parameterized to a single tone to detect.

Whatever spectral estimator you use, they will expose a specific variance that will typically decrease with the number of samples you take into consideration – relate your desired spectral resolution (1 Hz) to this variance, and you'll get a number of samples you'll need to get a sufficiently "exact" estimate. That number of samples will imply a sampling rate, which probably will be much higher than the aforementioned Nyquist bound.

Now, if your signal only contained the perfect primary modulated with the clean secondary signal, you'd be lost here – because if you found out your resulting signal was at let's say 150.3 kHz, you couldn't tell which part of that would be primary and which part secondary. It's mathematically impossible to tell a variation of the primary frequency from a variation of the secondary, lest your modulation doesn't just add up the frequencies. Hence, just fulfilling Nyquist for your "clean" signals won't help you – there's simply no information in these that allow you to discriminate between primary and secondary signal.

Luckily, you have noise. You don't tell us anything about that, but I'll assume the following:

It's additive to your secondary signal,
it has zero mean (which would be good from a measurement quality point of view; otherwise, we can model it as DC + zero-mean process),
it's observation's representation is symmetric in spectrum,
it's uncorrelated to both the primary and the secondary signal.

Being additive to your secondary signal implies that you'll get primary signal multiplied with noise. That means we get another, albeit "unclean" "tone"!

So we get two spectral components in our resulting signal:

A tone with frequency $f_\text{primary}+f_\text{secondary}$,
an area around the primary frequency, spectrally shaped by convolution of the spectral shape of the noise and the dirac of the primary tone

Now, 2. is a thing hard to estimate, right, because it's all "blurred out" in spectrum?

The good thing is that with sufficient oversampling and by subtracting the result of subtracting the tone estimated by ESPRIT we can find the "center" of the second component in spectrum – and that's where our "clean" primary signal was! Now, with that frequency, and our ESPRIT-estimated sum frequency, we can find the difference frequency – which is the frequency of the secondary signal.

However, your 1 Hz resolution is very unlikely to happen, for various reasons. In fact, I'm a bit too lazy to do the math now, but I don't think the time$\times$bandwidth product of your observation would allow for that. Also, the above considerations neglect any aspect of how much noise power is actually in your secondary signal, and how much additive noise you'll also see on your overall signal. My guess is that additive noise in your observation will be in the same order of magnitude as the multiplicative noise, making the estimation of a spectral center much less clear.

Also, uncorrelated multiplicative noise is usually a gross simplification – in fact, it's not how active modulators work at all.

So, really, before you try to build a system that does what you want, try to disprove the claim that you mathematically can do it at all, based on estimates of noise and signal powers, and noise bandwidths. If disproving the possibility is hard, or you get a proof that it should be possible, then go for an implementation.
If you find a proof that 1 Hz resolution is impossible, don't waste your time implementing such a system. It's natural that some phenomena aren't possible to observe with a given certainty, and that's one of the core realizations that digital communication engineers have to suffer through.


Side Lobes in the Magnitude Response of a Low Pass FIR Filter I have the impulse response (from the filter coefficients) of an FIR filter obtained from MATLAB's "interp" function using the command:
[y, b] = interp(x, 5); % where b contains the interpolation filter coefficients.
When I plot the magnitude response of this filter using "freqz" in MATLAB, I get the following

What do those side lobes between 0.5 and 0.7 (x pi rad/sample) and 0.9 and 1 (x pi rad/sample) mean for my FIR filter? Why are they there?

The filter is optimized to suppress signal components in the frequency ranges where the side lobes are low. The high side lobes in are just a consequence of optimizing the stopband suppression in those other frequency bands, where aliasing components are expected. The regions with high side lobes are basically "don't care" regions. If you wanted to increase stopband attenuation there, the stopband attenuation in the other - more important - regions would be decreased.


How to calcuate the beta metrics? I am working on Speckle reduction in ultrasound images. I read a paper (snapshot of that paper is included below) where they have used a metric called beta metric which indicates the edge preserving capability of the filters. I tried to implement it using the equation (7) provided in the below snapshot.

This is the code.

 %---Beta Calculation
    h = fspecial('laplacian');
    I1 = imfilter(Orig_Image,h);
    I2 = imfilter(Esti_Image,h);
    I_1 = mean2(I1);
    I_2 = mean2(I2);
    metrics.Beta = sum(sum((I1 - I_1).*(I2 - I_2)))./(sqrt(sum(((I1 - I_1).^2).*((I2 - I_2).^2))));

As you can see that the beta value is a single value [1x1] vector in Table 2 inside the snapshot. But when I run my code, I get beta value of size 1x256 for a 256x256 image.
Is my code wrong? How to implement the given equation (7) correctly in MATLAB? Please help.

The correct code is

metrics.Beta = sum(sum((I1 - I_1).*(I2 - I_2)))./(sqrt(sum(sum((I1 - I_1).^2))).*(sum(sum((I2 - I_2).^2))));

The reference you have shown is incorrect, reffr paper:

  2013 August Gregorio - ‘A suitable threshold for speckle reduction in
  ultrasound images’,  IEEE transactions on instrumentations and
  measurement


Huffman coding for $8\times 8$ blocks of an image I have a question regarding the Huffman coding in Image processing. In order to compress the signal, we usually operate on $8\times 8$ pixels (blocks) on which we perform DCT, zig-Zag scan, quantization then construct the alphabets using run-length encoding (RLE). After all of that, we do Huffman coding.

My question is: Do we have to compute and store the Huffman tree/table for each block separately and augment it to the bit stream of the block or we compute the frequency of the alphabets in all the blocks jointly and do Huffman coding only once at the end (this saves the memory needed to store different table/tree for each block).

In other words, isn't it almost useless to use Huffman coding on only $64$ coefficients ($1$ block) and save the table with data stream (which consumes memory)?

Huffman coding, as an example of statistical entropy coding, is used as the last step of still image compression to further reduce any statistical redundancy (possibly) present at the output representation of quantized DCT coeffcients of a given block of N by N pixels (usually N = 8)

When the statistical behaviour of those DCT coefficients are analysed throughout the whole image, it will be seen (for naturally generated images) that certain symbols occur more frequently to the extend that using a variable length code (VLC) will become benefical instead of a fixed length one.

Thus to utilize this further possibility of statistical compression, VLC is employed via Huffman code tables for each of the AC, DC, Gray and Color channel coefficients separately according to their respective statistical symbol frequency distributions.

Now, the proper way of generating a statistically optimized Huffman table per block, actually necessitates a per-block optimized huffman table as you outlined in your question. However due to several reasons, one of which is again mentioned as the complaint in your question, such a method would oppose the initial purpose of image compression; to reduce its size.

Mainly:
1- inlcusion of per block optimized huffman tables into each block will reduce benefits of compression (it can even increase the block size for high quality compression where there will be so many symbolls and the table will include as many data as the block itself)

2- The loss of the benefit from discarding a per-block optimized Huffman table utilization will be mostly insignificant (i.e. A single whole image optimized Huffman table will perform almost as good as a per-block optimized one - even when excluding the effect of the table storage per block)

3- It will be more robust to create the table from an ensamble of blocks rather than a single one. (However note that due to statitical nonstationarity of typical images, the statistical behaviour from one block to an adjacent one other will change and hence the ensamble is not always a qualified source of data for table optimization)

4- Nevertheless, the approach taken is to use a single Huffman table for the whole image and store this inside the image file headers.

5- Now this single Huffman table can be defined in one of two ways: as either a pre-defined standard table or as an image specific optimized table. Most often the former approach is taken as the benefit from using specific image optimized Huffman table is attenuated by the table processing time and the complexity of programming required for per-frame Huffman table optimization.


Losslessness of Laplacian Pyramid I was just reading up on the Gaussian pyramid and the Laplacian pyramid, used in compression applications of image. The source is here - Carnegie Mellon 16-385 Computer Vision - Spring 2019, Lecture 3 - Image Pyramids and Frequency Domain.

The slide claims that during Laplacian Pyramid based compression, the original image can be reconstructed back perfectly, because of saving the residuals images.

My understanding of the process is as follows

Convolve original scale image with a Gaussian
Compute the residual image, and save it
Downsample the convolved image,
dropping alternate row and column to get to the next level (lower
size image).
Repeat process at this scale, till you reach lowest resolution desired.

Now, it says that given the lowest size image, and the residuals, we can reconstruct back the highest size image.

However, since we are dropping alternate  rows and columns while downsampling, isn't this information lost? How can the reconstruction be perfect?

Actually the down sampling has no role here.
It is all based on a real simple equation:

$$ I = A + B $$

It is always enough to keep 2 terms of the 3 to restore completely and perfectly the information.

So let's look on this:

$$ {I}_{0} = \left( \left( {I}_{0} \downarrow \right) \uparrow \right) + {R}_{0} $$

So if we keep $ {R}_{0} $ and we have $ \left( \left( {I}_{0} \downarrow \right) \uparrow \right) $ we can recover $ {I}_{0} $ right?
But we can fully recover $ \left( \left( {I}_{0} \downarrow \right) \uparrow \right) $ form $ \left( {I}_{0} \downarrow \right) $ as it is deterministic known operation (Up sampling) right?

So instead of $ {I}_{0} $ we will keep $ {I}_{1} = \left( {I}_{0} \downarrow \right) $ and $ {R}_{0} $.

We now can repeat the same trick for $ {I}_{1} $:

$$ {I}_{1} = \left( \left( {I}_{1} \downarrow \right) \uparrow \right) + {R}_{1} $$

So instead of keeping $ {I}_{1} $ we can again keep $ {I}_{2} = \left( {I}_{1} \downarrow \right) $ and $ {R}_{1} $.
So in order to resotre $ {I}_{1} $ perfectly all we need is $ {I}_{2} $, $ {R}_{1} $ and $ {R}_{0} $.

Keep doing this and you will get the Laplacian Pyramid which says keep the residual of each level and the Low Pass image of the last image.

The trick here isn't the accurate definition of $ \downarrow $ or $ \uparrow $ but being able to get the same result every time you use them. Namley:

$$ \left( \left( \left( \left( I \downarrow \right) \uparrow \right) \downarrow \right) \uparrow \right) = \left( \left( I \downarrow \right) \uparrow \right) $$


In the context of transfer functions, what is the relationship between the terms "proper", "causal", and "realizable"? I am thinking about these terms in the context of linear control.

A transfer function is proper if the degree of the numerator is not greater than the degree of the denominator. I've read often that improper transfer functions are "not causal". I also often see the word "unrealizable" used often in this context.

If a control transfer function I've designed is improper, does that mean it is "causal" and/or "unrealizable"? What is the difference between these terms? What do they mean in practice?

Causality is a necessary condition for realizability. Stability (or, at least, marginal stability) is also important for a system to be useful in practice.

For linear time-invariant (LTI) systems, which are fully characterized by their transfer function, we get realizability constraints on the transfer function. For continuous-time LTI systems, if we work at frequencies for which the lumped element model is valid, we require the system's transfer function to be rational for the system to be realizable. Also for discrete-time LTI systems we require rationality of the transfer function, which implies that the system can be realized by adders, multipliers, and delay elements.

For an LTI system to be causal and stable, its poles must lie in the left half-plane (for continuous-time systems), or inside the unit circle (discrete-time systems). From this it follows that the rational transfer function of an LTI system must be proper, otherwise you would get one or more poles at infinity, causing the system to be unstable (or non-causal).


Fast Fourier Transform for varying period (Order analysis) I am trying to do a fast Fourier transformation on accelerometer data from a shaft rotating at varying speed in Python.

What I have done so far:

1: The original plot was in the time domain, and I therefore did a order analysis (resampled), and got the following plot:

This plot shows the angular rotation plotted against amplitude.

2: Now, an FFT was done with this code:

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

class FastFourierTransform:
    # Amplitudes is a row vector
    def __init__(self, amplitudes, t):
        self.s = amplitudes
        self.t = t

    # Plotting in the input domain before fft
    def plot_input(self):
        plt.ylabel("Amplitude")
        plt.xlabel("Shaft angle [Radians]")
        plt.plot(self.t, self.s)
        plt.margins(0)
        plt.show()

        '''
        The second half of this array of fft sequence have similar frequencies
        since the frequency is the absolute value of this value.
        '''
    def fft_transform(self):
        mean_amplitude = np.mean(self.s)
        self.s = self.s - mean_amplitude # Centering around 0
        fft = np.fft.fft(self.s)

        # We now have the fft for every timestep in out plot.

        # T is the sample frequency in the data set
        T = self.t[1] - self.t[0] # This is true when the period between each sample in the time waveform is equal
        N = self.s.size  # size of the amplitude vector
        f = np.linspace(0, 1 / T, N, )  # start, stop, number of. 1 / T = frequency is the bigges freq
        plt.ylabel("Amplitude")
        plt.xlabel("Frequency [Hz]")
        y = np.abs(fft)[:N // 2] * 1 /N

        # Cutting away half of the fft frequencies.

        sns.lineplot(f[:N // 2], y) # N // 2 is normalizing it
        plt.margins(0)
        plt.show()
        time = f[:N // 2]
        return fft, time

3. Result, with the normalized amplitudes plotted:

Questions:

Does this thought process look correct?
Is it correct to say that the final fft plot is in the frequency domain? From this link, http://zone.ni.com/reference/en-XX/help/372416L-01/svtconcepts/svcompfftorder/, it looks like the final plot domain should be in the order domain, but I am not sure since the fft was
done from the radian domain.

Vibration- Based Condition Monitoring of Wind Turbines by Tomasz Barszcz has this image

Thanks in advance.

There are a few things going on here:

  Does this thought process look correct?

Strictly, No.

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

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

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

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

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

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

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

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

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

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

So what?

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

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

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

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

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

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

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

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

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

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

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

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

And here is its spectrogram:

(Both images are from the Wikipedia article).

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

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

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

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

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

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

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

For more information please see here and here

Hope this helps.


Making the units of a analytically calculated PSD consistent with the units of an FFT One thing that has always puzzled me is making the units of an analytically derived Power Spectrum Density (PSD) consistent with the units of an FFT.
Let's say we record the output of a frequency generator, which we just set to output a signal of the form
$$v(t) = V_{\rm{pk}} \cos (2 \pi f_0 t)$$
for some finite amount of time, $T$, and take the FFT of this  acquired signal. We might then expect to see something like this:

Now as I understand it the FFT is a Linear Amplitude Spectrum (LAS) which will have units of Volts, or, a Power Spectrum (PS) which will have units of Volts-squared, depending on what we choose do with the resultant FFT. But it is clear the FFT can only have units of volts because that is what the units of the time-transient signal are.
Now lets say I want to fit this FFT with some realistic function. This is straight forward as I have the transient signal and I can calculate, analytically a PSD from a time truncated Fourier transform:
$$\frac{1}{T}\left|\int_{0}^{T} v(t) e^{-2 \pi j f t} \ {\rm{d}}t \right|^{2}$$
which will be in units of Volts-squared per Hertz ($\rm{V^{2}/Hz}$).
So its clear that the units of the FFT and the units of my function I want to fit to it don't match up! Multiplying by some quantity in units of $\rm{Hz}$ fixes the problem for example either the bandwidth of the FFT, $\Delta f = 1/ T$ or the Effective Noise Bandwidth (ENBW).
How can i reconcile the units of my analytic expression and the units of my FFT?

The objective is to simply fit my FFT data, so I given that all this should only really affect the amplitude, I could of course allow all this to be absorbed by an arbitrary constant -- but it would be nicer to have a consistent result.

You do this by doing two things:

First establish what the relationship of scale exists between your A/D converter and the numerical values that go into your DFT.  Find out how many volts corresponds to a value of $1.0$ in the samples that go into the DFT.

The next thing you gotta do is express your integral above without that $\frac{1}{T}$ factor (which is an error) as a Riemann summation having equal-width little rectangles of width equal to the sampling period.  That will become identical to the summation of the DFT.

That should spell it out.


Eliminate Signal A from Signal B i have 2 audio signals, both available in time and frequency domain:
A: Noise
B: Noise(A with minor variations) + Additive sounds (C)

i want to eliminate Signal A from B to leave only C.
C is for testing almost zero, so after removing A from B silence would be the best.

My signals to visualize the problem:
Signal A:

Signal B:

Now my questions:
What would be the best approach to eliminate A from B?

I have tried to implement power spectral subtraction with no results as the signals seem to bee not similar enough to just subtract them. (or i did it wrong, i this case just tell me if this is the best way to do it)

As a spectrum of C would be sufficient a more graphical approach may be better? To smooth my fft data, compute power spectra (basically what audacity does in the pictures) and subtract them to see whats left? (should be C)?

Additional question:
how can i change the overall amplitude of a signal in frequency domain without "colouring" the sound? i tried multiplying all real and imaginary parts by x but this changes the sound.

any help appreciated!
regards

The fact that you are plotting the FFT of the whole audio clip makes me think that you are looking at the wrong class of solutions. From the waveforms of signal A and B, it is clear that your noise and signal are not stationary. If you want to use a frequency-domain noise reduction technique, this should be done, in your case, on shorter overlapping windows across which you can assume that the properties of your signal and noise are stationary. You cannot avoid using the Short-term Fourier Transform.

In the rest of this post, I will switch notations to:

$X$ original (not noisy) signal
$N$ noise signal
$Y = X + N$ your observation, that is to say $B$. $X$ is what you want to recover, $N$ in your case is known (signal $A$)

Most denoising methods share the same structure: get a STFT representation of $Y$ and $N$ (matrices $Y(p, k)$ and $N(p, k)$, where $p$ is the frame index and $k$ the frequency bin index), process pairs of FFT frames from rows of $Y$ and $N$ to get an estimate $\hat{X}$ of a slice of the denoised signal spectrum, and use it to build a filter/mask applied in the frequency domain to $Y$, the conversion back to the time-domain being performed through overlap-add.

I suggest you to try to implement the Ephraim-Malah algorithm, which is a classic and often cited "baseline" denoising method. It is based on two intuitive tricks:

Estimate the SNR and attenuate the strength of the denoising filter in segments with low SNR (this prevents many artifacts in areas where the signal is close to the noise level).
Use temporal smoothing (if we have a high SNR in a frame, make the filtering/subtraction stronger in the next frame since it is also likely to contain signal...).

I have found a relatively self-contained set of slides on the topic - the author introduces increasingly sophisticated denoising rules, from spectral subtraction to full-blown Ephraim-Malah. The only thing he doesn't tell you to get to a practical implementation is how to get the estimate $\hat{S}_x(p,k)$ - since you don't observe $X(p,k)$, just $Y(p,k)$ and $N(p,k)$. The simplest approach is to subtract the power spectrum of the noise from the power spectrum of the signal, and zero the STFT cells in which the result is negative:

$\hat{S}_x(p,k) = \begin{cases} |Y(p,k)|^2 - |N(p,k)|^2 & when |Y(p,k)|^2 - |N(p,k)|^2 > 0 \\ 0 & otherwise \end{cases}$

You can then try the filtering rules shown on pages 5 and 10 of the slides.


Directly compare subpixel shifts between two spectra — and get believable errors I have two spectra of the same astronomical object. The essential question is this: How can I calculate the relative shift between these spectra and get an accurate error on that shift?

Some more details if you are still with me. Each spectrum will be an array with an x value (wavelength), y value (flux), and error. The wavelength shift is going to be sub-pixel. Assume that the pixels are regularly spaced and that there is only going to be a single wavelength shift applied to the entire spectrum. So the end answer will be something like: 0.35 +/- 0.25 pixels.

The two spectra are going to be a lot of featureless continuum punctuated by some rather complicated absorption features (dips) that do not model easily (and are not periodic). I'd like to find a method that directly compares the two spectra.

Everyone's first instinct is to do a cross-correlation, but with subpixel shifts, you're going to have to interpolate between the spectra (by smoothing first?) -- also, errors seem nasty to get right.

My current approach is to smooth the data by convolving with a gaussian kernel, then to spline the smoothed result, and compare the two splined spectra -- but I don't trust it (especially the errors).

Does anyone know of a way to do this properly?

Here is a short python program that will produce two toy spectra that are shifted by 0.4 pixels (written out in toy1.ascii and toy2.ascii) that you can play with. Even though this toy model uses a simple gaussian feature, assume that the actual data cannot be fit with a simple model.

import numpy as np
import random as ra
import scipy.signal as ss
arraysize = 1000
fluxlevel = 100.0
noise = 2.0
signal_std = 15.0
signal_depth = 40.0
gaussian = lambda x: np.exp(-(mu-x)**2/ (2 * signal_std))
mu = 500.1
np.savetxt('toy1.ascii', zip(np.arange(arraysize), np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth), np.ones(arraysize) * noise))
mu = 500.5
np.savetxt('toy2.ascii', zip(np.arange(arraysize), np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth), np.ones(arraysize) * noise))

I think using cross-correlation and interpolating the peak would work fine.  As described in Is up-sampling prior to cross-correlation useless?, interpolating or upsampling before the cross-correlation doesn't actually get you any more information.  The information about the sub-sample peak is contained in the samples around it. You just need to extract it with minimal error.  I gathered some notes here.

The simplest method is quadratic/parabolic interpolation, which I have a Python example of here.  It's supposedly exact if your spectrum is based on a Gaussian window, or if the peak happens to fall exactly on the midpoint between samples, but otherwise has some error.  So in your case you probably want to use something better.

Here is a list of more complicated, but more accurate estimators.  "Of the above methods, Quinn's second estimator has the least RMS error."

I don't know the math, but this paper says that their parabolic interpolation has theoretical accuracy of 5% of the width of an FFT bin.

Using FFT interpolation on the cross-correlation output does not have any bias error, so that's the best if you want really good accuracy.  If you need to balance accuracy and computation speed, it's recommended to do some FFT interpolation and then follow it with one of the other estimators to get a "good enough" result.

This just uses the parabolic fit, but it outputs the right value for the offset if the noise is low:

def parabolic_polyfit(f, x, n):
    a, b, c = polyfit(arange(x-n//2, x+n//2+1), f[x-n//2:x+n//2+1], 2)
    xv = -0.5 * b/a
    yv = a * xv**2 + b * xv + c

    return (xv, yv)

arraysize = 1001
fluxlevel = 100.0
noise = 0.3 # 2.0 is too noisy for sub-sample accuracy
signal_std = 15.0
signal_depth = 40.0
gaussian = lambda x: np.exp(-(mu-x)**2/ (2 * signal_std))
mu = 500.1
a_flux = np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth)
mu = 500.8
b_flux = np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth)

a_flux -= np.mean(a_flux)
b_flux -= np.mean(b_flux)

corr = ss.fftconvolve(b_flux, a_flux[::-1])

peak = np.argmax(corr)
px, py = parabolic_polyfit(corr, peak, 13)

px = px - (len(a_flux) - 1)
print px

The noise in your sample produces results that vary by more than a whole sample, so I reduced it.  Fitting the curve using more of the points of the peak helps tighten the estimate somewhat, but I'm not sure if that's statistically valid, and it actually makes the estimate worse for the lower-noise situation.

With noise = 0.2 and 3-point fit, it gives values like 0.398 and 0.402 for offset = 0.4.

With noise = 2.0 and 13-point fit, it gives values like 0.156 and 0.595 for offset = 0.4.


MFCC window size at different sampling rates The general recommendation for window size when calculating MFCC seems to be 20-40 msec. This is most often recommended in a context of 16000 samples per second, so leading to a window containing 320-640 samples.

What if my input is at 8000 samples per second, should the time frame of the window be targeted at 20-40 msec or should the amount of samples be targeted at 320-640?

I understand the window needs to have a minimal duration in order to include lower frequencies and that the maximum duration is limited to have stable sounds (phones). What would be the reason to still prefer 40 msec over 20 msec?

What determines the minimal amount of samples? A 20msec window at 8000 samples per second is 160 samples, is that sufficient?

Context is speech recognition

I will answer your question from my experience. For most of the time, I've been using window length which is a power of 2, trying a few overlap percentages, training the system and picking the one which performs the best.

In the case of 16000 Hz, I would try 256 samples (16 ms) and 512 samples (32 ms).
For 8000 Hz, my bet would be on 256 (32 ms). Don't forget to try a few overlap values, however, it shouldn't change much.

At the end of the day, what matters is the number of MFCC filter banks and the resultant amount of coefficients.


Is this the correct solution for finding the samples of phase delay at a given frequency of an IIR BiQuad filter? I want to calculate the number of samples of phase delay at a given frequency of an IIR BiQuad Filter. I believe I have solved this, but I would like to confirm I have done it correctly. My work is as follows.

Background:

The transfer function of a BiQuad filter is:

$$
H(z) = g \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{a_0 + a_1z^{-1} + a_2z^{-2}}
$$
ω is normalized by the sample rate (samples/sec) with units of radians/sample. The normalized Nyquist frequency is π radians/sample and the normalized sample rate is 2π radians/sample.

$$
\omega = \frac{2 \pi f}{f_s}
$$
The frequency response of a system can be calculated by letting $z = e^{j\omega}$, where ω is a normalized frequency in the range (−π,π):

$$
H(e^{j\omega}) = g\frac{b_0 + b_1e^{-j\omega} + b_2e^{-j2\omega}}{a_0 + a_1e^{-j\omega} + a_2e^{-j2\omega}}
$$

The value of $H(e^{j\omega})$ will be complex. The phase response at the frequency ω is just the phase angle of the resulting complex number. The magnitude response at the same frequency is likewise equal to the magnitude of the number.

Euler's Formula for solving complex exponents is:
$$e^{ix}=cos(x)+isin(x)$$

Solution:
$$
H(e^{j\omega}) = g \frac{b_0 + b_1e^{-j\omega} + b_2e^{-j2\omega}}{a_0 + a_1e^{-j\omega} + a_2e^{-j2\omega}}
$$
$$
H(e^{j\omega}) = g \frac{b_0 + b_1(cos(-ω) + jsin(-ω)) + b_2(cos(-2ω) + jsin(-2ω))}{a_0 + a_1(cos(-ω) + jsin(-ω)) + a_2(cos(-2ω) + jsin(-2ω))}
$$
$$
H(e^{j\omega}) = g \frac{b_0 + b_1cos(ω) + b_1jsin(-ω) + b_2cos(2ω) + b_2jsin(-2ω)}{a_0 + a_1cos(ω) + a_1jsin(-ω)) + a_2cos(2ω) + a_2jsin(-2ω)}
$$
$$
H(e^{j\omega}) = g \frac{[b_0 + b_1cos(ω) +b_2cos(2ω)] + j[b_1sin(-ω) +  b_2sin(-2ω)]}{[a_0 + a_1cos(ω) + a_2cos(2ω)] + j[a_1sin(-ω))  + a_2sin(-2ω)]}
$$
$$
H(e^{j\omega}) = \frac{[g(b_0 + b_1cos(ω) +b_2cos(2ω))] + j[g(b_1sin(-ω) +  b_2sin(-2ω))]}{[a_0 + a_1cos(ω) + a_2cos(2ω)] + j[a_1sin(-ω))  + a_2sin(-2ω)]}
$$

Simplify:

$$m = g(b_0 + b_1cos(ω) +b_2cos(2ω))$$

$$n = g(b_1sin(-ω) +  b_2sin(-2ω))$$

$$p = a_0 + a_1cos(ω) + a_2cos(2ω)$$

$$q = a_1sin(-ω))  + a_2sin(-2ω)$$

$$H(e^{j\omega}) =\frac{m+in}{p+iq}$$

Find the Angle (delay in samples?):

For a given complex equation of:

$z= x+iy$

Where φ is the angle between the x axis and the vector z measured counterclockwise in radians:

$φ = atan2(y, x)$

Additionally, the phase of the ratio of two complex numbers is the difference of the phases of the numerator and denominator (numerator - denominator), as noted by Matt L in his answer below:

$$φ(ω) = atan2(n,m) - atan2(q,p)$$

To get the final phase delay in samples, we must take the negative of the phase we have solved and then divide that by the normalized frequency:

$$\tau_p(\omega)=-\frac{φ(\omega)}{\omega}$$

Therefore the final equation for phase delay for any normalized angular frequency ω would be:

$$\tau_p(\omega)= \frac{-atan2(n,m) + atan2(q,p)}{\omega}$$

(Updated from Matt L's feedback.)

You almost got it, but the final result is wrong because the phase of the ratio of two complex numbers is the difference of the phases of the numerator and denominator, respectively, not the ratio of the phases:

$$z=\frac{u}{v}=\frac{|u|e^{j\phi_u}}{|v|e^{j\phi_v}}=\frac{|u|}{|v|}e^{j(\phi_u-\phi_v)}\Longrightarrow \phi_z=\phi_u-\phi_v$$

Furthermore, for finding the phase delay, you need to divide the phase by the (normalized) frequency in radians (as shown in this answer). The result is the phase delay in samples.


What are the keywords to find an algorithm that detects embossed text from a 3d model? Currently I am doing an acadmic research about an idea which came into my mind lately.  I have a 3d model and I want to extract embossed text .. anyone knows any papers? or if it is already proposed? . I tried to research different keywords but nothing popped up which is relevent.

The text are bulging out (by a few millimeters)

One approach would be to use 3D reconstruction and then use deformable ICP (iterative closest point) to recover the deformation map of the embossed text. The deformation map would transform the CAD model to the deformed scene. The inliers you compute during correspondence estimation would also allow you to segment these characters. (This might be similar to computing optical flow)

For reconstructing your scene, I would recommend photometric stereo techniques. They are pretty accurate in reconstructing small details. For large objects, they would not work. Some intro texts are here and here.

Also, in such cases, you could use an appropriate lighting to emphasize the text. A low angle ring light would be an example.


Have you seen these kind of DFT artifacts? we've been working with the Spiral DFT Implementation on an FPGA and managed to get it working.
Unfortunately we see these strange artifacts and have no idea where they come from. Our input is a perfect sine wave, generated by a python script and fed into the simulation as a text file.

We have 2 issues:

Additional peaks in the frequency spectrum that should not be there.
The main amplitude breaks down in the center of the peak.

--> we thought this is due to a value overflow but the phenomenon persists even if we decrease the amplitude of the input signal.

We should note that our DFT takes 16 bit values, which are the upper 16 bit of our 24 bit sample values.

Here are some examples:

Different Frequencies:

Frequency: 10kHz, Amplitude: (2**23-1)/14

20kHz, Amplitude: (2**23-1)/14

Frequency: 30kHz, Amplitude: (2**23-1)/14

Changing the Amplitude results in following:

Frequency: 6.5kHz, Amplitude: (2**23-1)/10

Frequency: 6.5kHz, Amplitude: (2**23-1)/13

Our Sample Rate is 65.1 kHz.

If you have seen this phenomenon before or have suggestions on how to approach this issue, we're happy about any ideas! :)

I think your sine have some sample(s) missing or is poorly sampled. Do you have a way to check that? You said that you use 65.1k to sample the sine, but what is its frequency?


How to find roll-off from its transfer function equation? I want to find roll-off (unit dB/decade) of highpass filter $H(z) = 1 - a \cdot z^{-1}$.

I derived magnitude response to be

$1+a^2 - 2a\cdot \cos(\omega)$

To find out roll-off factor from above equation, how should I proceed?

I thought of doing

$1+a^2 - 2a\cdot \cos(10\omega) \over {1+a^2 - 2a\cdot \cos(\omega)}$

But $\omega$ is a variable, how do I simplify this?

Your method will work, but you need to first find the cutoff frequency by setting the magnitude response to $\sqrt {2} /2$. Then plug that in for $\omega$.


Matched filter for recieved RF digital packet by LPF Im transmitting by cc1350 launchpad microprocessor digital RF Packets, transmitted modulation FSK. bit rate of transmitting is  50Kbouad, frequency of transmitting is 868MHZ . Im sampling the signal(RF PACKETS) that Im transmitting by a dongle SDR-RTL and it samples on frequency 2048000Hz , I want to decode the RF packets that Im transmitting, so I build demodulator and I see the data of my RF packet which it's: (I plot it in matlab plot)

So if you see here on the plot it shows data of my sampled RF packet, and in order to decode it to 0/1 (binary), I want to do LPF in order to smooth/smear much more my received Rf packet "to look like smoothly sinusoidal"
any help how can I do that?
I mean by demodulator is that I've done two matched filters (Mf1, MF2) as I shown in the photo below:

Well, you know that the output of your filters shouldn't change faster than a symbol period. You could low-pass filter these to a bandwidth only a bit bigger than the symbol rate.
Since a filter is a linear operation, you can commute it and have it before the matched filters, so you only need one low-pass filter.

smoothly and to do zero crossings in order to choose when I have 0 and 1 in binary

Don't. You should be using a threshold instead of going by zero crossings: the zero is your maximum likelihood decision boundary.
You'll need symbol timing recovery, anyways. Since this is a bursty system, not an "endless" stream of symbols, I'd guess the packet preamble is probably meant to give you that information – use it!
If you can correctly estimate symbol timing from the preamble, your Matched Filters can be coherently decimating down to the symbol rate.


Algorithm for Discrete Signal Decomposition I am building a sensor and I am trying to understand how to process the signal that it generates. The sensor has a library of reference signals. When 'event a' occurs, it produces signal A, when 'event B' occurs it produces signal B. EX:

It is possible for multiple (but few) events to occur, and their signals would add together linearly. I am looking for an algorithm for signal decomposition which takes advantage of the fact that the decomposition is discrete (It can give integer multiples of the library signals only). I also need some measure of the reliability of the identification.

Ex:

  Input:

     (I made this image by taking A + 2 * B + noise)

Example Output:
A + 2 B  -  95% confidence
A + B    -  3% confidence
A        -  1% confidence
None of the Above - 1% confidence

I've noticed that I can compare two different possibilities by taking an inner product of the normalized proposed solution and input signals. For instance:

Input Signal .  (A + 2 B) = 0.999

but

Input Signal .  (A + B) = 0.985

So the first solution is a closer match to the input signal. I thought about using a Newton-Raphson to maximize this inner product, but when the actual sensor is in use it could have 1000s of signals in it's library and that would not be realistic.

I'd just take a KIS (keep it simple) approach as a first step.

Define your (unknown) signal as
$$ s(t) = k A(t) + m B(t) + n(t)$$
and then just define the error:
$$ e(\tilde{k}, \tilde{m}) = \sum_{\forall t} \left| s(t) - \tilde{k} A(t) - \tilde{m} B(t) \right|^2 $$
and then your estimates $\hat{k}$ and $\hat{m}$ are just chosen as:
$$
(\hat{k}, \hat{m}) = \arg \min e(\tilde{k}, \tilde{m})
$$

The R code below does this, and the output seems to do the right thing. Here $n(t)$ is a Gaussian noise with zero mean and standard deviation of $0.1$.

The actual $k$ and $m$ values are at the green circle, the estimates are the blue diagonal bar beneath it (in the rightmost plot).

This appears to work for any values of $k$ and $m$, including zero.  You can probably extend it for more signals $C(t)$, $D(t)$, etc. too.

R Code Below

# 30067

prototype_length <- 40

A <- rep(0, prototype_length)
A[(prototype_length/4+1):(3*prototype_length/4)] <- sin(2*pi*seq(0,prototype_length/2-1)/prototype_length)

B <- rep(0, prototype_length)
B[(prototype_length/2+1):(prototype_length)] <- `^`(sin(2*pi*seq(0,prototype_length/2-1)/prototype_length),12)*3/8

par(mfrow = c(1,3), pty="s")
plot(A, type="l", col="blue", lwd=5)
lines(B, col="red", lwd=5)
title('Original Signals')

k <- sample(0:10, 1)
m <- sample(0:10, 1)

data <- k*A + m*B + 0.1 * rnorm(prototype_length, 0, 0.1)
plot(data, type="l")
title('Scaled & Added Signals')

err <- array(0,c(11,11))
kvals <- array(0,c(11,11))
mvals <- array(0,c(11,11))
for (khat in 0:10)
{
  for (mhat in 0:10)
  {
    err[khat+1,mhat+1] <- sum('^'(abs(data - khat*A - mhat*B), 2))
    kvals[khat+1,mhat+1] <- khat
    mvals[khat+1,mhat+1] <- mhat
  }
}

image(0:10,0:10,err, xlab="k hat", ylab = "m hat")
ix <- which.min(err)
lines(kvals[ix] + c(-0.2,0.2), mvals[ix]+ c(-0.2,0.2), col="blue", lwd = 10)
points(k, m, col="green", lwd = 10)
title('Error, actual values, and minimum')


Concept question - What is noise? I've been studying image processing methods (without signal processing background) and many things are interesting to me. I was think about noise. In many definitions noise is unwanted  and without meaning  electrical or electromagnetic energy that degrades the quality of signals and data. Considering noise a signal type. So my question is

Noise is undesirable signal, or incomprehensible signal?

Note: Sorry by the philosophical question.

Noise is definetly undesired, but not incomprehensible.

The fact that noise is comprehensible is the reason we talk about it. Cosmic Rays might be present, but they have no effect on an old mechanical clock. Would you consider them "noise" when you read the clock?

You might call that loud disturbing sound that you hear "noise". A deaf person would not consider that noise, because it's incomprehensible.


2D IIR filter design (x-axis width proportional to y-axis) I want to design a 2D IIR filter (that is, created line-by-line progressively, in this case from left-to-right).  What is the best way to go about this?

I would like the response to look something vaguely like this:

The actual shape in the x-axis (time) isn't important - what's important is that the response's width for the x-axis is proportional to the y-axis position.

Any ideas?

First, you might be looking for a FIR filter (2d convolution)? If that's the case, you gave the answer yourself.

If you're actually looking for an IIR filter, then I'd like to ask for the desired order.


Zero-pad before or after windowing for FFT What's the correct way. Should I zero-pad a signal before or after applying a windowing function?

If done correctly, the order should not matter. You apply windowing to remove the jumps at the end of the signal segment by fading the values there to zero. Zero-padding then continues with zeros. Irregardless if padded or not, the windowing should only be applied to the original signal segment.

The wrong way would be to zero-pad and then apply windowing to the now longer signal over its full length, i.e., without regard to the position of the jumps. This would not remove the jumps and thus keep the background in the spectrum that is not related to the frequency content of the signal but the segmentation method.


how do windowing functions not distort the signal? I'm trying to gain a better understanding of windowing functions, but so far I can't tell why they're used. As an example, I'm familiar with the Hamming Window:

$ w(k) = 0.54-0.46cos( \frac{2{\pi}k}{N-1})$ if $0<k<N$

But I don't understand how this can be done without distorting the signal. Are these specially chosen functions that just happen to retain the spectrum of the signal sample? Can it be proven mathematically that multiplying by this function retains the spectrum? And what is the advantage to using this, rather than simply windowing it and not multiplying it by any windowing function at all (i.e. a windowing function of 1)?

They do "distort" the signal, as you noted. Multiplying by a window function does not "retain" the spectrum of the original signal. Recall the multiplication property of the Fourier transform:

$$
x(t) w(t) \Leftrightarrow X(\omega) * W(\omega)
$$

where the $*$ denotes convolution. So, multiplying by a window function in the time domain yields convolution of the input signal and window function's spectra in the frequency domain. In this way, there is some "distortion" of the input signal's spectrum.

As you noted, however, there are cases where using window functions can be useful. There is a good illustration of the tradeoffs here. I will shamelessly steal the plot for reproduction here:

Summarizing the discussion briefly:

If you don't use a window, that is equivalent to using a rectangular window, assuming you are analyzing a finite-length signal (and if you aren't, then you better be patient). The rectangular window has a very slow decay rate on its sidelobes. You have to move away from the main lobe a long way in frequency in order to achieve a certain amount of attenuation.
Say you had a case where you wanted to analyze two tones in a signal, where one is 40 dB stronger than the other. Compare this to the plot of the rectangular window's frequency response in the linked answer above. You can see that for that case, the sidelobe response of the rectangular window never drops below ~30-35 dB or so (except for the very narrow nulls). That means that the sidelobe response from the stronger signal would "wash out" the other peak, thwarting your attempt to analyze both signals simultaneously.
The solution to such a problem is to use a window function: there are a few illustrated there, but there are a lot of them. Select one that has enough sidelobe attenuation for your application. For example, the Blackman-Harris window has greater than 90 dB of sidelobe attenuation, which provides a lot more dynamic range than the rectangular window.
You pay a price for such a choice, though. The "distortion" that you referred to, caused by the convolution of the input signal and window spectra, comes into play. Note that the greater the dynamic range of the window function (i.e. the lower its sidelobes are), the wider its main lobe is. This can affect your ability to resolve two tones that are very close to one another (if they are too close, then the width of the window's main lobe will cause them to smear together).

This is the inherent tradeoff when selecting a window function.


Nonuniform Background Subtraction I have multiple images like that shown below in which there are essentially two different brightness levels. I've looked into histogram equalization techniques in scikit-image as well as several background subtraction techniques in ImageJ, neither of which seem to be effective over the parameter ranges I've tested. I've thought of using edge detection to process the two regions separately, but my concern is that it will result in unnatural edge effects at the boundaries. There's also a lot of features in the image, making the simple edge detection techniques I've tried so far not work so well.

Are there better techniques for doing this sort of processing? I'm working in python with scikit-image, but would welcome any solution.

If it's not clear what kind of result I'm looking for, consider the below image that has an evenly balanced background.

Histogram equalisation (on its own) is not likely to work here because it corrects the brightness levels of the image as a whole.

What you are dealing with here is a uniform reduction in brightness akin to the effect of an optical filter.

What is really required is the correction of the mean brightness level over the "patch" that appears darker. Therefore, the problem now is to find the patch and correct it.

I am going to use Octave for the rest of this post but as you will see too it is not too difficult to transpose the idea to Python.

Let's have a look at the damage:

%Load the image and plot its average values by column.
a = imread('SomeImage.jpg');
imshow(a); %This should simply depict the image
%Let's take the cross section of the image by "squashing" all of its columns into their average.
q = mean(a);
plot(q);

This comes out like:

The most striking feature in this image is the sharp drop of the brightness level, about a third of the image width from the left.

Let's now look at the vertical cross section by "squashing" all of the image's rows into a mean column. All you have to add is:

q2=mean(a');
plot(q2);

And this comes out like:

In fact, if you have MATLAB somewhere handy, you might want to try to surf(a) and that will show you the "basin(s)" very clearly. For some reason, Octave's surf is getting its triangles tangled when I try to do it.

Right, so, what do we observe from this?

The image on either side of the discontinuities is predictable. It looks as if you simply took part of the waveform and slide it downwards...similar to the effect of a neutral density filter. There's simply a dip in the brightness that has a rectangular shape.
The problem on the first image is worse than the problem on the second image. On the second image we observe a sharp discontinuity which is approximately 4 pixels wide but it is very very sharp. On the first image, the discontinuity is about 8 pixels wide. Whatever data there was in those pixels has been "swallowed" by the gradient.
Working in 1 dimension (at first) you can see that just by either adding a value to the dip or subtracting a value from the "step", you can more or less match the means of the two regions on either side of the discontinuity. But for those pixels in the transition band, you are going to have to do something else (e.g. interpolation) to remedy the problem. The data in that region has been lost.

Where is the (big) patch? It starts approximately at 85,78 and extends all the way to the margins of the image. The X point is found by squashing all columns into one row and then looking at where the discontinuity starts and the Y point is found by doing the same but this time squashing all rows into one column.

Let's try to correct the X direction. First of all, we are going to get rid of the transition gradient. This is between pixels 85-93 (not across the whole image though, remember that this was a patch and that here we are trying to correct for the one dimensional case).

Q = q;
Q(85:93) = 0;
plot(Q);

Now, correct for the brightness. The difference is approximately 36. All I am doing here is to add a "step" (a constant value) from pixel 93 onwards, like so:

m=93;
Q2 = Q+36*[zeros(1,m),ones(1,length(Q)-m)];
plot(Q2);

This now looks like this:

What can you do for that gap of 8 values? You can fill it in with a constant value but that is going to show up as a "band" of 8 pixels and it will definitely be visible. The other thing you can do is to fill it with some Gaussian noise that has approximately pixel value 36 mean and enough variance so that you cannot distinguish it from the surrounding pixels. Unfortunately, because in this case the noise would be uncorrelated that would show up as banding again but this time much more "fuzzy". And finally, the other thing you could do is to train a model that would "generate" realistic values from the surrounding areas so that you "heal" the bands gracefully.

This is the basic idea. To automatically find the gradients you can either extract the gradient of the image with imgradient or do a similar calculation in a loop with a threshold on the local slope.

Your problems are not over here. As you can see in all images, there are more than 1 "basins" in the image. Your first image has a light band on the top right part which goes right into the darkness right below it.

Hopefully they are all rectangular so you can simply repeat the exact same process as many times as it is required until the image contains "no basins".

Hope this helps.


Comparison between watershed and grabcut I know that watershed and grabcut algorithms are both common tools for image segmentation. They both seem to rely on predefined markers to form something like a source-sink system and then uses the gradient of the image to find the boundaries.

The largest difference I observe is that grabcut uses a graph structure where the gradient translates to edge weights, while watershed intuitively works by slowly adding water to the gradient image. But to me, it seems like they would yield similar results?

My question is what are the relative strengths and weaknesses between these two algorithms? When should I use one over another?

Watershed-based segmentation will typically lead to over-segmentation, and is very sensitive to local image noise (e.g. see discussion here). Typically watersheds are pruned/merged by using thresholds on size (area/relief/volume).

I have not had much direct experience of GrabCut, mostly know it from reading  papers. However I can say that it is a more "modern" approach than watersheds (i.e. much more recently developed). Because it is based on global optimization (via graph cuts), it is less sensitive to local image noise, compared to watersheds.


Implementing Kalman filter or extended or unscented with only position information I am new to Kalman filter. Is it possible to use a Kalman filter with just the position information of the target?

I am so confused with states and measurement.

from the generality of your question then yes you can design a Kalman filter which would accept the target position as the only measurement possibly corrupted with noise.

Then the Kalman filter will try to estimate the true position of the target based on the noisy position measurement and the assumed motion model. Hence in effect the Kalman filter would be reducing the measurement noise of the target position measuring device.

The states of the filter would be defined based on the assumed motion model: a constant acceleration model would use for example the position, velocity and acceleration as the internal states which are being estimated (when the acceleration is unknown) or another model would use a deterministic input for the known acceleration and use two states instead.


What is Earth Movers Distance and how is it used in pattern recognition? I have obtained this from Wikipedia about the application of Earth Movers Distance to pattern recognition. Can anybody explain the concept of signature used in this Wikipedia article. Thanks.

A set of clusters. Each cluster is represented by a pair with (1) one single point (usually the center in the cluster) ,and (2) a weight coefficient that usually in proportional to the size of the corresponding cluster. Note that the point can be anything including the coordinates of the pixel, the color information, intensity...etc.


Image processing using Gabor filter I need to compare images using Gabor filter. I cant seem to understand how to find a quantifiable value to compare through Gabor.
Can someone give me steps to do the same. I got some code, but it doesnt compare two images.
Any help is appreciated.

Gabor can be used as a filter bank in order to extract texutre
features from an image (or image patch).
To compare images (or patches from image) using gabor I would follow the following pipline:

compute a filter bank from the Gabor (different scale and orientation)

post process to create new feature. you can compute the energy, and first statistics measures on the filtered images.
now that you have a vector of features for each image you need some measure to compare tham. You can use simple L2, or fisher criterior as some papers offer Comparison of Texture Features
Based on Gabor Filters

Please also take a look at the following question that is simmilar to yours.


Periodogram and Welch periodogram comparison I'm trying to compute the PSD using the periodogram method. I have a signal that is periodical at 5kHz. The frequency resolution sought for is 1.221Hz. It's a PRBS signal that has its energy drop to zero at the generation frequency, i.e. 5kHz in this case. The signals i prety much noise clean.

I tried computing the periodogram and I get what I'm after. Then I tried using Welch's averaging method with N = 500 segment length and overlap of 50% and I get something completely different. I would expect to see an exact copy of the periodogram since there is practicaly no noise. Did I completely misunderstood the method?

What I obtain is this:

I used the following lines of code to estimate the periodogram exactly at the frequencies I want and I scaled the periodogram. ( I have also added .mat files if anyone wants to give it a go)

https://www.dropbox.com/s/pe6qbi883kkocxf/data.mat?dl=0
https://www.dropbox.com/s/2fqmqrg1xpgbibi/freqs.mat?dl=0

load('data.mat');
load('freqs.mat');
1/DT = 5000;
[pxxp,fp] = periodogram(V,rectwin(length(V)),freq,1/DT);
[pxxw,fw] = pwelch(V,rectwin(500),250,freq,1/DT);

semilogx(fp,20*log10(pxxp*8775))
hold on
semilogx(fw,20*log10(pxxw))

Looks similar to me

>  clear all
>     close all
>     load data
>     samplerate=1
>     periodogram (V)
>     title("Periodogram",'FontSize', 14)
>     xlabel('Frequency (Hz)','FontSize', 14)
>     ylabel('Amplitude (dB)','FontSize', 14)
>
>     figure
>     pwelch(V)
>     title("Welch's method",'FontSize', 14)
>     xlabel('Frequency (Hz)','FontSize', 14)
>     ylabel('Amplitude (dB)','FontSize', 14)


Any guidance on converting Time Domain to Frequency Domain in Matlab Below is the graph of my load cell readings at 3 samples per second. I have collected 500 such samples. I want to stabilize the readings and read that first approach is to convert time domain data to frequency domain. I am new to DSP, so i seek your expert guidance on this.
My main aim is to stabilize the readings so that i can increase an accuracy of my measuring device.
Graph after fft is as shown below

@Nimit. I recommend that you study and learn to design and implement what are called 2nd-order IIR lowpass filters. Once you thoroughly understand such filters you should design a 2nd-order IIR lowpass filter, that has unity gain at zero Hz, and experiment by applying your noisy signal to that filter and see how much noise reduction is achieved. Then you should apply your noisy signal to a cascade of two identical 2nd-order IIR lowpass unity-gain filters and see how much noise reduction is achieved. Finally you should apply your noisy signal to a cascade of three identical 2nd-order IIR lowpass unity-gain filters and see how much noise reduction is achieved. Nimit, you have much experimenting to do.


How to extract prosodic cues from a wav file using Python How do I extract the following features as from wav files in Python -

Voicing rate (Voice rate) - frequency of phonemes while speaking, and represents the pace of a conversation
Number of autocorrelation peaks
Pitch, pitch bandwidth

This project Speech Signal Processing Toolkit (SPTK) provides several features you are looking for. Here is a good wrapper around it https://github.com/r9y9/pysptk
This wrapper is using a slightly different version of SPTK, but installation is straightforward.

pip install pysptk

1) I would use kaldi or another speech recognition package on phones output and a VAD system to compute the voicing rate.

2) For autocorrelation, convert the waveform to mel-cepstrum with the function mcep then compute the autocorrelation with c2acr

3) Here is an example to compute the Pitch http://nbviewer.jupyter.org/github/r9y9/pysptk/blob/master/examples/pysptk%20introduction.ipynb


How to identify and differentiate frequency and time in EEG data using python? I have an EEG data of 200 Hz and sampled at 4097. I have few doubts and questions related to this data:

What does it mean by sampling in EEG data?
I am getting the sample rate and the frequency when I plot the data as shown below. Shouldn't it be the time and frequency?

I don't understand what I'm missing here. Can anyone explain how am I supposed to identify the time and frequency of EEG data?

EEG data file : http://filesave.me/file/123174/eeg-txt.html

Thank you in advance.

Code for the above plot:

x = np.array([])
input = open('eeg.txt', 'r')
for file in input:
    file = file.replace(',','')
    x = np.append(x, float(file))
input.close()
plt.subplot(1, 1, 1)
plt.plot(x)
plt.ylabel('EEG signal')
plt.show()

FFT calculation:

import numpy as np
x = np.array([])
input = open('eeg.txt', 'r')
for file in input:
    file = file.replace(',','')
    x = np.append(x, float(file))
input.close()
plt.subplot(1, 1, 1)
ff = fft(x)
plt.plot(ff)
plt.ylabel('EEG signal')
plt.show()

Wavelet decomposition:

coeffs = wavedec(x, 'db4', level=6)
cA2, cD1, cD2,cD3,cD4,cD5,cD6 = coeffs

plt.subplot(7, 1, 1)
plt.plot(x)
plt.ylabel('Noisy Signal')
plt.subplot(7, 1, 2)
plt.plot(cD6)
plt.ylabel('noisy')
plt.subplot(7,1,3)
plt.plot(cD5)
plt.ylabel("gamma")
plt.subplot(7,1,4)
plt.plot(cD4)
plt.ylabel("beta")
plt.subplot(7,1,5)
plt.plot(cD3)
plt.ylabel("alpha")
plt.subplot(7,1,6)
plt.plot(cD2)
plt.ylabel("theta")
plt.subplot(7,1,7)
plt.plot(cD1)
plt.ylabel("delta")
plt.draw()
plt.show()

Results of wavelet:

It is hard to answer your question, since you do not seem to have experience with EEG data and/or general signal processing.

You say, your data is sampled at 200 Hz, which seems good to me for EEG data.

What does it mean by sampling in EEG data?
That means, that the (continous) EEG signal is choped up into discrete values. This is always necessary for further signal processing with computers.
I am getting the sample rate and the frequency when I plot the data as shown below. Shouldn't it be the time and frequency?

I don't know what you mean by I am getting the sample rate and the frequency when I plot the data. You certainly don't. The sampling rate/frequency is a fixed number given by your experimental hardware. You don't get it from your data - you have to know it from your experiment to get meaningful results from your data.

I also don't know what Shouldn't it be the time and frequency? means. Your EEG data is certainly a voltage signal over time, hence you have a signal in time domain.

Your sampling rate is given from your experimental settings and is 200 Hz as you say. Hence, each time point is separated by $$\Delta t = \frac{1}{f} = \frac{1}{200\,\textrm{Hz}} = 5\,\textrm{ms}$$

Your total signal that you have plotted is therefore $$4097 \times 5\times 10^{-3}\textrm{s} = 20.485\,\textrm{s}$$
long.

That is what you know from your data. Any further data extraction is up to you or your research question. If you'd like to get the frequencies that are in the signal, try to Fourier-transform it and extract the frequency region that you are interested in. As I see it, you have some apparent frequency components at approx. 13.4 Hz, 57.7 Hz, 81.0 Hz, and 88.0 Hz. Probably you are only interested in the lower spectrum around 13.4 Hz...


Getting hilbert transform of IQ signal Im new to communications so this may be a stupid question, but I am somewhat confused. I am trying to train an automatic modulation classification (AMC) algorithm on IQ data from the RadioML dataset. The data comes  in IQ format (eg [1024, 2] matrix).

I am reading papers now which use feature based methods and talk about getting the hilbert transform / DFT of the received signal and extracting features from there (for example, the std deviation of the instantaneous amplitude, see Table 1 here).

My question is what do I take the hilbert transform/DFT of? Would it be the IQ signal (I+jQ) or is there something more else involved (multiplying by sin/cos)?

Usually an IQ signal is generated from a real (RF, AF, etc.) signal by a radio (by quadrature heterodyne, Tayloe demodulation, complex vector arithmetic, or otherwise).  However, if you need IQ data, and the radio only generates I (or real) component samples, a Hilbert transform can be used to generate “fake” Q components (an approximation by making some assumptions about the original signal) to use in downstream IQ processing or analysis.
Thus, if you have an IQ signal, there is rarely any need to generate another fake similar Q signal component. Just use the Q component the radio provides instead of a Hilbert transform of the I component.  An IQ signal created by a Hilbert transform of the I component contains at most half the amount of information as contained by the original IQ signal (e.g. potentially aliasing the upper and lower sidebands).
And, to get a real signal back, just reverse whatever process (superhet, etc.) the radio used to generate IQ from RF.


Can Aliasing Happen in Laplacian Pyramids? I want to use the Laplacian pyramid framework for representing images. The The Laplacian Pyramid as a Compact Image Code's language seems to provide freedom in the choice of blurring/interpolation kernel, subject to certain constraints. However, as we are performing subsampling when constructing pyramids, would aliasing be a problem with an injudicious choice?

I read in the Toronto lecture notes on image pyramids(Slide 8) that a laplacian pyramid can reconstruct an image exactly, agnostic of the kernel. Can i take this to be a guarantee against aliasing, or can aliasing be "hidden" by the reconstruction process?

The reconstruction will remove aliasing, no matter which filter kernel is used.

The lower layer of the 2-layer Laplacian pyramid on page 8 of the Toronto lecture notes, adapted here:

Figure 1. A 2-layer Laplacian pyramid.

can be redrawn as:

Figure 2. An alternative presentation of the lower layer of the Laplacian pyramid of fig. 1.

This layer is an identity operation: what is subtracted is also added. This allows to draw the next (and in this case the only) layer above it the same way and it also can be identified to be an identity operation. And so on until the whole process becomes an identity operation.


How would Fourier and Cosine Transforms responds to summation of cosines with same frequency but different phases? For example, if I have two signals, $\cos(2\pi ft+\frac\pi4)+\cos(2\pi ft+\frac\pi3)$, what would be different in both transforms (Fourier and cosine) how would the spectrum changes?

And What would happened for a signal like $\cos(2\pi ft)+\cos(2\pi ft+0.99\pi)$.

These integral transforms are linear, which means the transform of a sum is identical to the sum of the individual transforms, i.e.

$\mathcal F\left\{\cos(a)+\cos(b)\right\}= \mathcal F\left\{\cos(a)\right\}+\mathcal F\left\{\cos(b)\right\}$,

and also, that when multiplied with a constant, that constant can be "inside" or "outside" the transform:

$\mathcal F\left\{\alpha g(t) \right\}=\alpha\mathcal F\left\{g(t) \right\}\, \forall \alpha \text{ const.}$

No magic happening here; since (Euler's formula! It's important!):

$$\begin{align}
\mathcal F\left\{\cos(a)\right\}&=\mathcal F\left\{\frac 12\left(e^{i a}+e^{-ia}\right)\right\}\\
&=\frac12\mathcal F\left\{e^{i a}+e^{-ia}\right\}\quad\quad\quad\quad\quad||\,a=2\pi f t + \varphi\\
&=\frac12\mathcal F\left\{e^{i(2\pi ft + \varphi)}+e^{-i(2\pi ft + \varphi)}\right\}\\
&=\frac12\mathcal F\left\{e^{i 2\pi ft}\cdot e^{i\varphi}+e^{-i 2\pi ft }\cdot e^{-i\varphi}\right\}\\
&=\frac12\mathcal F\left\{e^{i 2\pi ft}\cdot e^{i\varphi}\right\}+\frac12\mathcal F\left\{e^{-i 2\pi ft }\cdot e^{-i\varphi}\right\}\\
&=\frac{e^{i\varphi}}2\mathcal F\left\{e^{i 2\pi ft}\right\}+\frac{e^{-i\varphi}}2\mathcal F\left\{e^{-i 2\pi ft }\right\}\tag{*}\\
\end{align}$$

we see that for any fixed phase $\varphi$ of a cosine, for example $\frac\pi 3$ or $\frac\pi 4$, the Fourier Transform "consists" of two dirac impulses $\delta$ with positions that only depend on the cosine's frequency $f$ and with amplitudes that depend only on the phase offset.

Now, just insert your different $\varphi$ in $(*)$ and add up those sums of Fourier Transforms to see the very basic point:

What you're modelling here, mathematically, is the interference of two waves at the same frequency. Tada! Physical significance added.

Now: good luck with the rest of your homework.


Impulse Response from Frequency Response at Specific Values I'm at a loss trying to answer the following question:

The frequency response $H(j\omega)$ of a length 4 FIR filter with real impulse response has the following specific values: $H(j\cdot0) = 2$, $H(j\cdot \pi/2 ) = 7 − 3j$, and $H(j\cdot\pi) = 0$. Determine the impulse response function $h(n)$.

I can't think of anything to do but apply the inverse DTFT but how would that work at specific values?  Is there another method I'm forgetting?


Reason for >8 channels in an audio codec I'm new to audio (and signal processing in general), and I've been learning about the AAC codec. I've read that it supports 1-48 synchronous channels in a single stream. My question is: what need is there for 48 channels in a single audio stream? I understand that if you have a 7.1 speaker setup you can send a separate channel to each speaker, but what is the need for 48 separate channels?

Since it will most likely not hurt to have the possibility for more streams, this might for one be just to future-proof the system. Nothing would be more stupid than to specify for 8 channels and suddenly 9.1-systems become the next big thing and your codec is not suitable.

Other than that I can imagine bigger scenarios with much more channels such as recording all the voices in an orchestra. But I honestly don't know whether or not AAC is of any importance there.


Is there any relationship between the convolution and cross-correlation of two signals? There are two square 2-D signals with the same size $m \times m$, if the result of convolution of those two signals is first picture, And the result of cross-correlation of those two signals is second picture.How to prove it in mathematics?

okay, let's associate the symbols $u$ and $v$ with the dimensions of $x$ and $y$ with respect to each.

and let's assume periodicity in both the $x$ and $y$ dimensions for both pictures. $f(x,y)$ and $g(x,y)$.

$$ f(x+m,y) = f(x,y) \qquad \forall x,y $$
$$ f(x,y+m) = f(x,y) \qquad \forall x,y $$
$$ g(x+m,y) = g(x,y) \qquad \forall x,y $$
$$ g(x,y+m) = g(x,y) \qquad \forall x,y $$

then 2-D Cross-Correlation is:

$$ R_{f,g}(u,v) = \sum\limits_{x=0}^{m-1} \sum\limits_{y=0}^{m-1} f(x,y) \cdot g(x+u,y+v) $$

and 2-D circular convolution is:

$$ (f\circledast \hat{g})(u,v) = \sum\limits_{x=0}^{m-1} \sum\limits_{y=0}^{m-1} f(x,y) \cdot \hat{g}(u-x,v-y) $$

if $\hat{g}(x,y) = g(-x,-y)$, then the cross-correlation and the circular convolution are the same.  if $\hat{g}$ is an upside-down copy of $g$, the the correlation and convolution are the same.

or maybe they are upside-down of each other.  i think it's:

if $\hat{g}(x,y) = g(-x,-y)$, then

$$ (f\circledast \hat{g})(u,v) = R_{f,g}(-u,-v) $$

that looks about right.  someone needs to check this.


How to recognize an object from a small training set of images? Given a small training set of images (say, around 6 max) of the same object, how to measure how likely is another query image to contain the given object?

The training images are "clean", i.e. only the original object is present. On the other hand, the query images normally contain only an object of the same class (i.e. if the training images contain a mango, the query images will be definitely fruits, although maybe not mangoes). The query images have varied backgrounds, though.

I tried SURF and FLANN (from OpenCV), with various values of Hessian, but it's not very accurate: it works for certain images but not for others. Maybe it's because of the nature of the objects I'm training to detect? Any pointers would be appreciated. Thanks!

Few reasons i could think of are:

Size of your training set is very small. Larger training sets have always been the key for accuracy.
Each algorithm will have some drawback like SURF is not good at viewpoint change and illumination change. So if the test images vary from the training images in case of illumination and viewpoint, results on those data will be poor.I don't know how FLANN works.


Applying Image Filtering (Circular Convolution) in Frequency Domain To filter an image we can:

Use a 3x3, 5x5, 7x7, etc. filter, that is convolve the image and the filter in the space domain.
Use a FFT on both the image and the filter, multiply them together in the frequency domain, and then perform an IFFT on the image to get the filtered result.

These two approaches should produce equivalent result.

Questions:

Are there any conditions placed on the filter so that the results are equal?
Will there be any information loss in the second approach or will the images be exactly equal (neglecting floating point precision error)?
As the kernel size increases, will the second approach start outperforming the first?

In my StackExchange Signal Processing Q38542 GitHub Repository you will be able to see a code which implements 2D Circular Convolution both in Spatial and Frequency Domain.
Pay attention to the function CircularExtension2D().
This function align the axis origin between the image and the kernel before working in the Frequency Domain.
Remember that for Discrete Signals the implicit assumption on signals, In frequency Domain analysis, is being periodic (Circular).
In the discrete case one could indeed apply Circular Convolution by element wise multiplication in the Frequency Domain.
With proper padding one could apply linear convolution using circular convolution hence Linear Convolution can also be achieved using multiplication in the Frequency Domain.
See:

In depth description can be found in FFT Based 2D Cyclic Convolution.
Regarding your questions:

The filter is just an array of numbers. As long as you are after 2D Circular Convolution there is no constraints on the Filter. If it is valid for 2D Spatial Circular Convolution it is valid for Frequency Domain Circular Convolution.
Up to Floating Point Quantization errors both are mathematically equivalent (See Convolution Theorem).
If the Convolution Kernel is similar in its size to the image and both are large enough Frequency Domain Convolution becomes more efficient than Spatial Domain.


Implementing echo on adsp 21262 I'm working on an echo implementation on the board adsp 21262. I've added a buffer to have more samples to work with. I end up then with 4 buffers of size 2048 each.

The problem is that I need to access at least 7000 samples back. (An echo can only be heard with a delay of 1/15th of a second, we have 2*48K~100k samples per second. So we need 100k/15~7k samples).

The problem is that 2048 samples per buffer is the max we can have because afterwards we're out of memory.
I'm trying to access the ROM but I have no idea as to how to do that. I've spent a lot of time looking for this information finds nothing.

Thank you !

Assuming you have this board, there are 512kByte of SRAM on the board. That should easily accommodate a 7k delay buffer.


Maximum likelihood estimator for multiplicative Gaussian noise So I'm trying to derive an analytical solution for a MLE that should estimate a static value polluted by multiplicative Gaussian noise.

The vector of measurements $\tilde{\boldsymbol{d}}$ is given as $a[n]x_u$ where $a[n]$ is the $n$-th realization of the random variable $A \propto \mathcal{N}(\mu_A, \sigma_A)$ and $x_u$ is an unknown constant that should be estimated from $\tilde{\boldsymbol{d}}$.

I got as far as the log-likelihood function, which I now need to maximize with regard to $x$. To do this, I need to take the derivative of the log-likelihood , set it to zero and solve for $x$. However, the summation term at the end is giving me a headache as I can't figure out its derivative.

$\ell(x, \tilde{\boldsymbol{d}}) = -N \cdot \log x - N\cdot \log\sigma_A - \frac{N}{2} \log 2\pi - \frac{1}{2\sigma^2_A}\sum\limits_{n=1}^N\bigg(\frac{\tilde{d}[n]}{x}-\mu_A\bigg)^2$

OK, let's have a look at one of the problematic terms:
$$
\frac{\delta}{\delta x} \bigg[ \bigg(\frac{\tilde{d}[n]}{x}-\mu_A\bigg)^2 \bigg ]
= - \frac{2 \tilde{d}[n] \bigg (\tilde{d}[n] - \mu_A x\bigg) }{x^3}
$$
which can be verified by Wolfram Alpha.

The full derivative of the summation term is then just this summed over $n$.


Finding phase angle via FFT I am using the FFT to calculate the amplitude and phase of AC cos signal. Regarding amplitude it gives me an accurate result but I get wrong result for a phase: 6 degrees instead of 36 degrees.

My MATLAB code is as follows:

fs=3840;                 % sampling frequency
N=64;                    % number of sample
f=60;                    % signal frequency in hz
w=2*pi*f;
amp=220;                 % ampitude of signal
delt=[1:N]*2*pi/N;       % sampling time interval
noise=randn(size(delt)); % added noise

sig=zeros(size(delt));
for i=1:1:length(delt)
    sig(i)=noise(i)+amp*cos(w*delt(i)+pi/5);
end

FF=fft(sig,N);
stem(abs(FF))
grid on

REAL=real(FF);
SUM_REAL=sum(REAL)
IMAG=abs(imag(FF));
SUM_IMAG=sum(IMAG)

% *******************************PHASOR REPRESENTATION*******************

real_amp=sqrt(SUM_REAL)                  % Amplitude of real part
imag_amp=sqrt(SUM_IMAG)                  % Amplitude of imaginary part
sumation=real_amp+imag_amp;
amplitude_nominal=sqrt(sumation)/N
phase=atan(abs(imag(FF))/abs(real(FF))) % phase in radians
PHI=abs(phase)*180/pi                   % phase angle in degrees

Since the length of your FFT is exactly one period of your sinusoidal input signal (N = Fs/F0) , the signal will appear in bin 1 (the one above DC) of the FFT result (the conjugate mirror is redundant for strictly real input).  To estimate the phase of the sinusoid of frequency F0, referenced to the start of the data window, just compute the atan2() of the signed imaginary and real components from that single complex bin 1.  No need to abs() components or sum or window or estimate the amplitude or examine any other bin.  Unless you want to sanity check the magnitude of bin 1 to make sure it isn't zero or just numerical noise.


Causality of Linear Phase Filter If a digital LTI filter has phase response on the form $arg[H(e^{jw})] = -\alpha\omega$, what could one say about it's causality (and consequently, the causality of it's impulse response)? Also is Linear Phase a characteristic of only Causal FIR Symmetric Systems?

For a generic value of the group delay $\alpha$, the filter has a 2-sided infinite impulse, and is therefore not causal.

Causal linear phase systems are (almost) always FIR systems with any of the 4 types of symmetry.

In these cases, the group delay is multiple of 1/2 (0, 1/2, 1, 3/2, etc.).

You can have causal linear phase IIR filters, but these are not of finite-order (cannot be implemented with a difference equation), so they are of little use in practice.


What is the purpose of the log when computing the MFCC? The steps of computing the Mel-Frequency Cepstrum Coefficients (MFCC) are:

Frame blocking -> Windowing-> abs(DFT) -> Mel filter bank-> Sum coefficients for each filter-> Logarithm -> DCT

But what is the purpose of the logarithm step?

The logarithm serves to transform a multiplication into an addition. It is part of the computation of the cepstrum. The basic idea is as follows:

Assume a source signal $x$ is convolved by some impulse response $h$. The resulting magnitude spectrum is

$$|Y(\omega)| = |X(\omega)||H(\omega)|$$

By applying the logarithm we get

$$\log |Y(\omega)| = \log |X(\omega)| + \log |H(\omega)|\tag{1}$$

If we want to equalize / undo the effect of filtering by $H(\omega)$ we can hope that this task is easier if we transform the convolution into additive noise. This is exactly what happens by taking the logarithm in (1).


Are there any analytic approaches to the design of quadratic polynomial filters? The following paper describes an application of the Teager-Kaiser energy operator to x-ray image enhancement:

  Reinhard Bernstein, Michael S. Moore and Sanjit K. Mitra, "Adjustable
  Quadratic Filters for Image Enhancement" Proc. IEEE International
  Conference on Image Processing (ICIP), Santa Barbara, CA, vol. 1, pp.
  287-290, Oct. 1997. http://vision.ece.ucsb.edu/publications/view_abstract.cgi?52

The authors develop intuition for the behavior of the filter through analogy with a similar linear operator (i.e. "Thus the output of a Teager filter is approximately equal to a highpass filter response weighted by the local mean."). For the sake of precision, by quadratic polynomial filters, I mean non-linear, non-recursive filters that can be completely characterized by a truncated Volterra Series, as follows (for the 1D case):

$$
y(n) = \sum_{m_1=0}^{N_1-1}{ h_1(m_1)x(n - m_1) } + \sum_{m_1=0}^{N_2-1}{\hphantom{.}\sum_{m_2=0}^{N_2-1}{ h_2(m_1,m_2)x(n - m_1)x(n - m_2) } }
$$

It seems that most approaches to the design of low-order polynomial filters involve system identification frameworks, but without any deep understanding as to why the estimated filters work.  Are analytic approaches based on linear analogies currently the state-of-the-art, or are there any known mathematical tools that can be used?

Not really an answer (hence this is community wiki), but I thought we should capture @Mohammad's code:

%Mohammad Z

%Teager-Kaiser Non-Linear Energy Operator.
function [out] = TKEO(x)
    N = length(x);
    x = x(:).';
    temp = x(2:N-1).^2 - x(3:N).*x(1:N-2);
    out = [temp(1) temp temp(end)];
end


For a discrete LTI system, does "bounded memory" imply "rational transfer function?" Every LTI system with a $\mathcal Z$-domain transfer function that is a rational function - aka a quotient of two polynomials in $z$ - can be implemented using a bounded amount of memory.

Is the converse true? That is, are there any LTI systems that can be implemented using finite memory that do not have a transfer function that is a rational function?

Proving rational transfer function $\rightarrow$ finite memory is easy - any such system can be rewritten in the time domain as a linear constant-coefficient difference (LCCD) equation, so you end up needing only a finite number of delay lines to provide an implementation of the system.

Proving finite memory $\rightarrow$ rational transfer function seems less easy.

One way to look at it is, does finite memory $\rightarrow$ LCCD equation? That is, are there any finite memory LTI systems that cannot be expressed as an LCCD equation?

Non-trivially, the answer is no, a Bounded In Memory (BIM) LTI system do not imply the system has a rational Transfer Function (TF).
Proof (TL;DR):
Every discrete LTI system $H$ has an Impulse Response (by definition), and hence have a TF in the z-domain (by construction), which could be infinite (IIR) or finite (FIR).
A FIR LTI is trivially BIM, requiring a bounded number of operation for implement every convolution.
A IIR LTI could have a Rational TF (RTF) -such as most ARMA systems- or could have an not rational -an infinite series (IS)- TF -such as the counter-example 1 below.
If the system has an IS TF which can be decomposed as a rational TF -such as the Counter-Example 2-, then the system is bounded in memory, because the rational TF can be decomposed into a pair of finite series, which can be implemented and evaluated with finite memory and finite CPU.
Similarly, if the system has an IS TF with an Known algebraic expansion (KIS), -such as the Counter-Example 1-, then again the system is BIM, because, the general term is computed in each step, consuming an increasing amount of CPU -this is, an unbounded CPU system- but still BIM.
If this is not the case, that is, if the system has an IS TF without an unknown algebraic expansion (UIS) -which is indeed real for most LTIs- then is not BIM in practice.
The condition here for an IS to be not rational is having no poles and no zeros, that is, being bounded:
$$0<H(z)=\sum_0^\infty c_nz^{-1}<\infty$$ $|z|>z_0$ for some $z_0$.
If the series crosses by zero through $z_i$, or diverges for some $p_i$, then it has a rational set of component $(1-z_iz^{-1})$ or $(1-p_iz^{-1})$, which is nothing more than a polynomial that can be factorized, being this condition applicable in a recurrent way.
Hence our definitions:
ALL $\leftrightarrow$ FIR $\lor$ IIR
IIR $\leftrightarrow$ TF $\lor$ IS
IS $\leftrightarrow$ KIS $\lor$ UIS
And our known relations:
FIR $\rightarrow$ BIM (By construction)
TF $\rightarrow$ BIM (By construction)
KIS $\rightarrow$ BIM (Counter-Example 1)
UIS $\rightarrow$ not BIM (Infinite Terms)
And finally, a Bounded in Memory System is a FIR, a TF or a Known Expansion Infinite Series.
BIM $\leftrightarrow$ FIR $\land$ TF $\land$ KIS
Counter Example 1:
The LTI system, with IIR defined by its Z-transform, is BIM, but has not rational TF:
$$H(z)=log(1+az^{-1}), \\ h(t)=\frac{(-1)^{t+1}a^t}{t}\Gamma(t), t>0$$
We are using the exploit that:
$$log(1+x)=\sum_{n=1}^\infty \frac{(-1)^{n+1}x^n}{n}, x<1$$
which in general is not valid, thus this procedure will produce unbounded LTI systems except for the trivial power series functions $sin(x)$, $cos(x)$,$e^x$...
We are also assuming that our TF must be finite in terms, which is the purpose of having a TF, provided that all Impulse Responses can be converted into series in Z space.
http://www.comm.utoronto.ca/~dkundur/course_info/362/6_Kundur_Rationalz-Transf_handouts.pdf
Counter Example 2:
The IIR AR System has both a rational and an infinite series transfer function:
$$H(z)=\frac{1}{1-az^{-1}}=\sum_{n=0}^{\infty}a^nz^{-n}, |a|<|z|$$


When should I calculate PSD instead of plain FFT magnitude spectrum? I have a thirty-second speech signal that was sampled at 44.1 kHz. Now, I would like to show what frequencies the speech has. However, I'm not sure what would be the best way to do that. It seems sometimes one calculates the absolute value of a Fourier transform, and sometimes power spectral density. If I understand correctly, the latter works so that I divide my signal into parts, do FFT part-by-part and somehow sum these. Window functions are somehow involved. Can you clarify this a bit for me? I'm new to DSP.

Now, I would like to show what frequencies the speech has. However,
  I'm not sure what would be the best way to do that. It seems sometimes
  one calculates the absolute value of a Fourier transform, and
  sometimes power spectral density.

If you want to attach physical meaning to your analysis, then go with the power spectral density, (PSD). This is because this will simply give you the power of your signal, in each frequency band. On the other hand if you do not want/care about a physical meaning, but want to know how the fourier amplitudes of each band vary relative to each other, you can stick to absolute magnitude.

In practice, you can compute the PSD as simply the absolute magnitude of the fourier transform squared. For example, if your signal is $x[n]$, and its DFT is $X(f)$, then the absolute magnitude of the DFT is $|X(f)|$, while the PSD is $|X(f)|^2$.

  If I understand correctly, the latter works so that I divide my signal
  into parts, do FFT part-by-part and somehow sum these. Window
  functions are somehow involved. Can you clarify this a bit for me? I'm
  new to DSP.

No, this is not true. What you are talking about here refers to the Short Time Fourier Transform, (STFT). This is simply chopping up your time domain signal, widowing it, and then taking the fourier trnasform. At the end of the day though, you will still have a complex matrix. If you choose to take its absolute magnitude, you will have an absolute magnitude fourier transform matrix. If you take its absolute magnitude squared, you will have a power spectral density matrix.


Frequency-domain zero padding - special treatment of X[N/2] Suppose we wish to interpolate a periodic signal with an even number of samples (e.g. N=8) by zero-padding in the frequency domain.

Let the DFT X=[A,B,C,D,E,F,G,H]
Now let's pad it to 16 samples to give Y.  Every textbook example and online tutorial I have seen inserts zeros at [Y4...Y11] giving
Y=[2A,2B,2C,2D,0,0,0,0,0,0,0,0,2E,2F,2G,2H].
(Then y = idft(Y) is the interpolated signal.)

Why not instead use
Y=[2A,2B,2C,2D,E,0,0,0,0,0,0,0,E,2F,2G,2H]?

As far as I can tell (my math knowledge is limited):

It minimises the total power
It ensures that if x is real-valued then so is y
y still intersects x at all the sample points, as required (I think this is true for any p where Y=[2A,2B,2C,2D,pE,0,0,0,0,0,0,0,(2-p)E,2F,2G,2H])

So why is it never done this way?

Edit: x is not necessarily real-valued or band-limited.

Let's look at the frequencies of the bins in your 8-point DFT:

$$
\begin{array}{c}
\omega_A = 0,\\
 \omega_B = \pi/4,\\
 \omega_C = \pi/2,\\
 \omega_D = 3\pi/4,\\
 \omega_E = \pi = -\pi\ ({\tt mod}\ 2\pi),\\
 \omega_F = 5\pi/4 = -3\pi/4\ ({\tt mod}\ 2\pi),\\
 \omega_G = 3\pi/2 = -\pi/2\ ({\tt mod}\ 2\pi), \\
 \omega_H = 7\pi/4 = -\pi/4\ ({\tt mod}\ 2\pi)
\end{array}
$$
So when you interpolate by a factor of 2, point $E$'s frequency becomes $-\pi$ or $+ \pi$.

At a first glance, I can't see what the problem is with your approach as it's not clear whether $E$ should be put into the bin associated with $\pi$ or $-\pi$.

On Julius O. Smith III's page, he states a condition:

  Furthermore, we require $x(N/2) = x(-N/2) = 0$  when $N$  is even, while odd  requires no
  such restriction.

And his example there is for an odd $N$, which avoids the problem.

Not sure it's required, but here's the full reference to Julius's work:

Smith, J.O. Mathematics of the Discrete Fourier Transform (DFT)
with Audio Applications, Second Edition,
http://ccrma.stanford.edu/~jos/mdft/, 2007, online book,
accessed September 28, 2011.


I/Q plane plot understanding? I am experimenting with some signals. I ma generating a 2.47 Ghz signal from signal generator and receiving it by an SDR (LMS7002m Rf chip) with Rx frequency tuned to 2.47 GHz. Upon plotting the received I/Q data I get the following plot. Can anyone help me understanding why there is a second circle below the circle centered at (0,0). Ideally I should be seeing a single circle..

IQ modulators, demodulators, and test equipment (oscilloscopes, etc.) can have channel imbalance in either or both gain and offset.  It looks like the Q channel of either the source or instrumentation changed or adjusted its offset during your measurement.


How to filter uniform noise in a microcontroller? By uniform noise I mean a band of a fixed width, no spikes, etc.  This has to be one of the simplest examples. I'm looking for the appropriate theoretical approach as much as something relatively computationally efficient.

The true signal changes slowly (slower than 1Hz, a battery voltage).  The noise band remains fairly constant (~50mV noise on a ~4V signal).

Initial ideas:

Long moving average (e.g. a 32-point moving average).
Leaky-integrator/low-pass filter/FIR filter:
"Low-pass filter" in non-EE, software API contexts
Recursive arithmetic average:

y_new = (new_sample + y_last) / 2

I'd love to learn more about what I think is a geometric average.  I'm not sure I can quantify the time-to-settle/time-to-seed for these options either.

If there are more standard approaches, I'd love to learn about those as well.

You need to check the noise bandwidth using an oscilloscope. By saying noise I mean the hardware noises before digitalization by ADC.
If the noises have high frequency component, remove it by analog filters. High frequency is defined as frequency higher than 1/2 of the sampling rate. And in practice, usually a sampling rate of 8-10 times of the pass band of the analog low pass filter is used. You may try simple RC filters.
Use FIR in a MCU to manipulate your data. Depends on the speed requirement, you can choose number of taps to use and etc. You should be able to find libraries provided by your chip manufacture, since FIR is very common in use.

Hope it helps!