Is there a method to combine speech signals recorded by different microphones? I have 3 speech files of the same speech recorded by 3 different microphones at the same time.

I'm only concerned to improve the speech quality (in terms of SNR and intelligibility) by combining the 3 speech files that I have.

Is there any method to solve the problem? What is the term used for it?

The techniques that I found involving multiple microphones was related to microphone array, but I'm interested to work with randomly placed microphones.

The easiest approach would just be averaging – since noise on these three microphones should at least partially be uncorrelated, that would improve SNR.

However, that requires that the three microphones are

synchronously sampled (i.e. the $n4th sample from each microphone happened "exactly" at the same time) and
in the same distance from the desired source of sound.

Since both 1. and 2. are usually not fulfilled, you'll have to look into adapting to that – the result is a beamformer. This might be a term often linked to regularly spaced microphone arrays, but that's not a mathematical requirement – you'd just start by detecting the phase/time offset between the strongest (or, better: the desired) source of sound on each of your microphone channels, and then align (read: fractional delay) the channels in a manner that maximizes SNR for that. This inherently leads to a geometric beam shape, though it usually is but a singular value decomposition of a receive signal matrix (that's where "beamforming" and "MIMO processing" meet).


What is the difference between phase noise and frequency noise? I'm reading Audoin and Guinot's The Measurement of Time: Time, Frequency and the Atomic Clock, and ran across a confusing an interesting passage (Section 5.2.5, pp 72-73):

  The physical origins of white phase and frequency noise are
  understood. They are associated with thermal noise in the radio region
  or shot noise from a particle flux.

There is also a chart that explicitly lists out phase noise and frequency noise as different components that lead to time and frequency instability.

I have always understood phase noise to be what is described in the Wikipedia article:

$$v(t) = A\cos(2\pi f_0t + \varphi(t))$$

where the $\varphi(t)$ portion is the phase noise of the sinusoidal signal.

I can't determine from the context of this book what frequency noise means, nor do any of my books mention frequency noise, just phase noise.

So, what is the difference between phase noise and frequency noise?

Phase Noise and Frequency Noise are not two different noise sources, they are artifacts of the same noise, it is just a matter of what units you want to use.  Frequency and Phase are directly related as frequency is phase changing with time, so if you have one you will always have the other; frequency and phase are related by derivatives and integrals: the derivative of phase equals frequency  (as a change in phase versus time is frequency).

For a waveform in the time domain that represents the frequency changing with time as $f(t)$, the equivalent function as phase changing with time would be:

$$\phi(t) = \int f(t) dt$$

And likewise to go from phase to frequency:

$$f(t)= \frac{d\phi(t)}{dt}$$

Note this can be very confusing as we are talking about a waveform whose value represents frequency changing with time; so to be very clear this waveform I am describing is indeed in the time domain (as the independent variable is time), and we are observing how the frequency changes versus time. We can take the Fourier Transform of this waveform, and then observe the rate of change of frequency (which is the frequency of our frequency variable....even more confusing!).

However getting that straight is important to understanding the difference between Phase Noise and Frequency Noise:

Integration in the time domain is a low pass filtering response in the frequency domain, as indicated by the Laplace Transform property for integration:

$$L\{ \int x(t)dt\}=\frac{1}{s}X(s)$$

And similarly the derivative in time results in a multiplication by s in Laplace, which represents a high pass function (as you can get the frequency response (or Fourier Transform) for both by setting $s=j2\pi f$)

Therefore a Phase power spectral density and Frequency power spectral density are related as follows:

$$S_f(f) = (2\pi f)^2 S_{\phi}(f) $$

Where $S_f(f)$ is the power spectral density due to frequency fluctuations, and $S_\phi(f)$ is the power spectral density due to phase fluctuations.

For example, if you had a white frequency noise process (meaning the FREQUENCY power spectral density was flat for all frequencies), the PHASE power spectral density would be going down at the rate of $1/f^2$ or 20 dB/decade.

Here is a plot I have showing the relationship between Phase Noise and Frequency Noise as time domain plots. Notice how the phase is smoother but wanders off relatively slowly with varying random offsets (low pass with random walk), while the frequency plot has higher frequency content but does not deviate from 0 for long durations (high pass filtered).


Gentle request for explanation on LPC and LPCC coefficients :) If I understand what LPC does, I think. It predicts the next sample by taking into account previous $p$ samples. I don't understand why it is useful in DPS and speech recognition.

I also came across LPCC that can be derived from LPC and are related to sound frequency. I also find it hard to follow the idea why cepstrum coefficients are useful or even better for speech recognition.

I did read about MFCC, though. I understood that before DCT, the coefficients (no MFCC yet) are simply energies at particular frequency range of a signal that is filtered by a triangle band pass filter. I do not have a problem with MFCC apart from the last step when DCT is applied. I think it also has something to do with cepstrum.

Why do you think cepstrum is so useful for speech, speaker and music recognition?

I'd be grateful for some words on that :)

EDIT. I have found a typo. Instead of "I have a problem with MFCC apart from at the last step when DCT is applied" I really meant "I do not have a problem with MFCC apart from the last step when DCT is applied".

Could someone please add a tag "LPCC" to this topic? It is one of the techniques in speech I think and deserves to be tagged. I cannot do it at the moment. Cheers!

As you mention, linear prediction attempts to estimate the next signal sample from a linear combination of the previous P outputs.  Mathematically, we can then express the $n^{th}$ sample of our signal $x[n]$ like this:
\begin{equation}
x[n] = \sum^{P}_{k=1} a_k x[n-k]+e[n]
\end{equation}
where $e[n]$ is the error signal.  Applying linear prediction to a frame of our signal of, say, length $N$ samples ($N\gg P$) gives us the linear prediction coefficients ${a_1,a_2,\cdots,a_P}$.

For signals where this model is valid, the error signal $e[n]$ contains very little information compared to $x[n]$.  The linear prediction operation then has given us $P$ numbers which represents most of the signal, as opposed to $N$ of them; in other words, we've been able to code the signal, hence linear predictive coding.  This coding of the signal is enormously useful for many different applications.  (Your mobile phone determines the linear prediction coefficients of your voice hundreds of times a second while you call your friends!)

Moreover, the linear predictive coefficients are essentially the coefficients of the optimal all-pole filter which describes $x[n]$. (In what sense of optimality depends on the algorithm used to determine them, but it's generally some form of time-domain least mean squares).  This means that we have an estimate of the frequency-domain envelope of $x[n]$, which we can derive from our coefficients.  Speech recognition generally operates from a frequency-domain perspective, and so linear prediction coefficients can be useful for that application too.

For improved recognition performance, it is useful to represent the frequency spectrum in other forms more closely linked with the way that we hear, i.e. some kind of perceptual domain. MFCCs and other similar signal representations (I assume that LPCCs do something similar, as I'm not familiar with them) do just that and warp the linear frequency domain in a way such that perceptual differences will be more apparent.


should I apply low-pass filter when calculating central derivative? Let us say that we have a discrete signal $I_n$, $n=0, 1, 2, ...$. According to Nyquist theorem the maximum frequency for such discretization is $f_{max} = 0.5$.

Now imagine that I want to calculate the derivative $D_n$ for this signal. The simplest approximation is the right hand side derivative

$$D_n = I_{n+1} - I_n$$

But what if I want to use a central derivative instead?

$$D_n = \frac{1}{2}(I_{n+1} - I_{n-1})$$

Just looking at this relation, it seems like I am increasing the sampling time from $1$ to $2$. Does this mean that before calculating central derivative I should apply low-pass filter to kill all the frequencies $f > 0.25$? Or am I being wrong and the simple fact that I am using the spacing = 2 for calculating derivative does not mean that I increase the sampling distance, because I can still calculate this derivative for each value of $n$?

No matter what the correct answer is (yes or no), can you please explain in more detail why? If the low-pass filter should be applied, can you explain if there is a method, how given discrete time relation I can infer to what degree I should smoothen the signal before using such relation?

No, you are not increasing the sampling time.  The sample rate stays the same.  The Central difference derivative filter is simply longer than the First difference derivative filter.  A longer filter just means more delay before you get your result.

No, you should not low pass filter without a clear reason for why it is needed to meet your signal processing objectives.  Your reasoning is flawed in this case, so don't apply a LPF.

Note that the Ideal derivative filter multiplies a signal's spectrum by $j\omega \cdot e^{j\omega \mu}$ in the frequency domain.  In the time domain, this ideal filter is an infinitely long filter with taps:

$$h[n] = \dfrac{1}{n+\mu}\left[\cos(\pi \left[n+\mu\right])-\mathrm{sinc}(n+\mu)\right]$$

Where

$n$ is the filter tap index in $(-\infty, \infty).$

$\mu$ corresponds to an inter-sample shift of the filter taps in the range $[0.0, 1.0]$ samples.

Given the above expression for $h[n]$, note that:

The Central difference is $h[n]$ truncated to 3 taps with $\mu = 0$.

The First difference is a scaled version of $h[n]$ truncated to 2 taps with $\mu = 0.5$.

Update to add filter magnitude responses plot


gnuradio: Best way to sync to simple preamble I never figured out a clean way to sync to a simple preamble in GNU Radio 3.7/3.8, so I built my own block to do so. It was a variation on the Correlate Access Code - Tag Stream block that simply cut out the code relating to the header length fields and supplied a property for the user to hardcode the payload length in the flowgraph. This is an extremely useful thing to do when reverse engineering RF systems. You can find the blocks here, along with a few other custom blocks:
https://github.com/paulgclark/gr-reveng
With GNU Radio 3.9 out, I was about to start the process of porting these blocks over to the new OOT scheme, but I wanted to make sure I wasn't duplicating features already in the tool somewhere. I haven't found them anywhere, but I know there's a bunch of great new stuff in there like the new Symbol Sync block, so I thought I'd ask. I have not found a way to get the Protocol Parser to implement this without customization, but I may be missing something.
Any other ideas, or should I just start porting my block?
Thanks
PC

After a bunch of searching, I discovered that Daniel Estévez has a block in his gr-satellites module that does exactly what I'm looking for, and even pulls in a few downstream blocks to clean up the flowgraph:
https://github.com/daniestevez/gr-satellites
The key blocks are:
Satellites->Packet Operators->Sync and create PDU
Satellites->Packet Operators->Sync and create packed PDU
The first is helpful if you have non-mod-8 length payloads, and generates PDUs with unpacked bytes (everything is 0x00 or 0x01) corresponding to your payload bits. The second packs each 8 payload bits into a packed byte (0x00-0xFF).
Previously, it was necessary to route things as follows:
Correlate Access Code -> Repack Bits -> Tagged Stream to PDU
Now, however, the gr-satellites blocks take care of all three operations.
I know the bulk of the gr-satellites module is dedicated to things much more advanced than this, but thanks to Daniel for also providing several nice, general-purpose blocks as well!


How to compute magnitude and phase response from transfer function in Z-domain? I have a transfer function $$H(z)=\frac{1+1.2z^{-1}+0.8z^{z^-2}}{1-0.9z^{-1}}$$ from which I'm supposed to sketch the magnitude and phase response. I know that you can transform $z=e^{j\omega}$ to get the frequency response, but I don't really understand what it means. The magnitude response is then $$\left|H(e^{j\omega})\right|=\left|\frac{1+1.2e^{-j\omega}+0.8e^{-j2\omega}}{1-0.9e^{-j\omega}}\right|,$$ but how do I calculate the values from here?

For example, if we set $\omega=0$, then the result is $30$. But if we set $\omega=\pi/4$, then the result is $$\left|\frac{1+1.2e^{-j\pi/4}+0.8e^{-j\pi/2}}{1-0.9e^{-j\pi/4}}\right|.$$ How does one solve this?

I also know that the formula for phase response is $$\arg(H(e^{j\omega}))=\frac{\Im(H(e^{j\omega}))}{\Re(H(e^{j\omega}))}$$ but I'm having a hard time to understand how to use it and what phase response actually is. Can you give me an example of how to calculate the phase response of a certain angle?

Use euler identity to expand the exponential to cos and sin. Then rationalize the denominator, and the last step is to group it to real part and imaginary part. To compute the magnitude just use the pythagoras theorem.

For an example I'll use simple transfer function like: $$H(z) = \frac {z}{z-0.5}$$
Then substitue the $z$:$$H(e^{jw})=\frac{e^{jw}}{e^{jw}-0.5}$$
Expand to cos and sin:
$$H(e^{jw})=\frac{cos(w)+jsin(w)}{cos(w)+jsin(w)-0.5}$$
Rationalize the denominator, group into real and imaginary part:
$$\begin{eqnarray*}H(e^{jw})&=&\frac{cos(w)+jsin(w)}{cos(w)+jsin(w)-0.5}\cdot\frac{cos(w)-jsin(w)-0.5}{cos(w)-jsin(w)-0.5}\\
&=&\frac {e^{jw}\cdot (e^{-jw}-0.5)}{(e^{jw}-0.5)\cdot (e^{-jw}-0.5)}\\
&=&\frac {e^{jw-jw}-0.5e^{jw}}{e^{jw-jw}-0.5e^{jw}-0.5e^{-jw}+0.25}\\
&=&\frac {1-0.5e^{jw}}{1.25-0.5cos(w)-j0.5sin(w)-0.5cos(w)+j0.5sin(w)}\\
&=&\frac {1-0.5cos(w)-j0.5sin(w)}{1.25-cos(w)}\\
&=&\frac {1-0.5cos(w)}{1.25-cos(w)} - j\frac {0.5sin(w)}{1.25-cos(w)}\\
\end{eqnarray*}$$
Compute the magnitude:
$$ \begin{eqnarray*} \left|H(e^{jw})\right| &=& \sqrt{\Re(H(e^{jw}))^2+\Im(H(e^{jw}))^2}\\
&=& \frac {\sqrt{(1-0.5cos(w))^2+(0.5sin(w))^2}} {1.25-cos(w)}\\
&=& \frac {\sqrt{1-cos(w)+0.25cos^2(w)+0.25sin^2(w)}} {1.25-cos(w)}\\
&=& \frac {\sqrt{1-cos(w)+0.25(cos^2(w)+sin^2(w))}} {1.25-cos(w)}\\
&=& \frac {\sqrt{1-cos(w)+0.25}} {1.25-cos(w)}\\
&=& \frac {\sqrt{1.25-cos(w)}} {1.25-cos(w)}\\
\end{eqnarray*}$$
And you're done!


multiplication of a function with a Fourier-transformed equals to Fourier-transformed with a function I already showed b item using the fact that it is $h\left(0\right)=\int \:f\left(t\right)g\left(0-t\right)dt$
I struggle a lot of hours trying to find the trick in item C.
Can anyone help please ?

Hint:
You only need to use the definition of the Fourier transform and its inverse transform to show that if
$$G(f)=\mathscr{F}\{g(t)\}$$
then
$$g(f)=\mathscr{F}\{G(-t)\}$$
holds. Then just use the result from $b$.


Create images with sub-pixel shift I want to test some motion estimation algorithms which compute motion with sub-pixel accuracy. I can create images shifted by the order of pixel (e.g. [2,2]) using MATLAB or OpenCV but I am not clear on how can I shift images with sub-pixel (e.g. [3.57,4.83]) accuracy. Any help regarding this will be much appreciated.

The science of finding values in between those points where you actually have samples is called Interpolation.

It's a rich field, and every signal processing toolkit (be it Matlab or OpenCV, or …) has functions for it.

Picking the right interpolator has a lot to do with the nature of your input, and the purpose of you doing the interpolation  – so, no general advice can be given on what specific interpolator to use.


How to create a sine wave that transitions through frequencies and returns to phase For an auditory psychological experiment, I need to create a sine wave, such that at a given point the frequency of the wave will rise smoothly, and descend back to base frequency. I need it to return in phase with a second sine wave which is a continuous sine wave of the base frequency:

ie: one ear will get a continuous 300 Hz sine wave,
    the other will get a 300 Hz sine wave, that rises to 305 for a short while, and descends back to 300 Hz in phase with the second ear.

I generated the frequency modulation with a complex phasor recursion, as suggested here: https://dsp.stackexchange.com/a/1087/18658

ie, each point in the sine wave is given by

z(n+1) = z(n) * exp(j* 2pi* frequency).

The parameter I can play with is the length of modulation. Any ideas about how to find the nearest length that gives me the in-phase return?

For the first continuous $f_0=300\mbox{Hz}$ signal, after $2N$ samples (for a duration of $2N/f_s$ at a sampling rate $f_s$) the phase will have changed by:
$$
\begin{align}
 \Delta \phi_1
   &= 2\pi \frac{f_0}{f_s} (2N-1)
\end{align}
$$

On the other hand, as the second signal ramps up linearly (in small quantized steps as per your phasor implementation) from $f_0=300\mbox{Hz}$ to $f_1=305\mbox{Hz}$ then down from 305Hz to 300Hz over the same $2N$ samples, the phase will have changed by:
$$
\begin{align}
 \Delta \phi_2
   &= \sum_{n=1}^{N} 2\pi \left(\frac{f_0}{f_s} + \left(\frac{f_1}{f_s}-\frac{f_0}{f_s}\right)\frac{n}{N}\right)
      + \sum_{n=1}^{N-1} 2\pi \left(\frac{f_1}{f_s} + \left(\frac{f_0}{f_s}-\frac{f_1}{f_s}\right)\frac{n}{N}\right) \\
%   &= 2\pi \frac{f_1}{f_s}
%      + 2\pi \sum_{n=1}^{N-1} \left(\frac{f_0}{f_s} + \left(\frac{f_1}{f_s}-\frac{f_0}{f_s}\right)\frac{n}{N} + \frac{f_1}{f_s} + \left(\frac{f_0}{f_s}-\frac{f_1}{f_s}\right)\frac{n}{N}\right) \\
%   &= 2\pi \frac{f_1}{f_s}
%      + 2\pi \sum_{n=1}^{N-1} \left(\frac{f_0}{f_s} + \frac{f_1}{f_s}\right) \\
   &= 2\pi \frac{f_1}{f_s} + 2\pi \left(\frac{f_0}{f_s}+\frac{f_1}{f_s}\right) \left(N - 1\right)
\end{align}
$$

For those two signals to be in phase after $2N$ samples, their phase difference must be a multiple of $2\pi$. Thus the constraint is that
$$
\begin{align}
  \Delta \phi_2
    &= \Delta \phi_1 + 2\pi k \\
  2\pi \frac{f_1}{f_s} + 2\pi \left(\frac{f_0}{f_s}+\frac{f_1}{f_s}\right) \left(N - 1\right)
    &= 2\pi \frac{f_0}{f_s} \left(2N-1\right) + 2\pi k
\end{align}
$$
for some arbitrary integer $k$, or equivalently:
$$
\begin{align}
  N &= k \frac{f_s}{f_1-f_0}
\end{align}
$$
Which, provided that $kf_s/(f_1-f_0)$ is an integer greater than 0, yields a signal duration of
\begin{align}
  \frac{2k}{f_1-f_0}
\end{align}

In your specific case for a frequency excursion of $5\mbox{Hz}$, the signals will return in phase for modulation durations of 0.4, 0.8, 1.2, ... seconds.

Note that if you need to sustain the signal for some $M$ samples at $f_1$ (instead of immediately going back down after the rise), then the constraint becomes:
$$
\begin{align}
  N+M &= k \frac{f_s}{f_1-f_0}
\end{align}
$$
In other words the duration of the rise and plateau together would, in your case, need to be a multiple of 0.2 seconds. So for example if you want to have a plateau of 1 second at 305Hz, you might choose a rise + plateau of 1.2 seconds, giving 0.2 seconds for the ramp up and another 0.2 seconds for the ramp down (for a total signal duration of 1.4 seconds)


chord recogination by multiple fundamental frequency estimation I'm gonna implement multiple frequency estimation by summing harmonics amplitude!
I've read Anssi Klapuri paper but i can't understand him.
first of all he said spectral whitening i did it by i don't know whether it's true! after that i have to implement this formula:
$$s(\tau)=\sum_{m=0}^Mg(\tau,m)max_{}|Y(k)|$$
i can't understand what this $\tau$ is!
please give some hint to implement it or some other method for implementation.

$g(\tau,m)$ is the weighting function used in the paper. Equation (7) tells you how to calculate its value.

$\tau$ is the fundamental period which is the reciprocal of the fundamental frequency, $f_0$. If the fundamental frequency is 440 Hz then $\tau$ is $1/440 = 2.27$ ms.


Is there any way that we can perform speech recognition without using Fourier transforms? I am trying to research about speech recognition and why everyone uses Fourier transforms in going about the topic. I know that we get information related to the frequency of each sound uttered which is extremely important in the process, however I wanted to know if the amplitude of a sound wave can be used in order to recognize phonemes and the information related to a spoken word.

All recognition tasks (doesn't even have to be speech recognition) are reductions of a very high-dimensional signal (your speech recording's dimension is the number of audio samples!) to a low-dimensional signal.

As such, it is generally advisable to transform the input signal through an easy operation to a representation where the dimensionality can be more easily reduced. For speech, the frequency domain is such a representation – physically, only a very limited set of discrete frequencies is sufficient to make up speech.

Hence, having a frequency-domain transform (typically, in speech recognition) early on in the process makes a lot of sense: after the transform, only very few coefficients contain information, and need to be evaluated further.

Hence, since speech models are based on the physics of speaking, and that is a time-frequency physical problem, it's very likely that an algorithm that actively tries to stay away from frequency domain will not perform good.

Nowadays, however, many classification problems are solved using Deep Neural Networks; these basically are just a network of nonlinear functions with a constant factor and offset for each of them. Training such a Neural Network to do something useful is just finding the right factors and offsets. There's heuristics that allow us to do that automatedly at large scale.

Such networks can of course be trained directly on discrete time-domain speech signals, and might have good recognition rates after sufficient training has happened.

Notice that training such Neural Networks is not a trivial task in itself; it's usually an intentionally or inherently stochastic procedure, and might or might not lead to a functioning classifier. There's a lot of experience and luck involved building a good Neural Network classifier instead of just "something that does something close to what I want".

However, the size and hence cost, and the effort spent on training, and the chance that learning doesn't lead to a working result, can be influenced by the way the signal is preprocessed. Think of it like this: Imagine you'd have to learn a language by listening to someone through a bad phone connection, compared to having a good phone at each end of the line that does some error correction and sound enhancements.

A typical preprocessing step here would be some sort of frequency-domain transform. Simply because the result of that more compactly describes the speech content than the time representation, it'll also make the Neural Network training's job of finding important and unimportant parts of the signal easier.

So, to conclude:

Yes, you can avoid the extremely simple Fourier transform (or one of the many other frequency-domain representations, many of which are based on the Fourier transform), it would probably make your life much harder than necessary.


Uniform noise signal Good afternoon. I have some heart monitor data that I am attempting to analyze. I have been given time between the peaks (in milliseconds), and have attempted to reverse-engineer a signal from this using the following code.

clc;
clear all;
%Sampling frequency
Fs = 250;
%total data time (s)
ts = xlsread('12_RR_1.xlsx','G:G');
%total data time (ms)
%total_time = xlsread('12_RR.csv.xlsx','I:I');
%Point in time where there is a max point
xn = xlsread('12_RR_1.xlsx','H:H');
xnt = size(xn);

t = 1:ts;
time = transpose(t);

for i = 1:ts
    for j = 1:xnt(1,1)
        z(i,j) = xn(j,1);
    end
        i = i+1;
        j = j+1;
end

for i_2= 1:97484
    for j = 1:xnt(1,1)
   if time(i_2,1) == z(i_2,j)
       z_raw(i_2,j) = 1;
   else
       z_raw(i_2,j) = 0;
   end
    end
    i_2 = i_2+1;
   j = j+1;
end

 peaks = sum(z_raw,2);
 Peaks = abs(peaks);
  N = length(peaks);

 subplot(2,1,1)
 plot(t,Peaks)
 title('Raw Data')
 subplot(2,1,2)
 plot(t/Fs,fft(Peaks,ts))
 title('FFT')

This produces a binary signal where y=1 if a peak is detected and y=0 if it is not. However, when I attempted to take the FFT of this it returns a signal with weirdly uniform noise, as shown below.

What's going on--why does a seemingly low-frequency signal have a uniform frequency response? What can I do to my original data (the time between heart impulses) to extract more meaningful information and get a more realistic FFT (high amp low freq data with small amp high freq data)

Your signal is a superposition of fairly widely spaced impulses. Each impulse by itself has a flat magnitude spectrum and the intervals between impulses are only pseudo periodic so the phase response superpositions are effectively incoherent. The problem isn't your FFT, it is your signal model.

An actual ECG is much smoother. If you expect to see a predominant frequency at around 1Hz, the signal would look like a sine wave of around 1 Hz, not a sequence of isolated spikes.

You probably really want to measure the time intervals between spikes directly in the the time domain.


What is the filter to get a heavy metal film ilustration? How I can transform an image to get a same image style like the Heavy Metal film?.

I remember that I could read a paper about this, but I lost it.

There are two techniques that can achieve similar results, but if you caught that in a recent paper then you are probably after Style Transfer.

It is based on "teaching" a neural network how to combine specific elementary features in one context (e.g. Renaissance paintings) and then applying them in a different context (e.g. some random picture).

The other technique is "Colour Transfer" (or "image re-colouring") and that one attempts to copy just the colour tones and overall "feel" from one image to another but without changing the structure of its content. Here is a good overview and a particular method of applying it.

Hope this helps.


Recovering time domain signal from sum of independent and identically distributed signals Background: I observe a sample of a variable z that is the sum of two independent and identically distributed variables x and y. I'm trying to recover the distribution of x, y (call it f) from the distribution of z (call it g), under the assumption that f is symmetric about zero. According to Horowitz and Markatou (1996)  we have that the Fourier Transform of f is equal to sqrt(|G|), where G is the Fourier transform of g.
Example:
import matplotlib.pyplot as plt
import numpy as np
from scipy.stats import gaussian_kde, laplace

# sample size
size = 10000
# Number of points to preform FFT on
N = 501
# Scale of the laplace rvs
scale = 3.0

# Test deconvolution
laplace_f = laplace(scale=scale)
x = laplace_f.rvs(size=size)
y = laplace_f.rvs(size=size)
z = x + y
t = np.linspace(-4 * scale, 4 * scale, size)
laplace_pdf = laplace_f.pdf(t)
t2 = np.linspace(-4 * scale, 4 * scale, N)

# Get density from z. Kind of cheating using gaussian
z_density = gaussian_kde(z)(t2)
z_density = (z_density + z_density[::-1]) / 2
z_density_half = z_density[:((N - 1) // 2) + 1]
ft_z_density = np.fft.hfft(z_density_half)
inv_fz_density = np.fft.ihfft(np.sqrt(np.abs(ft_z_density)))
inv_fz_density = np.r_[inv_fz_density, inv_fz_density[::-1][:-1]]
f_deconv_shifted = np.real(np.fft.fftshift(inv_fz_density))
f_deconv = np.real(inv_fz_density)

# Normalize to be a pdf
f_deconv_shifted /= f_deconv_shifted.mean()
f_deconv /= f_deconv.mean()

# Plot
plt.subplot(221)
plt.plot(t, laplace_pdf)
plt.title('laplace pdf')

plt.subplot(222)
plt.plot(t2, z_density)
plt.title("z density")

plt.subplot(223)
plt.plot(t2, f_deconv_shifted)
plt.title('Deconvolved with shift')

plt.subplot(224)
plt.plot(t2, f_deconv)
plt.title('Deconvolved without shift')

plt.tight_layout()
plt.show()

Which results in

Issue: there's clearly something wrong here. I don't think I should need the shift, yet the shifted pdf seems to be closer to the truth. I suspect it has something to do with the domain of the IFFT changing with the sqrt(abs()) operation, but I'm really not sure.

The answer was given in a comment in a cross-post. I wasn't sure if it was a computational or conceptual issue I was having. Turns out it just an issue with the way FFT expects the input to be ordered and how the output is sorted.


Multichannel input/output evaluation board I'm working on a project at school and the first step is to choose a dsp/evaluation board but I'm having trouble. The project is to create a system that can be used for a variety of audio applications (such as sound localization or digital effects processing). Because of this the evaluation board should have at least 3 inputs and 1 output. Some other requirements at 16 bits/sample and at least a 16KHz sampling rate and the ability to operate in real time. The final result should be a setup where to do a new audio project, one would only have to do simple coding (preferably C, or maybe python). Thanks for any recommendations or general help, feel free to ask more questions about the project.

Ryan

This is a little bit of a tricky question because it is homework and a direct answer would simply equate to doing that homework which is not what this website is about. Apart, from the purely technical bit of the right decisions, which you should do entirely on your own, I am interpreting this question as "There are so many options out there! Where do I start? I am lost!" and I am providing the following answer on a general framework that can be used in these cases.

Very briefly: Use a requirements matrix to rank your options.

The requirements matrix lists each of the requirements of the project. These are already available (Minimum sampling frequency, C compiler, number and type of I/O, real-time requirements, etc). Some of the requirements of course will need further unpicking. For example, the requirement "Capable of doing sound localisation" immediately leads to the question of "How accurate does this have to be?" which tightens the requirements of sampling frequency and data processing. Is floating point processing required or fixed point? Is the board fast enough to emulate floating point? And other similar considerations.

For each design option, add one column and provide a simple answer of how well does it cover that requirement. For example, suppose that you have an audio design project and your possible implementation options are Product 1, product 2, product 3, ... and each one has different ways of providing an audio output (Simple filtered Pulse Width Modulation, audio chipset, virtual sound card, etc).

At the very minimum, this matrix will help in selecting the one option that covers the most requirements for the design. More elaborate scoring scales are possible but also application dependant.

For more information, please see this link.

Hope this helps.


Allan Variance vs Autocorrelation - Advantages I am currently studying oscillator stability and have come across the Allan variance. I gather that it was developed as an alternative to the standard variance as it doesn't necessarily converge for power law noise. However, I have a hard time seeing how this niche tool differentiates itself (or how it is superior) to other more conventional statistical tools like the autocorrelation or PSD.
Does anyone have any insight?

I currently work with the inner design details of Atomic Clocks so use the Allan Variance and Allan Deviation (ADEV) extensively. The primary point is that it can be used for non-stationary processes (which frequency noise is). For non-stationary signals where the autocorrelation or power spectral density can't be used with consistency, the Allan Deviation will converge to a consistent metric for comparison of the frequency accuracy of two different clocks (and generally for the comparison of the noise in non-stationary signals).
Non-stationary signals will often have a shorter time duration over which the signal can be sufficiently assumed to be stationary. Oscillator phase noise is an excellent example of this: for larger frequency offsets (typically 1 Hz and above for lower frequency oscillators such as 10MHz) a power spectral density (PSD) can be computed and given as the oscillator’s phase noise performance. The issue is this PSD cannot be computed consistently at lower frequency offsets as the non—stationary contributions start to dominate and this is where the ADEV computation will shine: in computing the longer term frequency accuracy of the clock over time durations where non-stationary contributions dominate (1/f noise, drift). ADEV is computed as an rms error given by the difference between the average frequency accuracy over interval tau seconds compared to the average frequency accuracy tau seconds ago. The rms of this error is computed over a very long duration (this error difference is a stationary signal and thus the standard deviation can be consistently computed).
One specific utility for me that the ADEV has provided well beyond the world of clocks is to determine the time constraints under which a random process can be sufficiently assumed to be stationary: if the random process is stationary we can continue to average to get a better estimate of the underlying mean— therefore we can use the ADEV to determine the optimum averaging time after which no further improvement can be gained (and through further averaging our estimate can get worst). This applies to everything from channel estimation to financial markets, wherever we may be observing past results to make an estimate of an underlying process.
I explain further details of the Allan Deviation and its potential use and application well beyond oscillator stability in these other posts:
How to interpret Allan Deviation plot for gyroscope?
What determines the accuracy of the phase result in a DFT bin?
Averaging data from 2 sensors
Is it possible to add some kind of noise to a non stationary (having unit root) signals to make them stationary?


Error introduced by smoothing filter at the end of a processing block I'm running analysis on a multichannel audio signal and due to the size have decided to process it in blocks (my computer doesn't have enough memory to process it in one go). Some of the data produced by the analysis needs to be smoothed.
The smoothing function means that each sample was an average of the samples contained within a window and the window being centered around sample n. At the edges this obviously isn't possible so the at the start of the signal the window remains static until such point that sample n is in the centre of the window, and then the window can slide along the signal - this also means at the end of the signal there will come a point where the window can't move along and will remain static. So essentially there will be a number of identical values at the start, and then at the end of the signal. This happens at the end of each block (as it can't continuously slide in the same way it would if processing the whole signal in one chunk). I know in order avoid discontinuties at the start of each block an overlap add/save method can be used, but I'm not sure how to fix this problem that occurs at the end of a block.
Just to explain further how the smoothing window moves along the signal, if we have a window length of 5 samples, the first 3 averaged values would be identical, at the 3rd sample, our current time step would now be located at the centre of our window, when moving to sample 4 the window would also slide along by one, keeping our current time step at the centre of the window, and now producing a different value. Once we get to the end of our signal the window will remain static for the last 3 samples, again producing identical values. This is fine when processing something in one big block, but produces discontinuties when trying to process in multiple blocks.

A moving average filter is a FIR filter of equal weights 1/N and a FIR filter is a degeneralization of a IIR filter.
In a non realtime implementation, it is functionally equivalent if you consider the filter to be N samples delay, N samples lookahead or N1 samples delay and N2 = N-N1 samples look ahead. Pick whatever time philosophy makes sense in your head and minimize tedious indexing.
Matlabs filter() and filtic() functions provide the means for carrying the  result from block to block for IIR filtering. Sadly only single channel support.
If you use something like:
conv2(ones(5,1), 1, X, ‘full’);
You will get an output vector that is larger than the input vector. Those surplus samples are exactly what you need to add to the initial samples of the next block.


Can Temperature Data be Predicted Using Adaptive Filter (Such As LMS) Algorithm? I am working on a project which requires me to implement adaptive filter as a predictor.
I have just started on adaptive filter and I intend to use least mean square algorithm for weight adjustment.

How can I predict future values from this system ?

Any help would be beneficial for me.
Thanks.

Yes you can predict future temperatures, based on past temperatures, using adaptive filtering as well.

The optimal linear estimation of a WSS random process from its past values, which is known as linear prediction, is given by a Wiener filter structure where the desired response to be estimated is the current sample of the input (current temparature in your case) and the filter input is the $N$ past samples of the input, (assuming one step forward prediction of order $N$).

The LMS adaptive filtering algorithm simply approaches this optimal Wiener predictor coefficients for WSS signals and for non WSS signals tries to continue to be optimal by tracking it.

This prediction mechanism does not depend on the physical origin of the signals but on their statistical characterisation. As long as your temperature data posses reasonable degree of correlation within it, then the filter will do its best to predict it.


Confusion about subtle difference between discrete-time and continuous-time In Alan Oppenheim's book Signals and Systems a comparison is made between the properties of discrete-time and continuous-time complex exponential signals in section 1.3 pg. 26.  Specifically it says:

  The continuous-time complex exponential $e^{j\omega_0t}$ has two
  properties:
      1.  the larger the magnitude of $\omega_0$, the higher is the rate of oscillation in the signal
      2. $e^{j\omega_0t}$ is periodic for any value of $\omega_0$

  Consider the discrete-time signal $e^{j(\omega_o+2\pi)n}= e^{j2\pi n}
e^{j\omega_0 n}=e^{j\omega_0 n}$. From this we see that the
  exponential at the frequency $\omega_0+2\pi$ is the same as that at
  frequency $\omega_0$.Thus we have a very different situation from the
  continuous-time case, in which signals $e^{j\omega_ot}$ are all distinct for
  distinct values of $\omega_0$

The part I don't understand is the one in bold . It seems to me I can do exactly the same thing for continuous-time complex exponential by replacing $n$ with $t$:

$e^{j(\omega_o+2\pi)t}= e^{j2\pi t}e^{j\omega_0 t}=e^{j\omega_0 t}$

So there we have it.  I have produced the same signal in the continuous-time complex exponential.  So what are they talking about? I see no difference.

The problem is that $$e^{j 2 \pi t} \neq 1,$$ except in the case when $t$ is an integer. However, in that case, you have the discrete-time formula: $n$ is assumed to be an integer. In most DSP literature, wherever you see $n$ used as in this case, you can safely assume it's an integer.


Scikit-Image, Numpy, and Selecting Colors (python) I'm trying to mask colored features from a photograph so that I can do some other processing on them. I've played with a few packages (scikit-image, mahotas, and openCV) and have settled on Scikit-image because it plays well with scikit-learn.

I'd like to accomplish two things, eventually: (1) get colored features so that I can compute things like length, and (2) remove colored features from the image while retaining the white spidering-looking veins (the lattice-looking stuff around the margin of the image will need to go, but that can be a later endeavour).

I'm new to image processing in Python, and I think I've mislead myself by using strange color spaces. Perhaps it would have been better to define an upper and lower boundary for each colour, and mask the image array? I'm definitely open to any flexible solution to the problem. Eventually I'll be applying this to hundreds of photos..

Original Photo

Code

from skimage import io, img_as_float, color, exposure

img = img_as_float(io.imread('./images/testimage2.JPG'))

# Isolate paint marks

# Put image into LAB colour space
image_lab = color.rgb2lab(img)
img = exposure.rescale_intensity(img)

# Colours of interest
color_array = np.array([
        [[[255, 255, 0.]]],    # Yellow stuff
        [[[255, 190, 200.]]],    # Pink stuff
        [[[255, 165, 0.]]],    # Orange stuff
        [[[255, 0, 0.]]],    # Red stuff
        ])

# Loop through the color array and pick out the colored features
for i in range(0,color_array.ndim):
    # Compute distance between the color of interest and the actual image
    # http://scikit-image.org/docs/dev/api/skimage.color.html#skimage.color.deltaE_cmc
    # "The usual values are kL=2, kC=1 for “acceptability”"
    distance_color = color.deltaE_ciede2000(color_array[i], image_lab, kL=2, kC=1, kH=0.5)

    # Normalise distance
    distance_color = exposure.rescale_intensity(distance_color)

    # Mask image
    image_filtered = img.copy()
    image_filtered[distance_color > 0.5] = 0

    # Plot it up
    print ("Filtered to: ", color_array[i])
    f, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(20, 10))
    ax0.imshow(img)
    ax1.imshow(distance_color, cmap='gray')
    ax2.imshow(image_filtered)
    plt.show()

Results

You can look at your picture in a different color space, e.g. the HSV space. In this space, each pixel has 3 components: hue, saturation and value. Hue defines the color's tone, saturation the color's strength and V the brightness of the pixel. For your particular problem, you can search for pixels that have e.g. a decent saturation (i.e. are reasonably colored and not just gray). The following code snippet uses OpenCV, but I believe you will be able to translate into scikit-image:

I = cv2.imread("/tmp/XApdz.jpg")  # Load the image

hsv = cv2.cvtColor(I, cv2.COLOR_BGR2HSV);  # Convert to HSV space
h, s, v = cv2.split(hsv);  # extract the single components

plt.figure(1, figsize=(10, 10))
plt.subplot(3,2,1); plt.imshow(h, cmap='gray'); plt.title('Hue')
plt.subplot(3,2,2); plt.imshow(s, cmap='gray'); plt.title('Saturation')
plt.subplot(3,2,3); plt.imshow(v, cmap='gray'); plt.title('Value')

mask = (s > 80).astype(np.uint8);  # Thresholding in the Saturation-channel
plt.subplot(3,2,4); plt.imshow(mask); plt.title('mask')

disk = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(3,3))  # remove very small regions
opened = cv2.morphologyEx(mask, cv2.MORPH_OPEN, disk);
plt.subplot(3,2,5); plt.imshow(opened); plt.title('Opened mask')

I2 = I.copy()
I2[opened.astype(bool), :] = 0; # Set the pixels to zero, where
plt.subplot(3,2,6); plt.imshow(I2)


How does one calculate a steering matrix for DOA given real data and only antenna positions I am currently trying to implement the MUSIC algorithm for a real dataset which captures a sound signal in a room using a cubic 8 sensor array setup.
Every example I've found for the MUSIC algorithm talks about the steering matrix and calculates it using a set of DOAs which are, well, the actual DOAs that the algorithm is trying to estimate? This part deeply confuses me.
All I have as inputs are:
X = [M x N] matrix, where M is the number of antennas, N is the number of timesteps of the signals. So essentially an 8 channel signal snapshot.
Y = [M x 3] matrix, which holds the 3D coordinates of each sensor.
It is worth noting that the sensors in my dataset are not actually a perfect cuboid. It's essentially two rectangles where the upper rectangle is above the other, but also slightly shifted in the X axis.
It is my understanding that a steering matrix should be computed using only Y. X is used later, for the eigenvector computation. So how would I go about generating this matrix if all I have are mic positions?
Any advice, literature would help me greatly.

Figured it out reading https://arxiv.org/pdf/2008.12323.pdf
Basically, you construct the matrix for all directions. so 360x180 for all elevations and azumoths.


What are the advantages and disadvantages of the Ridgelet transform? I used the continuous wavelet transform (CWT) to analyze the behavior of a 1-D smooth signal (extracting its modal property, damping and frequency components) in my research. Now someone asked:

  Why the presented CWT algorithm is better than Ridgelet transform?

I thought it is better to ask in this forum about disadvantages of Ridgelets.
For example:

it is said that Wavelet can handle point discontinuities whereas
Ridgelet is better to deal with line discontinuities. So, do we have
any line discontinuities in a 1-D signal?
what is computational complexity of Ridgelet in term of big $O$ i.e.
(like $O(N)$), for both continuous (in scale) and discrete versions?

Because of their ability to better represent signals, Ridgelet transforms are preferred over Wavelet transforms in general. However, to be used in practice it has to be discretized. Now, the discretization of Ridglet transform is a challenging task on its own as it involves interpolation in polar coordinates which makes perfect reconstructions (inversions) difficult. For example, Minh N. Do and Martin Vetterli propose an orthonormal version of Ridgelet Transform to achieve better invertability [1].

Ridgelet transform enhances the idea of point-to-point mapping of singularities to point-to-line mappings, which is more effective in handling directions. Yet, even Ridgelets might be ineffective for certain applications such as image processing, where edges are curved rather than being straight lines. In such cases, better transforms such as Curvelets are already proposed and proven to be more powerful.

In this regard, the method of choice really depends on your application. As wiki quotes:

  Continuous Wavelet Transform (CWT) is very efficient in determining the damping ratio of oscillating signals (e.g. identification of damping in dynamical systems). CWT is also very resistant to the noise in the signal.

If your problem is suitable to CWT, then you will be perfectly fine with using it. The Academia involves a lot of works in comparing and demonstrating the strengths of all of these transforms. Some of these would be [2,3,4].

Under the shade of the comments below, I now elaborate more on the theory.

First of all, singularity is simply the concept of point where the function or its derivative is not defined, the point of cusp, double point or any other simlar undefinedness, which I won't go into details here. You can check the wiki "Mathematical Singularity". A good measure of singularity strength would be Lipschitz exponent, which can well be characterized by Wavelet transforms. For more information on that take a look at [5].

But, as you ask, why are you interested in the singularities anyway? Well, because such irregularities basically play an important role in carrying the information that characterizes the signal [6].

Having explained the singularities a bit, I certainly agree that singular lines are a concern of multidimensional structures, more specifically $R^2$ ($L^2$ etc.). This is why transforms like Ridgelets and Curvelets are good at image denoising than Wavelets. Yet, this doesn't violate the fact that 2D Wavelets would also be useful in image analysis.

Regarding the computational complexity, even before analyzing the equations, you could rely on your knowledge that computation of the Ridgelet transform involves a Wavelet transform in the Radon domain. That is why it should definitely be heavier in terms of computational complexity, compared to a Wavelet transform.

.  REFERENCES

http://www.ifp.illinois.edu/~minhdo/publications/FRIT_trans.pdf
http://jstarck.free.fr/curvencyclop09.pdf
http://authors.library.caltech.edu/1381/1/STAieeetip02.pdf
http://www.iro.umontreal.ca/~mignotte/IFT6150/ComplementCours/CurveletTransform.pdf
http://www.iis.sinica.edu.tw/papers/whwang/5114-F.pdf
http://www.cmap.polytechnique.fr/~mallat/papiers/MallatSingDet92.pdf


When is cubic spline interpolation better than an interpolating polynomial? The following plot is a slight variation of an example in a text book. The author used this example to illustrate that an interpolating polynomial over equally spaced samples has large oscillations near the ends of the interpolating interval. Of course cubic spline interpolation gives a good approximation over the whole interval. For years, I thought high order polynomial interpolation over equally spaced samples should be avoided for the reason illustrated here.

However, I recently found many examples of bandlimited signals where a high order interpolating polynomial gives less approximation error than cubic-spline interpolation. Typically an Interpolating polynomial is more accurate over the entire interpolating interval when the sample rate is sufficiently high. This seems to hold when the samples are equally spaced with a sample rate at least 3 times greater than the Nyquist frequency of the signal. Furthermore, the advantage over cubic spline interpolation improves as (sample rate)/(Nyquist frequency) increasees.

As an example, I compare cubic-spline interpolation with an interpolating polynomial for a sine wave with a Nyquist frequency of 2 Hz, and a sample rate of 6.5 Hz. Between the sample points, tthe interpolating polynomial looks exactly the same as the actual signal.

Below I compare the error in the two approximations. As with the first example, the polynomial interpolation does worst near the beginning and end of the sample interval.  However, the interpolating polynomial has less error than a cubic spline over the whole sample interval. The interpolating polynomial also has less error when extrapolating over a small interval. Did I discover a well known fact? If so, where can I read about it?

The phenomenon being discussed is Runge's phenomenon.

The maximum absolute value of the $n$th derivative of $\sin(\omega t)$ is $\omega^n$. For Runge's function $\frac{1}{25t^2+1}$ the maximum absolute value of the $n$th (even) derivative is $5^nn!,$ where $n!$ denotes factorial. This is much faster growth. Only if the derivatives grow too fast by increasing $n$, then it is possible that interpolation error diverges as the interpolation order is increased. Exponential in $n$ is not yet too fast. Have a look at: James F. Epperson, On the Runge example, The American Mathematical Monthly, vol. 94, 1987, pp. 329-341.

If a function has only continuous derivatives, then the competing approach, piece-wise polynomial spline interpolation always converges if a small fixed number of its early derivatives are bounded over the interval of interest, see Wikipedia article on linear interpolation as an example.

If both methods converge, then (non-piecewise) polynomial interpolation has the advantage of a higher polynomial degree if many samples are used, and may provide a better approximation, as you saw in your sine example. You may also be interested in L.N Trefethen, Two results on polynomial interpolation in equally spaced points, Journal of Approximation Theory
Volume 65, Issue 3, June 1991, Pages 247-260. Quote:

  [...] in band-limited interpolation of complex exponential functions
  $e^{i\alpha x}\, (\alpha \in \mathbb{R}),$ the error decreases to $0$ as $n \to \infty$ if and only if $\alpha$ is
  small enough to provide at least six points per wavelength.

You have 6.5 samples per wavelength.


MATLAB's designfilt vs butter function I am relatively new to signal processing, and have always used MATLAB's designfilt option for my signal processing needs.  Recently, however, I have seen people just using MATLAB's butter function.  We all primarily just analyze EEG signals.  I am just curious as to what makes them do this, as I currently don't see a major difference between the two aforementioned functions.
Again, I am pretty new to signal processing.

The function butter, as the name indicates, is used to construct the Butterworth IIR filter.
The function designfilt can be used to design the Butterworth filter as well, but you can design the other type filters with this function.
What type filter you need depends on your application. For example, the Butterworth filter rolls off more slowly around the cutoff frequency than the Chebyshev filter or the Elliptic filter, and it may decide your selection.
From the point of view of a coder, the use of designfilt can be preferable even when constructing the Butterworld filter. The function uses a name-value syntax and enables one to create a digitalFilter object:
d = designfilt(resp,Name,Value)

, d is the digitalFilter object. The digitalFilter object bundles all the parameters of a linear filter into a single container. The object digitalFilter is immutable, and after it is created, you can edit it only with the Filter Design Assistant that you open with the command designfilt(d).
On the other hand, the syntax of butter is simpler, the filter type stands out, and the function explicitly returns the matrix of coefficients.
But the most prominent feature of butter is that, using the syntax
[___] = butter(___,'s')

with the option 's' indicated (contrast it with default 'z' used for the digital filter design), butter returns the matrix of coefficients for the analog filter design -- which, by definition, the digitalFilter object cannot readily provide.


Dependence of $\tt n$ in $\tt fft(data,n)$ in auto correlation of data The followings are results of autocorrelation of an 'image'(1*) (subplot 221) using xcorr, n point fft, and 2n-1 point fft. I observe that I get different results. Why does it happen?

PS: result at subplot 223 is the expected image.

(1*) In fact it is identified as 'integrated image' in the paper doi:10.1038/nature11578 by Dr. Jacopo Bertolotti, which can not be considered as a general image that has both magnitude and phase information.

Doing an $N \times M$ point FFT for that sized image to calculate the autocorrelation will result in time-aliasing.

To multiply two signals of length $N$ in the frequency domain, the FFT length must be $2N-1$.


Can I calculate a half-zero data's FFT from its non-zero front part's FFT? I have data whose latter half are all-zero. For example, data=[a, b, c, d, 0, 0, 0, 0], and front_half=[a, b, c, d]. Is there any way to obtain FFT(data), by some compute-efficient transform on the FFT(front_half)?

I made some work on the equations. The even items in FFT(data) are just equal with FFT(front_half). But the odd items in FFT(data) seems twisted with dot-multiply with a time-domain weight vector, which looks hard to be calculated from the FFT(front_half).

The computationally efficient way is to do an IFFT/zero-pad/FFT.  It's usually less efficient to do a circular Sinc interpolation (to create the effect of a rectangular window in the other domain) for any large number of extra points.  Any other interpolation is less accurate.


How to predict next image in a sequence, or predict the deblurring of an image? I have 3 versions of an image at lower resolutions - 480p, 720p, and 1080p. Is it possible to use those 3 images to predict the next highest resolution image (the 4k image)?

My first guess was that I could calculate the difference (or momentum or acceleration?) of each pixel between images and use that to predict the new 4th image. I assume this would require me to give all pictures the same number of pixels. I also assume that the differences in each pixel between images would need to be smooth. Would this work?

My other guess was that maybe there is some kind of statistics or machine learning technique that I could use, but I don't know anything about that stuff.

Thanks!

Ps - this man is not my real data. My real images are an intensity map (temperatures over an area) at different resolutions. I suspect that the introduction of colors might complicate things, so should we pretend that the images below are black-and-white images?

Here are the pictures:

There is no extra information in the 720p and 480p images that is not already in the 1080p image. You can do the interpolation on the 1080p image. A recent work that may be of interest to you is Ledig et al. Photo-Realistic Single Image Super-Resolution Using a Generative Adversarial Network, arXiv:1609.04802v5, 15 Sep 2016 (v1), last revised 25 May 2017. But if your work is for scientific purposes rather than for looks, then it is probably best to use traditional interpolation methods that do not try to increase the resolution by making up new data.


Given a set of lines, find only those who are parallel (perspective) I have an image which includes straight borders which I need to detect:

The borders are thick, which causes hough transform to detect multiple lines for each border, as can be seen in the following image:

In order to find the 'true' lines, I want to find all parallel lines (i.e. find the lines which converge into a vanishing point). Does OpenCV (or any other library for the sake of this manner) contains such a functionality? (I'm using python, though I can port c++ scripts as well)

Thanks in advance

Parallel lines in the image do intersect at a vanishing point. Therefore simply hypothesizing lines (a gradient direction at a point suffices to describe it) and voting (see Hough voting) would suffice to identify this point. One could then record all the lines that casted votes to this very point and identify them. Care must be taken as it is difficult to vote directly on a vanishing point, i.e. the space is not bounded and explodes to infinity. Instead, it is possible to use a Gaussian Sphere as Barnard suggested.

Vanishing point extraction is a well studied topic and many algorithms exist, some of which are:

  A Novel Line Space Voting Method for Vanishing-Point Detection of
  General Road Images, Zongsheng Wu, Weiping Fu, Ru Xue and Wen Wang
  https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4970002/

  Vanishing point detection for road detection, Hui Kong, Jean-Yves
  Audibert, Jean Ponce http://www.di.ens.fr/sierra/pdfs/cvpr09c.pdf

  Vanishing point detection in corridors: using hough transform and
  K-means clustering, R. Ebrahimpour,  R. Rasoolinezhad, Z.
  Hajiabolhasani,  M. Ebrahimi
  http://ieeexplore.ieee.org/document/6135447/

  Hough Transform for Directional Orientation, Ivan Papusha, Matthew
  Ho
  https://stacks.stanford.edu/file/druid:rz261ds9725/Ho_Papusha_RealTimeHoughTransform.pdf

Finally, when multiple lines cause a problem, one could always cluster them using the line parameters $[a,b,c]$. When number of lines are not known, one can use a hierarchical clustering, such as single-linkage or agglomerative clustering.


Recommended approach for extracting specific edges from image? I have several sets of images with a typical image being:

and from each image, I would like to extract the top and bottom curves of the object, shown in red here:

What are some possible tools I could use to solve this problem?

Some things I've looked in to:

Edge detection

As seen in the pictures, there are more edges in the image than just the two edges I'm interested in. I've implemented some naive edge detection in Python with scipy, however, this approach returns a lot of edges I don't want and the data is noisy. If I were able to clean up the data, is there a way to automatically select the edges I'm interested in (the red lines in the picture above) from the whole set of edges?

Background Subtraction

The image sets are actually taken from a video of the object moving over time. Using background subtraction, I'm able to create images like this:

But since the object is fixed on the upper-left end, there is little movement and this approach doesn't provide enough data in that region to create a complete line. Additionally, some sets have very little movement anywhere in the image so this approach wouldn't register the object in these cases.

Possibility of using computer vision: An approach that I haven't looked into is implementing some form of computer vision. I see examples of object detection where a bounding box for an object is created. However, I would like to extract specific curves.

Any feedback or help pointing me in the right direction would be greatly appreciated, thanks.

You might find my, or the other answers, from here helpful.

After advice about detecting focus quality of objects in a photo detected using YoloV3

I would have put this in a comment, but I can't include a pic there.

Perhaps you could come up with a better verbal description of how to define the edges of interest.


Separate two measured combinations of two signals with different time delays I have two signal s1(t) and s2(t) that I want to extract but I can only measure:

y1(t) = s1(t) + s2(t-d1)

y2(t) = s1(t-d2) + s2(t-d3)

The time delays d1, d2 and d3 are unknown, though I have a rough estimate for d2.

Maybe (but only if necessary) I could approximate this by shifting y2 so that

y1(t) = s1(t) + s2(t-d1)

y2(t) = s1(t) + s2(t-d3)

Is there any way to separate the signals?
What happens if I have three measurements and three signals?

Elaborating a little bit on what I discussed in the comments: Let's assume the delays $d_1$, $d_2$ and $d_3$ are integer multiples of the sampling time interval $t_0$ so that we can treat it as an integer problem.

Then, as a first step, I would try to use the cross correlation between $y_1$ and $y_2$ to determine the delays. Ideally, it should have a peak at an offset of $d_2$ (when the shifted copies of $s_1(t)$ and $s_1(t-d_2)$ align) and another one at $d_3-d_1$ (when the shifted copies of $s_2(t-d_1)$ and $s_2(t-d_3)$ align).

Now compute $z[n] = y_1[n] - y_2[n-(d_3-d_1)]$ (using $[n]$ as a discrete time / sample index for clarity). If your estimate was correct, this cancels $s_2(t)$ and we are left with $z[n] = s_1[n] - s_1[n-\Delta]$ where $\Delta = d_3-d_1-d_2$.

Now, your signals ought to start somewhere (causality) so that we can assume $s_1[n]=0$ for $n<0$ (if you have different initial conditions, you need to work them in here). Then, we have $$z[n] = s_1[n]-s_1[n-\Delta] = s_1[n]\quad \forall n<\Delta.$$ Therefore $$s_1[n] = \begin{cases} 0 & n<0 \\ z[n] & n<\Delta \\ z[n]+s_1[n-\Delta] & \mbox{otherwise}\end{cases}.$$ This allows you to reconstruct $s_1[n]$ sample by sample.


Discrete time system signal flow graph problem I'm trying to solve the following question, but I've been told I'm doing it incorrectly. Could anyone here give me a hint as to where I'm tripping up? Below follows my current solution:

I do the following to solve it:

a)
$ y[n] = 3 ( x[n] + \frac{1}{3}x[n-1] + 2x[n]) + 6x[n-1] = 3x[n] + x[n-1] +6x[n] + 6x[n-1] = 9x[n] + 7x[n-1]$

b) I find the system function: $H(z) = \frac{Y(z)}{X(z)} = \frac{9X(z)+7z^{-1}X(z)}{X(z)} = 9 + 7 z^{-1}$

Then I apply the inverse z-transform and find the impulse response:
$h[n] = 9\delta[n] + 7\delta[n-1]$

Thanks!

There is clearly recursion in there and your solution does not reflect it since its purely FIR.

I recommend doing this as a coupled difference equation: Define a new variable $v[n]$ at the point where Peter drew the circle and than derive the transfer function from (1) X->V and then (2) X and V -> Y. Then you can eliminate V in eq (2) by using equation (1).


Why does subbing $s = j\omega$ into the Laplace transform of a cosine wave yield a purely imaginary result? The Laplace transform of a cosine starting at $t=0$ is given by

$$F(s) = \frac{s}{s^2 + \omega_0^2}$$

If I sub in $s = j\omega$, I get the Fourier transform of a cosine starting at $t=0$:

$$F(j\omega) = \frac{j\omega}{\omega_0^2 - \omega^2}$$

As this is a purely imaginary result it shows that the input function is made up of only sine waves and must therefore be odd.

This doesn't make sense as a cosine wave starting at $t=0$ is neither odd nor even, so $F(j\omega)$ should have both real and imaginary parts.

I find a similar issue when using a sine wave starting at $t=0$ which returns a purely real function of frequency:

$$F_2(j\omega) = \frac{\omega_0}{\omega_0^2 - \omega^2}$$

What is it that I'm misunderstanding?

If I sub in $s=j\omega$, I get the Fourier transform of a cosine wave starting at $t=0$.

No, you don't. You can't just set $s=j\omega$ in an expression for the Laplace transform and expect the result to be the Fourier transform of the function. This only works if the imaginary axis is inside the region of convergence (ROC) of the Laplace transform. Since the given Laplace transform has poles on the imaginary axis, the imaginary axis is not part of the ROC (the ROC is $\textrm{Re}\{s\}>0$).

The Fourier transform of the function

$$x(t)=\cos(\omega_0t)u(t)\tag{1}$$

exists in a distributional sense, i.e., if we allow Dirac impulses (and, possibly, its derivatives) in the expression of the Fourier transform.

Using the Fourier transform of the step function

$$\mathcal{F}\{u(t)\}=\pi\delta(\omega)+\frac{1}{j\omega}\tag{2}$$

and the modulation property of the Fourier transform

$$\mathcal{F}\left\{e^{j\omega_0t}u(t)\right\}=\pi\delta(\omega-\omega_0)+\frac{1}{j(\omega-\omega_0)}=Y(\omega)\tag{3}$$

we can easily derive the Fourier transform of $(1)$ by noticing that $x(t)$ is the real part of the time-domain function on the left-hand side of  $(3)$, so its Fourier transform must be the even part of $Y(\omega)$:

$$\begin{align}X(\omega)&=\frac12\left[Y(\omega)+Y^*(-\omega)\right]\\&=\frac{\pi}{2}\left[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)\right]+\frac{j\omega}{\omega_0^2-\omega^2}\tag{4}\end{align}$$


"Capacity" concept in Shannon Theorem I have a question about the Shannon Theorem. According to this theorem the maximum capacity of a link can be calculated with the formula:
$$C=B\log_2(1+S/N)$$
where $B$ is the bandwidth and $S/N$ is the signal to noise ratio. But I am not very sure which is the meaning of the "capacity" term in this equation, since it should include the probability of bit error in same way.
What I mean is that the capacity might not be defined as "number of bits transmitted per second" but as something like "number of bits transmitted and correctly received per second". But it still makes no sense, since you can transmit with low $S/N$ and symbol length tending to $0$, getting a capacity that tends to infinite and a probability of bit error close to 0.5 (it is theoretically imposible to have more than 50% of error rate), so the "number of bits transmitted and correctly received per second" would tend to infinite and it would be higher than the Shannon limit. I was thinking that the right definition of the $C$ in the formula could be "number of bits transmitted per second $\times$ (0.5- Probability of bit error)", which could have sense but I am not very sure which is the real meaning of this measure.

The formula you wrote corresponds to the error-free capacity. Namely, Shannon is telling us that $C$ is the maximum (theoretical) limit one can achieve.

Shannon managed to prove that if the information rate on the channel being used is less than $C$, then it is theoretically possible to find a coding system for the signal such that the transmission will have an arbitrarilly small probability of error. Unfortunately, the theorem does not tell us how to find that code, just that it exists.

That's why there is no reference to error rate in the formula for capacity. The formula itself is telling you the maximum information rate you can have with an arbitrary small amount of errors, and that value corresponds to $C$.

Sidenote: remember that the formula is only valid for white Gaussian noise.


What's the impact of aliasing in the time domain? I've been studying digital audio and come across something I can't understand. There appears to be something like a consensus (among those capable of understanding such things) that the impact of aliasing on reconstructed time-domain signals is to ... smear(?) the timing. (Is that the right way to put it?) Here's a quote from a paper by Bob Stuart and Peter Craven, two prominent, respected figures in the audio industry: "Aliasing in the frequency domain is equivalent to the time-domain phenomenon of an impulse response that depends on where, relative to the sampling instants, the original stimulus was presented: see footnote 8." Footnote 8, which refers to a different passage in the article, says, "The complication is that because of the sampling, the total system is not time-translation invariant and so does not have a unique ‘impulse response’ – the response is slightly different according to the position of an original impulse relative to the sampling points."

I did some thinking, and some simple modeling, and came to the conclusion that the effect of aliasing on the reconstructed time-domain signal is to alter the amplitude in a quasi-random way. At the actual sampling instants, the deviation of the reconstructed signal from the original signal resulting from aliasing is zero: At precise sampling instants, aliasing has no effect. But in between those instants there are amplitude errors in the reconstructed signal. This is my conclusion from my own analysis. This is broadly true; I did not do an analysis specifically to determine the effect of aliasing on impulse response.

I cannot see how my conclusion--that of added amplitude noise superimposed on the time-domain signal as a result of aliasing--is equivalent to "the time-domain phenomenon of an impulse response that depends on where, relative to the sampling instants, the original stimulus was presented."

Insights anyone?

I read Footnote 8,

  "The complication is that because of the sampling, the total system
  is not time-translation invariant and so does not have a unique
  ‘impulse response’ – the response is slightly different according to
  the position of an original impulse relative to the sampling points."

as a cumbersome way of saying that aliasing corrputs an LTI system such that its response can no longer be trusted to be linear and time-invariant.

If I were to rephrase things informally, I'd say that:

When a signal is undersampled in the frequency domain or has inadequate frequency resolution (aka freq domain aliasing), its spectrum has overlapping tails.

That means that it is no longer possible to correctly reconstruct the time-domain signal. The spectrum tails do not go to zero, but are folded back. This inversion of the tails are what we know as aliasing and manifest themselves in the time-domain signal as impulse responses (a short-duration time-domain signal) that mix with the original signal.

These mixtures are perceived as distortion/roughness/etc. Their constructive/destructive effect give you your "quasi-random" amplitude modulations.

The magnitude & position of aliasing in the time domain is related to the degree of frequency domain undersampling. See https://ccrma.stanford.edu/~jos/sasp/Example_2_Time_Domain.html

Finally, regarding "time smearing", the impulsive (quickly rising) parts of a waveform contain many high-frequency harmonics that when undersampled in the frequency domain, fold back as described above. However, I wouldn't call this distortion smearing, in the transient-smearing, or pre-ringing sense.

Hope that helps!


Quantized error in scatter plots for a hamming windowed LTE signal I have a time domain signal based on LTE standards.
Carrier spacing = 15KHz, Symbols per slot = 6, extended Cyclic Prefix = 512 samples
FFT size = 2048, No. of sub carriers = 1200, channel BW = 20MHz, Samp Freq = 2048*15KHz = 30.72 MHz, Total Bits = 14400 (Symbols * slots * bits used to MOD * total subcarriers).

It is QPSK modulated (using qammod ("IntegerSymbols",2^bits,binary) in matlab);
This QAM modulated data is IFFT'ed after zero padding on the sides and DC carrier insertion is done and then CP is added (512 samples to a 2048 sample bin).

This signal now is windowed (Point wise multiplied by a window function(only in the TRANSMITTER side) (FOR DIFFERENT REASONS) in time domain) "BY A HAMMING WINDOW FUNTION" Keeping all the overlap add constraint in mind.

On the receiver side, no noise is added, no specific channel is modeled. Simple reception. On the receiver side, CP is removed and then the stream is FFT'ed and only the central info carrying subcarriers are taken. This stream at this point is still modulated. However when I try to plot a "scatter plot figure" of the QAM modulated data from the txer side and the FFT'ed data in the receiver, I get a pattern as shown in the figure above which seems to be some sort of quantized thing. I am not sure what it is or what is it caused by. Could you please help me on it?
On decoding Error/BER is 0 but EVM is -18dB.

P.S: I've asked this ques before but I reckon the details weren't clear. So fingers crossed this time.

That's the picture you'd get usually if you had a single-carrier system with a slight IQ imbalance. (but that would usually be kind of slanted; but we don't know what your OFDM receiver does, so we can't tell you.)

In a multi-carrier system like OFDM, I'd look into whether I accidentally added inter-carrier-interference (ICI), for example by having an uncorrected fine frequency offset. Or, and that's not that uncommon for simulations, by having a control loop trying to achieve frequency synchronization, but goes wrong. (But that would usually be less "nice and orderly looking". However, we don't know what your OFDM receiver does, so we can't tell you.)

Since this is in simulation: Test. Send a stream of constant symbols, from only one subcarrier. That way, you can see where the effect comes from. Iterate. Debug.

You're the engineer sitting in front of that problem – go and investigate!


MatLab - Can not detect all lines in a "simple" image using houghlines I try to detect all lines in an image with some common geometric figures. The obvious choice, hough of MATLAB wouldn't print out the desired results.

I used the following MATLAB code:

[H,T,R]=hough(BW);
P=houghpeaks(H,10,'threshold',ceil(0.2*max(H(:))));

LINES=houghlines(BW,T,R,P,'FillGap',2,'MinLength',3);
figure, imshow(BW),hold on;

for k=1:length(LINES)
    xy=[LINES(k).point1; LINES(k).point2];
    plot(xy(:,1),xy(:,2),'LineWidth',1,'Color','Green');

    plot(xy(1,1),xy(1,2),'x','LineWidth',1,'Color','yellow');
    plot(xy(2,1),xy(2,2),'x','LineWidth',1,'Color','red');
end

EDIT 1: Including the original RGB image here:

And the image stored in the BW:

The result of this code is shown in the next image.

I think it is visible that hough method did not detect every line leaving quiet a few still to be detected.
Does someone have a solution as to why it did not detect those lines and how to optimize the code so that it does detect them?

houghpeaks function in MATLAB accepts a parameter named NHoodSize that specifies the size of the neighborhood where a non-maximum suppression (NMS) is carried out. In simple terms, NMS removes all duplicate peaks which lie in a certain window. This is what explains your double-lines with a single line. If you set the window size lower, such as $[7,7]$, you will see that it will start detecting them. For instance, modify the call as:

  P=houghpeaks(H,20,'threshold',ceil(0.3*max(H(:))), 'NHoodSize',[7 7]);

You will then see the following result:

You should play with those parameters till you find your ideal setting: Not many lines, and correct ones not missed.


Regularization for Inverse Problems using the Singular Value Decomposition (SVD) I am reading these lecture notes on Optimisation and Inverse Problems in Imaging, and I have difficulties understanding how figures on page 20 (Figure 3.2) or page 21 (Figure 3.3).
Precisely, I don't understand what numbers on horizontal and vertical axes mean. I would appreciate if you could explain me this. Here is the code for Figure (3.3).

Similar to The Concepts Behind SVD Based Image Processing the horizontal axis are the samples index of the SVD basis.
The idea in the chapter you linked is generalizing the Wiener Filter.
While the Wiener Filter uses the Fourier Transform as a basis the SVD uses the data adaptive basis.


Is there a way to further optimise my FFT? I've written my own FFT but wanted to know if there was perhaps a way to combine the Cooley and Turkey methods. So far I've got N/2 log N real multiplies for transforming a real signal.

My current implementation has N log N multiplies (well a little more because I've not optimised out the case of the complex sinusoid at 0 radians), but I have an optimisation that reduces the number of multiplies by 50%.

Note: this is absolutely not homework, I've written this algorithm because I found the explanations of the Cooley/Turkey method far too convoluted, and at the same time I wanted the ability to select which bins I processed in both the FFT and IFFT as well as being able to output the IFFT to multiple channels. So I just decided to figure out my own radix-2 while at the same time understanding how and why the transform works (which I have).

ORIGINAL QUESTION (edited for clarity) ---

Example code: kfft.c

I wanted to know whether there is a similar FFT algorithm to mine.

I basically perform FFT by factorising calculations by recursively subtracting the second half of a signal from the first, meaning x[n] only needs to be multiplied by the complex sinusoid for only half of the duration, resulting in N/2 log N real multiplies (shown below). The summed signal is for multiplying with complex sinusoids that are multiples of the sampling period. For example groups {2, 6} are multiples of 2, and therefore have a periodicity of 4 samples.

// For bins { 1, 3, 5 and 7 }

x1[0] = x0[0] - x0[4]
x1[1] = x0[1] - x0[5]
x1[2] = x0[2] - x0[6]
x1[3] = x0[3] - x0[7]
// Calculate sum (required for the next level of bins)
x0[0] = x0[0] + x0[4]
x0[1] = x0[1] + x0[5]
x0[2] = x0[2] + x0[6]
x0[3] = x0[3] + x0[7]

// For bins { 2 and 6 }

x2[0] = x0[0] - x0[2]
x2[1] = x0[1] - x0[3]
// Calculate sum (required for the next level of bins)
x0[2] = x0[0] + x0[2]
x0[3] = x0[1] + x0[3]

// For bin { 4 }
x3[0] = x0[0] - x0[1]
x3[1] = x0[0] + x0[1] // <-- This is the 0hz bin (sum)

Indeed, this is more or less the usual butterfly used in radix-2 FFT (can be extended to any prime number): http://en.wikipedia.org/wiki/Butterfly_diagram


Why would a time shifted signal have a higher maximum frequency than an unshifted one? I´m reading the book Numerical Simulation of Optical Wave Propagation by Jason Schmidt.
In section 2.1-2.2, the author performs the FFT of a centered and a shifted Gaussian pulse and makes the statement that the time shifted signal has a higher maximum frequency than the unshifted one. For that reason more samples are required in getting the shifted signal spectrum.
Why does the time shifted image have a higher highest frequency? I don´t understand.
How much I should increase my sampling rate for a given time shift?
Thanks for your attention.
Carolina Rickenstorff

I was able to take a look at the book.
You are right it says: "The highest significant frequency in the shifted Gaussian signal is higher than that in the centered Gaussian".  I think this is misleading.
Looking at Figs 2.2 and 2.3 I can only think that the author is saying that in Fig 2.3 the phase is "significant" across the whole spectrum, while in Fig 2.2 phase is zero everywhere.
However, the magnitude curves are identical. If magnitude is zero then phase is irrelevant. (i.e. atan (small/small)).  The required sample interval is governed by the magnitude spectrum. That is the required sample interval is not affected by a simple time shift, as stated by the above commenters. If this wasn't the case, then real-data analog-to-digital conversion would be pretty difficult!


How to best separate my signal into frequency regions? I have a signal with frequencies betwwen [0.5,500] Hz and i want to create a filter to best separate the signal in the following regions:

[0.5-3]Hz
[3-8]Hz
[8-13]Hz
[13-30]Hz
[30-500]Hz

My questions are:

How can i find the sampling frequency if i know that the signals are
'downloaded' in Europe (energy transport network of 50Hz)
Whats the most appropriate computational method to separate the signal in
the regions.

I think that in order to find the sampling frequency i need to find the max frequency of my signal and then fsampling=2*fmax.Does the information about the 50Hz network affects my sampling frequency?

Computational effort is probably totally irrelevant considering the low maximum frequency, implying a low required sampling rate, and the fact that this all isn't done in real-time but after the fact (and computers are fast). So, just use your favourite math toolkit (sounds like you'd like to use Matlab, Octave, or Python/Jupyter Notebooks/Scipy), and design five bandpass filters with the characteristics you need, and apply those.

About these characteristics: You don't give us any info about what you want to do with that signal, what that signal is or applicable information about how flat you need your passband to be, how attenuated your stopband or how wide your transition between. So, that's up to you to figure out. I don't think this is out of reach for you – these frequency ranges sound suspiciously much like EEG, so, maybe, read another five pages in a book about such signals.

Now, no, inherently, the grid frequency doesn't have anything to do with the sampling rate. In the presence of strong narrowband (as in: single-tone) interference, it can be advantageous to use a sampling rate that puts intermodulation products into irrelevant bands, but a) you don't have irrelevant bands and b) considering the low highest frequency, you'll probably significantly oversample, anyway, so that this will be inherently be minimized. There's more than just clock/interferer intermodulation, but that'd be a very hardware-specific aspect, and can't be analyzed with much more detailed knowledge about the hardware you're using.

That grid frequency was given to your for another reason. I'll avoid speculating on this here - because honestly, for a "last year exam in DSP", your questions are a bit too basic, and I don't want to take the excitement of getting things right yourself from you.

So, sampling rate: 2·500 Hz is the theoretical lower limit. You can't use that, because your highest filter band needs some transition width. Not that it matters - any cheap soundcard can do 48 kHz or more, and I'd assume whatever device you're using for sampling isn't worse. So: sample first, then reduce bandwidth digitally with decimating filters.

Oversampling (that means, sampling with a rate factors higher than your necessary rate) has a dynamic range and SNR advantage, anyway, so you'd do that in a real system if possible.


Conceptual question on FFT and chirp signal If I take the FFT of a sinusoid I will get a plot whit all the energy of the signal concentrated at the sinusoid frequency. But what happens if I have a signal in which the frequency keeps changing?(like the chirp: linear and exponential)

Is it possible to compute the FFT of a chirp signal? And what should I expect as result?

I tried to do it in Matlab creating a chirp using the matlab function and then passing it to the FFT function but I do not know if this makes any sense and how to 'read ' the result I get.

Could you also suggest me some readings (papers, books) on the topic?

Your question is deeper than you think. For signals whose spectra is time varying,like your chirp signal, you lose this spectral variation information in simple Fourier transforms. You must consult time-frequency analysis to retain the full information of both time variation and spectral variations.

Consider using "specgram()" function of matlab which computes and plots time-frequency spectograms of signals.

Consider the following code from matlab:

  t=[0 0.5 1.0 1.5 2.0]; % time breakpoints
  f=[0 200 100 150 300];  % instantaneous frequency breakpoints
  p=polyfit(t,f,4);      % fit 4th order polynomial over time
  t=0:0.001:2;           % 2 secs @ 1kHz sample rate
  y=chirp(t,p);
  figure,specgram(y,128,1E3,128,120); % 2D time-frequency display

On the contrary a single FFT of the complete signal, which includes all the frequencies, nevertheless, loses their time variation information.


How to characterize this whitening filter? A little bit of background: I am interested in whitening in the least-square estimation.

In this sense, consider a univariate signal $y$.  Assuming the covariance of $y$ is $\mathbf{A}\sigma^2$.  If we knew matrix $\bf A$, then, whitenning of this signal could be done by:
$\mathbf{A}^{-1/2}\bf y$.  Since the  $\mathbf{A}^{-1/2}\bf y$ has a covariance matrix of $\mathbf{A}^{-1/2} \bf A \mathbf{A}^{-1/2}\sigma^2 = I\sigma^2$.  Thus, $\mathbf{A}^{-1/2}\bf y$ is whitened.

In practice, the $A$ is unknown but can be estimated using different methods such maximum likelihood or Yule-Walker.  Here is what I am interested to figure out.

To be more specific, for this purpose one can use Yule-Walker equation in particular and assume an AR(p) model, say order $p=1$ and estimate the model parameter of this autoregressive.

This is indeed a filter, but my question is how to characterize this filter in frequency domain?  How can I find the $H(z)$ of this filter?  Is there a way to find its zeros and poles and further know the characteristics of this filter?

I would appreciate if you could help me.

Hi: The material you are asking about is on page 7 of the document at the link below.

http://www.ee.ic.ac.uk/hp/staff/dmb/courses/DSPDF/DSPDF_2010.pdf


Why do signal have mirror signal after pure frequency shift? In my Matlab code I have a simple simulation of FPGA with the ability to simulate doppler shift, and after I and Q components plugged into IQmod with some carrier freq.

In my Matlab code when I multiply primordial signal $s_{prim}(t)e^{(jwt)}$
and have mirror signal after? And what's strange when I carrying that signal onto carrier frequency this carried signal doesn't have mirror signal.

totally working Matlab code included, will be glad to any ideas.

t = 0:0.000001:0.1;
fs = 1/(t(2) - t(1));
%ss = exp(j*2*pi*100*t) + exp(j*2*pi*200*t) + exp(j*2*pi*300*t) + exp(j*2*pi*900*t) + exp(j*2*pi*600*t); %Generating I-Q signal *.bin file for AWG

N=8;
bit_stream = round(rand(1,N));
f = 1000;
P1= 0;
P2= pi;
PSK_signal = [];

for ii = 1:1:N
PSK_signal = [PSK_signal (bit_stream(ii)==0)*sin(2*pi*f*t + P1)+...
        (bit_stream(ii)==1)*sin(2*pi*f*t + P2)];
end

ss = PSK_signal;
figure(1); plot(ss);

sp = fft(ss);
N = length(ss);
fp = ((0:N-1)/(N-1))*fs;
fp = fp(1:N/2);
sp = abs(sp(1:N/2))/(N/2);
figure(2); hold on; plot(fp,sp); %Spectre of primordial signal

w = 3000*pi*2;
t = (1:length(ss))/fs;

%Real outputs from AWG
I = real(ss); % I Input on FPGA
Q = imag(ss); % Q Input on FPGA

%Adding frequency SHIFT in Simuation FPGA plate
ss_i_real = I.*cos(w*t) - Q.*sin(w*t); % I component from FPGA
ss_i_imag = I.*sin(w*t) + Q.*cos(w*t); % Q component from FPGA

sp = fft(ss_i_real+1j*ss_i_imag);
N = length(ss_i_real);
fp = ((0:N-1)/(N-1))*fs;
fp = fp(1:N/2);
sp = abs(sp(1:N/2))/(N/2);
figure(2);  plot(fp,sp);

w = 40000*pi*2;
t = (1:length(ss_i_real))/fs;

%ss_c = ss_i.*exp(1j*w*t); %Modulation that carrying signal from FPGA on Carrier Frequency (40KHz)
ss_c_real = ss_i_real.*cos(w*t) - ss_i_imag.*sin(w*t);
ss_c_imag = ss_i_real.*sin(w*t) + ss_i_imag.*cos(w*t);

ss_c = ss_c_real+ss_c_imag*1j; % Real part of Signal onto carrier <- this is the signal that will be transmitting from modulator

sp = fft(ss_c);
N = length(ss_c);
fp = ((0:N-1)/(N-1))*fs;
fp = fp(1:N/2);
sp = abs(sp(1:N/2))/(N/2);
figure(2); plot(fp,sp); %Spectre of signal onto carrier freq

% figure(4); plot(ss_c_real); %I - component of signal onto carrier freq
%
% figure(5); plot(ss_c_imag); %Q - component of signal onto carrier freq

edit1:
What causes my confusion, when I'm using the signal:

 %ss = exp(j*2*pi*100*t) + exp(j*2*pi*200*t) + exp(j*2*pi*300*t) + exp(j*2*pi*900*t) + exp(j*2*pi*600*t); %Generating I-Q signal *.bin file for AWG

I have this nice picture:

1 primordial signal, 1 shifted signal and 1 carried onto 40KHz shifted signal.
When I'm using the BPSK signal I'm shifting this signal I have mirrored signal and shifted (i mean two signals with $F_{shift} \pm F_{prim}$)

To be simple I just want to add shift on intermediate frequency and after carry that shifted signal onto carrying frequency without clone signal at $f_{car}-f_{shifted}$ as I have now.

Your code is (probably) ok. But your way of plotting is slightly wrong.

The complete spectrum of the baseband PSK signal would include two peaks one at positive frequency $f_0=1000$ and the other at negative frequency $f=-1000$. When you IQ modulate this signal with $f_m=3000$ you will create two twins one at $f_{11}=2000$ and $f_{12} = 4000$ and the other at the negative frequencies at $f_{21}=-2000$ and $f_{22} = -4000$.

Since you are taking only the first half of the FFT you do not observe the peaks at negative ones, which causes the confusion. If you plot the full spectrums, you can see the missing components at the second half that you omit...

t = 0:0.000001:0.1;
fs = 1/(t(2) - t(1));
%ss = exp(j*2*pi*100*t) + exp(j*2*pi*200*t) + exp(j*2*pi*300*t) + exp(j*2*pi*900*t) + exp(j*2*pi*600*t); %Generating I-Q signal *.bin file for AWG

N=8;
bit_stream = round(rand(1,N));
f = 1000;
P1= 0;
P2= pi;
PSK_signal = [];

for ii = 1:1:N
PSK_signal = [PSK_signal (bit_stream(ii)==0)*sin(2*pi*f*t + P1)+...
        (bit_stream(ii)==1)*sin(2*pi*f*t + P2)];
end

ss = PSK_signal;
%figure(1); plot(ss); title('The PSK baseband signal');

ff = linspace(-fs/2,fs/2,length(ss));
figure,plot(ff,abs(fftshift(fft(ss))));
title('The PSK baseband SPECTRUM ');
% sp = fft(ss);
% N = length(ss);
% fp = ((0:N-1)/(N-1))*fs;
% fp = fp(1:N/2);
% sp = abs(sp(1:N/2))/(N/2);
% figure(2); plot(fp,sp); %Spectre of primordial signal
                     %title('The PSK baseband SPECTRUM ');

w = 3000*pi*2;
t = (1:length(ss))/fs;

%Real outputs from AWG
I = real(ss); % I Input on FPGA
Q = imag(ss); % Q Input on FPGA

%Adding frequency SHIFT in Simuation FPGA plate
ss_i_real = I.*cos(w*t) - Q.*sin(w*t); % I component from FPGA
ss_i_imag = I.*sin(w*t) + Q.*cos(w*t); % Q component from FPGA

sp = fft(ss_i_real+1j*ss_i_imag);  % computing the Spectrum of the ANALYTIC signal here! Hence the NEGATIVE frequencies are deleted.
figure,plot(ff,abs(fftshift(sp)));
title('IQ modulated signal FSK baseband');

% sp = fft(ss_i_real+1j*ss_i_imag);
% N = length(ss_i_real);
% fp = ((0:N-1)/(N-1))*fs;
% fp = fp(1:N/2);
% sp = abs(sp(1:N/2))/(N/2);
%figure(2);  plot(fp,sp);
%figure(3),plot(fp,sp); title('IQ modulated signal FSK baseband');

w = 40000*pi*2;
t = (1:length(ss_i_real))/fs;

%ss_c = ss_i.*exp(1j*w*t); %Modulation that carrying signal from FPGA on Carrier Frequency (40KHz)
ss_c_real = ss_i_real.*cos(w*t) - ss_i_imag.*sin(w*t);
ss_c_imag = ss_i_real.*sin(w*t) + ss_i_imag.*cos(w*t);

ss_c = ss_c_real+ss_c_imag*1j; % Real part of Signal onto carrier <- this is the signal that will be transmitting from modulator

sp = fft(ss_c);
N = length(ss_c);
fp = ((0:N-1)/(N-1))*fs;
fp = fp(1:N/2);
sp = abs(sp(1:N/2))/(N/2);
figure(4); plot(fp,sp); %Spectre of signal onto carrier freq
title('IQ modulated and up-mixed FSK signal');

% figure(4); plot(ss_c_real); %I - component of signal onto carrier freq
%
% figure(5); plot(ss_c_imag); %Q - component of signal onto carrier freq


what is a Gaussian signal? I am learning Wiener filtering for 2D images by myself. From one book, it reads

It should be noted that the study of the more general problem of signal denoising dates
back to at least Norbert Wiener in the 1940s. The celebrated Wiener filter provides the optimal solution to the recovery of Gaussian signals contaminated by AWGN. The derivation
of Wiener filtering, based on the so-called orthogonality principle, represents an elegant solution and the only known situation where constraining to linear solutions does not render
any sacrifice on the performance.

It mentions Gaussian signals. I know Gaussian noise. What is Gaussian signal? I cannot find the definition of Gaussian signal via google. Does it mean one signal is contaminated by Gaussian noise?

The stochastic description of signals is realized by random processes. The book you cite actually speaks of a Gaussian random process. By definition, every random variable drawn from that process has a Gaussian probabilty density function. Mathematically speaking: let $\mathbf X(n)$ be the random process representing the signal, then for every $n_1$ the random variable
$$
X=\mathbf X(n_1)
$$
is Gaussian distributed. This in turn means, that $X$ has the probability density function
$$
p_\mathrm X (x) = \frac{1}{\sqrt{2\pi\sigma^2}}\mathrm e^{-\frac{(x-\mu)^2}{2\sigma^2}},
$$
where $\mu$ and $\sigma$ are the mean and standard deviation of $X$, respectively.


Spectral Effeciency of BPSK The null-to-null bandwidth of the main lobe of a BPSK signal is twice the symbol rate.

BPSK transmits one bit per symbol, so the signal bandwidth is twice the bit rate.

That has been shown on a few stack overflow questions. This implies the spectral efficiency is 0.5 bits/second/Hz.

However, a quick Google for the spectral efficiency of BPSK gives an answer of 1 bits/second/Hz. Why?

The spectral efficiency depends on the pulse shape. The basepand BPSK signal can be written as $$s(t)=\sum_k a_k p(t-kT_b),$$ where $a_k$ is equal to either $\sqrt{E_b}$ or $-\sqrt{E_b}$, $E_b$ is the bit energy, $T_b$ is the bit interval (so that the bit rate is $R_b=1/T_b$), and $p(t)$ is a Nyquist pulse. The bandwidth of $s(t)$ is equal to the bandwidth of $p(t)$.

Unfortunately, many textbooks insist on limiting the discussion to square pulses, where $$p(t)=\begin{cases}1,\text{ if $0\leq t <T_b$}\\0,\text{ otherwise.}\end{cases}$$ This means that the spectrum of $s(t)$ is infinite, strictly speaking. This would seem to preclude any discussion of spectral efficiency.

However, when bandwidth is understood as an engineering (not mathematical) concept, we see that we can truncate the spectrum to a desired bandwidth $B$ and, if the pulse $p(t)$ is not distorted too much, communication is still possible.

Different authors make different choices for $B$. As you've seen, some say $B=2R_b$ and others say $B=R_b$. Personally, I use $B=5R_b$ (I'm not as optimistic as most book authors). Fortunately, you can make up your mind: using Matlab (or similar), you can create a close-to-ideal BPSK signal $s(t)$ by sampling it much faster than required, and then filter it using different values of $B$. You can see the resulting distortion and intersymbol interference, and decide what value of $B$ you're comfortable with.

If you have the time and inclination, you can do the experiment with an analog signal generator and analog low-pass filters; it's more work but it feels more "real" than a Matlab simulation.


Practical application of complex FFT Is there any application where one may need to compute the FFT of a complex-valued signal? I am interested in cases where the FFT of the complex-valued signal is of primary interest and is not an intermediate step in a larger computation, so cases like using analytic signals in the computation of the Pseudo-Wigner-Ville distribution is not of interest to me.

An SDR app (for example one for an RTL-SDR device plugged into a Raspberry Pi's USB port) which displays the RF spectrum of an IQ baseband signal uses the FFT of a complex input.  This allows showing the RF spectrum both above and below the PLL frequency.