spectral factorization method So I'm trying to figure out a method to do spectral factorization, where we take some function $X(D) = (1+cD)(1+cD^{-1})$. I believe this is used in equaliation I can do this if I just plugging in numbers until i get something that works, but that doesn't seem like the best way to go about it. Lets look at the following example:

$$X(D) = \frac{.19}{(1+.9D)(1+.9D^{-1})}+\frac{1}{10} = \frac{.371+.09D^{-1} +.09D}{(1+.9D)(1+.9D^{-1})}$$

The bottom part is easy, because it's already factorized into the form we want, so it's the top part that's the hard part. I can come up with the following answer but trying random numbers until I start getting close to what I want:

$$ X(D) =  \frac{.3477(1+.2588D)(1+.2588D^{-1})}{(1+.9D)(1+.9D^{-1})} = cY(D)Y^*(D^{-*})$$
where $Y(D) = \frac{(1+.2588D)}{(1+.9D)}$ and $c= .3477$.

I've tried searching for a method of doing this, but I just can't find one that makes sense and that works every time.

I've tried getting it in the form of $D(D^2 + bD + c)$ where the $D^2$ term has a coefficient of 1, and $b$ and $c$ are constants, and then try factoring that into two parts, but it just doesn't seem to be working out for me. I'd appreciate any advice!

The question is not so clear. I am not sure what other types of examples you consider. In this particular case, if you didn't get the hint notice that:

$$\begin{align}
c(1+\alpha D)(1+\alpha D^{-1}) &=  (c+\alpha^2c)+\alpha cD^{-1}+\alpha cD \\&=0.371+0.09D^{-1} +.09D
\end{align}$$

which gives you
$$\left\{\begin{matrix}
c(1+\alpha^2) = 0.371\\
\alpha c=0.09
\end{matrix}\right.
\Rightarrow \frac{1+\alpha^2}{\alpha}=\frac{0.37}{0.09}$$
and you will actually have two solutions

$$\left\{\begin{matrix}
\alpha=3.8634, c=    0.0233\\
\alpha=0.2588, c=  0.3477
\end{matrix}\right.
$$

See if the second answer is also valid...


Acoustic Echo Cancellation I have the following diagram for acoustic echo cancellation

I'm having a hard time figuring out what's present on the top and bottom lines in terms of reflected echos, original signal etc. Here's what I think

I am assuming a situation in which I am talking in my microphone at home via some software to a meeting room in which a teleconfrencing setup is present. The microphone and speaker in the diagram repersent the teleconfrencing setup. My voice goes through the teleconfrencing speakers, hit's the walls of the meeting room and back into the microphone.

So, what I think is that the right after the microphone symbol in the diagram is my echoed voice x_echo[n]
On the bottom line at the speakers, only my original voice should be present since that's what's meant to be going through the speakers - so x[n].
Now, if that's the desired signal for my LMS adaptive filter, then how would my LMS adpative filter work, there would be no correlation at all.

The acoustic echo canceler (AEC) is supposed to cancel any signal in the microphone path that is correlated to the speaker signal. The microphone doesn't only pick up your voice but it also picks up the signal from the speaker, i.e., the signal from the far end. If there were no AEC, that signal from the speaker would be fed back to the far end, which is very disturbing. Since the speaker signal is available at the near end, it can be used as the input to the adaptive filter, which tries to remove any signal in the outgoing path that is correlated with the speaker signal. Basically, the adaptive filter will converge to the room impulse response between loudspeaker and microphone.


Confusions regarding meaning of "kernel" in signal processing and image processing? While studying filters in DSP and DIP, the term "kernel" is often encountered
What is the meaning of "kernel" in context of filters??Please kindly explain in simple words with an example

I have also attached a snapshot from gonzalez

What is meant by "neighborhood of operation" here?

and in last line,author says that mask,template & window terms are synonyms for kernel.

Is it true even in case of FIR windows such as Hanning window,Blackman window etc? Those windows are also kernels?

What is the meaning of "kernel" in context of filters?

A "kernel" is the thing you convolve by.  So for a FIR filter that's defined by the taps $k = \begin{bmatrix}a_0 & a_1 & a_2 & \cdots\end{bmatrix}$, $k$ is the kernel.

For a 2D filtering problem, the kernel would, in general, be 2D (and you tend to see the term much more often in image processing, but I've also seen it in really highfalutin' stochastic signal processing literature).

  What is meant by "neighborhood of operation" here?

For an output pixel at position $\begin{bmatrix}x, y\end{bmatrix}$, all of the pixels in the input image that affect the output.  In general, if you have an $n \times m$ kernel, the neighborhood of operation is $n \times m$.

  Is it true even in case of FIR windows such as Hanning window,Blackman window etc? Those windows are also kernels?

Different disciplines, different terminologies.  No, because those don't define filter kernels -- the "windowing" operation is different from the "windowing" in 2D image processing.


What is the difference between the PSD of a deterministic and stochastic signal? I am learning about stochastic processes and I don't get one thing:

What is the advantage of calculating the PSD of a signal using the Wiener-Khinchin theorem $\Phi(\omega) =\mathcal{F}\{R_{xx}\}$ and calculating its PSD with $\Phi(\omega) =X(\omega) X^*(\omega)$.

I mean, if I record some random signal and I want to get its spectrum then I do $\mathcal{F}\{x(t)\}$.

What's wrong about using this way?
Why should I deal with autocorrelation and Wiener-Khinchin-Theorem?

I've created in MATLAB a cosine with random amplitude, so in my understanding this is a stochastic signal

t = 0:0.1:1000
x = rand() * cos(t)
plot(abs(fft(xcorr(x,x))))
plot(abs(fft(x)).^2)

Both plots seem to be very equal, I can't see the difference.

Generally, one cannot compute the Fourier transform of any random process, because of integrability reasons. The PSD is one way to define it, based on autocorrelation.

And now, this mathematically sound definition of a periodogram is consistent in some way: if you have a realization of a random process, on a limited time interval, then its Fourier transform can be computed, and is consistent with the PSD, which you'd expect.


Conditions for unique (not exact) recovery from noisy compressed measurements I am looking for theory on whether compressed sensing reconstruction via ℓ1-minimization is unique and under which conditions.
I have looked through:

  Tropp, J. A., "Just relax: convex programming methods for identifying sparse signals in noise," IEEE Transactions on Information Theory, 2006, 52, 1030-1051.

This paper states that the minimizer is unique but I can't quite boil down exactly what is required for this uniqueness (Eq. (ℓ1-Error) and Theorem 14) to hold.

I also looked at:

  Candès, E. J.; Romberg, J. & Tao, T., "Stable signal recovery from incomplete and inaccurate measurements," Communications on Pure and Applied Mathematics, 2006, 59, 1207-1223.

However, they do not seem to claim that the minimizer is unique.
Since both of these papers are from the earlier days of compressed sensing, I suspect that there may be newer results on this uniqueness of the solution. Do any of you have some hints?

"Recovery of Exact Sparse Representations in the Presence of Bounded Noise", by J. Fuchs, deals with the case you ask. From the abstract:

  The purpose of this contribution is to extend some recent results on sparse representations of signals in redundant bases developed in the noise-free case to the case of noisy observations. [..] We consider the case $b = Ax_0 + e$ [...] and seek conditions under which $x_0$ can be recovered from $b$ ...

I think theorem 2 and 3 are what you're looking for (at least in part).


Fourier Transforms and Series for the NON mathematically inclined. This would most likely be the opposite of this question ( Mathematically inclined Signal and Systems/Signal Processing book? ) I figured I'd ask here if there are any good books that while, necessarily, focus on the math are more oriented to what these concepts can DO instead of the mathematics behind it, so I guess more how it changes the graphing or the function etc instead of the functions behind them.

The book "Blip, Ping & Buzz", by M. Denny (2007, John Hopkins University Press), explains signal processing, as used in sonar and radar, at the level of a supermarket science magazine.  Some algebra required.  More on the use of Fourier analysis than any theory.


Is the typical implementation of low pass filter in C code actually not a typical low pass filter? Quite often in C code you see something along the lines of:

$$ y = (1-a) y_{prev} + input * a $$

..where $a$ is some small number. This equation then runs in discrete intervals at the loop rate.

Let's assume $a = 0.1$.

If I convert this difference equation into z domain I arrive at the following:

$$ Y(z) = 0.1 / (1 - 0.9 * z^{-1}) = 0.1z/(z - 0.9)$$

This function has a horrific phase response:

If I instead start with a continuous time low pass filter and convert it to difference equation I get a slightly different filter.

$$ H(s) = 1 / (s + 1)$$

I convert it into z domain using bi-linear transform:

$$ H(z) = (0.04762 + 0.04762 z^{-1})/(1 - 0.9048z^{-1}) $$

I then convert this into difference equation:

$$ y_n = 0.04762x_n + 0.04762x_{n-1} + 0.9048 * y_{n-1} $$

The result compared to the "simplistic" filter:

The result of the "proper" C implementation looks much more like a low pass filter than the "simplistic" implementation. The simplistic implementation is more like a lag regulator (I realized this while writing this post so it's almost that I have answered my own question). While the "proper" lpf is more like an analog LPF. Obviously this can negatively impact the performance of the system if the different behaviors are not taken into account. Are there situations when either one is better suited than the other?

The "simplistic" implementation that you referenced above is often called a leaky integrator, a special case of a first-order IIR lowpass filter. They are discussed in more detail in several previous questions:

Is a leaky integrator the same thing as a low pass filter?
Is there a technical term for this simple method of smoothing out a signal?
Designing a practical integrator
"Low-pass filter" in non-EE, software API contexts

To answer your question specifically, yes, it is a lowpass filter. You can certainly construct more complicated filters with a better lowpass frequency response. The main advantages of the leaky integrator construct are its very simple implementation, and the fact that you can adjust its cutoff frequency by simply changing the value of $a$ (as $a \to 0$, the cutoff of the filter decreases).


Fourier antitransform using scaling property? I'm trying to calculate the antitransform of:

$\frac{1}{2\cdot(1+5w)^2}$

Now I know the antitransform of $\frac{1}{(1+5w)^2} = t \cdot e^{-5t} u(t)  $

But in this case I got that divided by 2. I assumed I had to use the scaling property which says:

$F[f(ax)] = \frac{1}{|a|} \hat{f}(\frac{w}{a})$

Now I'm not really sure how to apply this. Could anyone help?

If $h(t)$ is the inverse Fourier transform of $H(\omega)$, then by linearity the inverse Fourier transform of $aH(\omega)$ is simply $ah(t)$. This has nothing to do with the scaling property you mentioned, because the latter refers to the scaling of the argument of the function.


Matlab: Difference between flat fading and Rayleigh multipath fading, terminologies and concept The information that I put down are from the site http://www.gaussianwaves.com/2011/05/ebn0-vs-ber-for-bpsk-over-rayleigh-channel-and-awgn-channel-2/ . My QUestions are due to some misconception and will really appreciate if some help is provided to clear the concepts.

AWGN channel:

If $y$ is the received data and $n$ is the input, then signal is received  as $y = x + n$. This is an FIR model.

Rayleigh channel

The received signal here is $y = hx +n$ where $h$ can be $h_1$ or $h_2$. This is also an FIR model.

Case1: When $h = h_1$:

h_1 = 1/sqrt(2)*[randn(1,N) + j*randn(1,N)]; % Rayleigh channel creates the coefficients $h$ of length $N$ representing the channel, where $N$ is the number of data samples. This is a flat fading Rayleigh channel because the number of taps = 1. Please correct me if I am wrong in saying that the number of taps =1, so it is a flat fading. In general, number of taps = number of delays.

Case2: When $h=h_2$

If $h_2 = [1, 0, 0.5, -0.2]$ then also we call it Rayleigh but now it is also called multi-path. So, number of multipaths = length(h_2) = 4 ?

How is $h_2$ a Rayleigh random variable? What is the proper way to represent $h_2$ as rayleigh channel and the output $y$? Explanation with the help of code statement will be useful.
Theory says that there is Clarke's and Young's model for channel representation. Clarke's model uses FIR representation using sinusoids and cosines. Why do we need this model when we can simply apply the technique used to generate $h_1$ ?

The fundamental idea to keep in mind is that in a wireless channels with reflections, if you transmit $s(t),$ you'll receive $$r(t)=\sum_{i=1}^Na_is(t-\tau_i).$$

Another important idea is that whether the channel is flat or not depends only on $s(t)$.

For instance, let's say that the symbol time is $T_s$. Let's call the longest delay in the channel $\tau_{L}$. If $\tau_L > T_s$, then one symbol will appear at the receiver "on top" of other symbols. This is called "intersymbol interference" and the channel is called a "wideband" channel or a "frequency-selective" channel.

On the other hand, if $\tau_L << T_s$ (that is, if the longest delay is much shorter than the symbol time), then there is no intersymbol interference. However, each symbol is still "interfering with itself": many replicas of the same symbol, each with different amplitude $a_i$, are received at the same time. If $N$ is large, then you can invoke the central limit theorem and claim that the overall effect of all these replicas is equivalent to multiplication by a single Gaussian random number. This channel is called "flat".

Note: the same channel can behave as a "wideband" channel or as a "flat" channel to two signals with different symbol rates.

Now, to the question about how to simulate this. It all depends on the sampling frequency used to sample the channel impulse response. When you say that a channel model consists of a single random number, what you're saying is that you're sampling the channel impulse response at the symbol rate. In other words, the channel gain $h_1$ is the aggregate effect of all the $a_i$ and $\tau_i$.

Even for a flat channel, if you sample fast enough, you'll be able to resolve all the $a_i$ and $\tau_i$.

In your question, vector $h_2$ can only be interpreted properly if you know the sampling rate. If the sampling rate is the symbol rate, then the channel's longest delay is of the order of 4 symbol times, and you'll have quite severe intersymbol interference. If the sampling rate is much faster than the symbol rate, the channel is flat.


Oversampling and Aliasing I am really confused about oversampling and upsampling. I know that upsampling means stretching our signal in time domain by a factor of K, then interpolating some values or using just zeros. As I found out stretching in time domain is shrinking in frequency domain by the same rate. For example, stretch by factor of 2 in time domain mean shrinking by factor of 2 in frequency domain.

Let talk about oversampling, we just increase the sampling frequency more than Nyquist rate. Then we have too much information. Is there any chance to have aliasing? I mean in time domain, lets say I am sampling my signal (frequency is 70 Hz) by sampling frequency equals to 280 Hz (Nyquist*2) so, is there aliasing in frequency domain? Sorry if my explanation is not clear. It is because I am bewildered with these concepts; oversampling, upsampling, downsampling, undersampling, critical sampling and aliasing.

theoretically, I know them but when it comes to practical view, I have serious problem.

When upsampling, you don't really stretch the signal in time. You insert new samples between the existing ones, without modifying the times at which those samples were taken. One property of upsampling is that the waveform remains exactly the same before and after the process; you just increase the number of samples. Note that the main difference between oversampling and upsampling is that the former occurs at the time of sampling, and the latter occurs after sampling has already been done. If there is no aliasing, in theory both produce the same result.

Aliasing is only present when there is a signal at the input of your sampler whose frequency is higher than the Nyquist frequency. In your example, if you sample a 70 Hz signal at 280 samples per second, you will not have aliasing. However, in practice you will not always have precise knowledge or control over the signal you're sampling. For example, you may want to sample a signal coming from an antenna. You don't know in advance what signals are going to be picked up by the antenna. In a case like this, one approach is to low-pass filter the signal before sampling, to ensure no aliasing will happen.


Using fft to calculate vsdft We are currently implementing the velocity synchronous discrete Fourier transform (vsdft) for order
tracking rotating machinery in Python.

The formula is similar to dft. Simplified, it is given by

$$ \text{VSDFT}(\Omega) = C\sum_{n=0}^{N-1} x[n] \omega[n] e^{j\Omega \theta[n]} $$

where $C$ is a constant and

$$ \omega = \frac{d\theta}{dt} $$

Since dft is given by

$$ F_k = \sum_{n=0}^{N-1} f_n\cdot e^{-i 2 \pi k n / N} $$

and we wanted to use the optimized fft implementing of numpy, and was thinking if there is a way to write the input f so that we could
use fft directly. However, my math skills are failing me.

Can any of you see a way or any pointer to implement the above VSDFT without having to reimplement the Cooley-Tukey FFT algorithm?

Refrence to the original paper can be found in here.

I believe the short answer is no.

The paper is behind a paywall so I didn't look at it.  In the abstract is the following blurb: " use interpolation techniques to resample the signal at constant angular increments, followed by a common discrete Fourier transform to shift from the angular domain to the order domain."  In which they are describing other techniques.

Looking at your equations, that is precisely what needs to be done to $\theta[n]$ in order to convert it into a regular DFT, suitable for a FFT.

The FFT exploits symmetries in the even spacing of the exponents in order to save computations.  With a varying $\theta[n]$, those symmetries aren't there so you aren't going to be able to implement a FFT algorithm.

You don't tell us what $x[n]$ is either.

How large is N?

Hope this helps at least a little.

Ced


Can anyone help with cross correlation error in Simulink? I am sorry if I say any thing silly. I have a vector x below and trying to implement correlation with itself (lets say x2) in Simulink, but getting error. I am getting all zeros in result of my output vector. Kindly can some one please help me out in it. Thanks in advance.

ts=1*10^-9;
t=0*10^-9:ts:150*10^-9;
fc=1*10^6;
yi=154*10^-9:ts:2000*10^-9;
y=yi*0;
xi=2*sin(2*pi*fc*t);
x=[xi y];
x2=x;
z=xcorr(x2,x);

The edited code in your question should not result in an all-zero vector. Of course, since x contains many zeros, the autocorrelation also has many zeros, but it's non-zero in its center portion (i.e., around index $2000$).

Also note that you don't need to define a second identical vector to compute the autocorrelation. You can use the command z = xcorr(x);


How to get the DCT energy distribution image in MATLAB,or we may call it the DCT coefficients distribution image? The picture below is drawn by the DCT coefficients after making DCT to an image :

I want to get an image like the picture above. Here are some problems:

Should I do the DCT to the whole image instead of using the block DCT? And what MATLAB function should I use?
How to deal with the coefficients? Make it in the range of $[0, 255]$? Or make the coefficients matrix to be a Binary Image?

I would really appreciate if someone can help me!

The general shape of this  image can (probably) be obtained with Matlab function DCT2 applied on the whole image (problem 1). I suspect this  image comes from the absolute of the  actual DCT matrix which generally has positive and negative values (problem 2). The absolute value was probably modified by a non-linear mapping, such as a logarithm or a power. The reason is you almost mostly black values on the bottom-right, probably zeroes at high frequencies. Would you have positive and negative values, this area could appear in mid-gray colors.

The PNG file has 251 colors, so the  absolute DCT matrix  with real values was somehow rounded to 8 bits, possibly scaled, and not binarized.

On your question: How to deal with the coefficients? you should tell us what you would like to do apart from displaying  a grayscale absolute 2D DCT image. You can compute energy concentration in the DCT domain, or identify the two bright lines as two main orientations with high gradients on the original image.

The following code produces a similar image, depicted below.

parBoxCox = 0.15;
data = imread('gantrycrane.png');
dataTrans = (1+(abs(dct2(rgb2gray(data))))).^parBoxCox;
dataTrans = floor(dataTrans*255/max(dataTrans(:)));
imagesc(dataTrans);colormap gray


What are the pros and cons of wavelet for filtering compared to conventional filters? Wavelets have been widely used in denoising or extracting one specific frequency band of a signal nowadays. However, these can also be done through conventional filters (e.g. butterworth, Chebyshev). So what are the pros and cons for these two methods for filtering?

Let us first consider the orthogonal DWT: this transform is build with several constraints: orthogonality (and invertibility) of course, and the discretization of some continuous wavelet transform, with a specific dyadic structure, yielding a frequency decomposition whose cut are approximately:

$$[0\;1/2^L\;1/2^{L-1}\ldots 1/2 \;1]$$

Discrete wavelets can be implemented as a bank of filters (or a filterbank) dedicated  to one of the above intervals, followed by the appropriate subsampling. And only some specific filters (albethey in infinite number) obey the above constraints.

Hence, with the DWT, only dyadic boundaries are available (*), only  a subset of filters can be used, and outputs are decimated. And orthogonality, scales are not used at all.

Whereas with standard filters, you only have linearity constraints: invertibility is not required, you can specify the bounds, ripples and decay you need (at the price of design issues).

So for mere filtering, I cannot see an advantage of using discrete and critical wavelet schemes. However, DWT can do some crude band-pass filtering while being used for other processing needs.

Wavelets now play a more interesting role for data sparsification, along with non-linear analysis, for compression, restoration, which standard filters generally cannot achieve.

As a side note, if you allow redundancy in the transformation, choosing a wavelet with a sufficient number of oscillations in a sufficiently long "window" could do a good job, but that would be overkill.

(*) with $M$-band wavelets and wavelet packets, you could get more generic $M$-ary integers $k/M^n$, but not all fractions or reals


Continuum Fitting with Python using specutils it's my first time processing spectra from a black body radiation experiment, I'm using Python and having some troubles...
I have this spectra with 2 peaks and uneven background noise which I want to normalize, I tried using specutils 1.1 "continuum-fitting" (documentation here: https://specutils.readthedocs.io/en/stable/fitting.html#continuum-fitting )
My implementation is this:
import matplotlib.pyplot as plt import numpy as np
from astropy.modeling import models
from astropy import units as u
from specutils.spectra import Spectrum1D, SpectralRegion
from specutils.fitting import fit_generic_continuum

...

x = data[:, 0]
y = data[:, 1]
plt.plot(x, y, label = 'My Data', c='C0')

spectrum = Spectrum1D(flux=y*u.Jy, spectral_axis=x*u.um)
g1_fit = fit_generic_continuum(spectrum)
y_fit = g1_fit(x*u.um)

plt.plot(x, y_fit, label = ' Specutils Continuum Fit', c='C1')

plt.legend()

But the result is somewhat disappointing

Am I doing something wrong? Do you have any other way to normalize the background noise in python? any suggestion at all?

At the end I settled on using a median filter from scipy (medfilt) and with a high enough kernel it removed the spikes completely.
The problem on my specutils implementation was that I needed to exclude the peaks in the spectrum window by doing this, althought I found this too hard to implement on 100-120 data measurement I had:
...
from specutils import SpectralRegion

spectrum = Spectrum1D(flux=y*u.Jy, spectral_axis=x*u.um)
g1_fit = fit_generic_continuum(spectrum, exclude_regions=[SpectralRegion(2 * u.um, 6 * u.um), SpectralRegion(58 * u.um, 65 * u.um)])
y_fit = g1_fit(x*u.um)


Downsampling/Decimation - When is well done? I am working on decimation of signal, and I want to know which is the best way to understand if the downsampling is well done or not. At the moment I am comparing the FFT of the source signal and the downsampled signal and I observed a downward shift of it (I think it is due to the lesser quantity of samples), I had also a look into the time behavior. Do you have other advices/suggestions?

EDIT

From the comments/answers:

one method to understand if the downsampling process is well done is by comparing the FFT of the decimated signal and the original signal. The original signal must be attenuated by the decimation factor as suggested in the answer from @user2934229. Otherwise the signals can be normalized by the length of the signal (not suggested in my case because the signal is too long);

Ignoring your original question, and referring to your comment:

  I am observing the FFT of the downsampled signal presents an attenuation. I don't understand why.

The output of FFT (or DFT in general) isn't normalized - it is relative to the signal length (the total amount of samples being transformed). When decimating a signal (say, by a factor of M), you end up with a signal with 1/M length the original signal's length.

So, because DFT/FFT results are relative to signal length, a decimated signal with a 1/M length would have DFT/FFT components "attenuated" by a factor of M as well, even if the original signal was properly downsampled (i.e. bandlimited before being decimated).

To properly compare between the DFT of different signals, you should divide the output of DFT by the length of the transformed signal, to normalize the transform.


Applying duality property to fourier transform of unit step function For Continuous time aperiodic signals, the duality property of Continuous Time Fourier Transform (CTFT) is following

$$\mathscr{F}\Big\{x(t)\Big\} = X(f), \qquad\text{then} \quad \mathscr{F}\Big\{X(t)\Big\} = x(-f)$$

Now we know while Dirichlet conditions are not satisfied for unit step function $u(t)$, so its CTFT analysis and synthesis cannot be done. However, we can still do it provided we are willing to accept occurrence of singularity functions like Dirac delta impulse in its Fourier transform equation.

i.e.

$$\begin{align}
\mathscr{F}\Big\{u(t)\Big\} &= \mathscr{F}\Big\{\tfrac{1}{2} + \tfrac{1}{2}\operatorname{sgn}(t) \Big\} \\
 &= \frac{\delta(f)}{2} + \frac{1}{j2\pi f} \\
\end{align}$$

where the signum function,

$$\operatorname{sgn}(t) \triangleq \begin{cases}
-1 \qquad & t<0 \\
0 \qquad & t=0 \\
+1 \qquad & t>0 \\
\end{cases}$$

However if I apply duality property to the above result, then I should get following:

$$\begin{align}
\mathscr{F}\Big\{\frac{\delta(t)}{2} + \frac{1}{j2\pi t}\Big\} &= u(-f) \\
 &= \tfrac{1}{2} + \tfrac{1}{2}\operatorname{sgn}(-f) \\
 &= \tfrac{1}{2} - \tfrac{1}{2}\operatorname{sgn}(f) \\
\end{align}$$

However when I read at least some books on Fourier transform, I find that result is $u(f)$ and not $u(-f)$. Question is why ?

The result you got is correct and it is also expected according to the two first formulas in your question. If $X(f)$ is the Fourier transform of $x(t)$, then the Fourier transform of $X(t)$ equals $x(-f)$. If $x(t)=u(t)$ is the unit step function then

$$X(f)=\frac12 \delta(f)+\frac{1}{j2\pi f}\tag{1}$$

and the Fourier transform of $X(t)$ is given by

$$\mathcal{F}\left\{\frac12 \delta(t)+\frac{1}{j2\pi t}\right\}=\frac12+\frac{1}{2 j}\mathcal{F}\left\{\frac{1}{\pi t}\right\}\tag{2}$$

where we recognize $1/\pi t$ as the impulse response of an ideal Hilbert transformer, the Fourier transform of which is given by

$$\mathcal{F}\left\{\frac{1}{\pi t}\right\}=-j\;\textrm{sgn}(f)\tag{3}$$

Combining $(2)$ and $(3)$ gives

$$\mathcal{F}\left\{X(t)\right\}=\frac12\big(1-\textrm{sgn}(f)\big)=u(-f)\tag{4}$$

just as expected.

Maybe you can clarify which books say otherwise.


Even and odd components of a constant function Let $x(t)= 1$ for all values of $t$. The even and odd components of $x(t)$ are ____ , _____
respectively.

You may ask yourself first whether the function is already even or odd. If not, check whether it contains an obvious odd/even part, and then whether the remainder is itself even/odd.
Finally, you can use the hammer of the parity decomposition. A function $f(t) $ can be split into the even and odd parts:
$$f(t) =\frac{1} {2} \big[f(t)+f(-t)\big]+ \frac{1} {2} \big[f(t)-f(-t)\big] $$
but it is often better to look at the function first.


BIBO Stability and the convergence of the frequency response of a system It is my understanding that an LTI system is BIBO stable if and only if its impulse response $h(t)$ is absolutely integrable. This also happens to be one of the Dirichlet conditions for the convergence of the Fourier Transform, from which we can conclude that $H(j\omega)$, the frequency response of the system (i.e. the Fourier transform of the impulse response) converges if $h(t)$ is absolutely integrable.

Now think of an ideal low-pass filter. Its impulse response $h(t)$ is a $\mathrm{sinc}$ function which is not absolutely integrable. In other words, it is not BIBO stable. However, it does have a Fourier transform despite violating the condition of BIBO stability (or equivalently one of the Dirichlet conditions). Why is that? Is it because the Dirichlet conditions are a set of sufficient but not necessary conditions? Or is it because the Fourier transform of a $\mathrm{sinc}$ function doesn't "exist" in the usual sense unless we allow step functions?

The same happens with sinusoids: they have a Fourier transform even though they are not absolutely summable in the time domain. These are cases in which the Fourier transform doesn't converge uniformly, but it does converge in some other sense.

The signal $x(n)=\mathrm{sinc}(n)$ is not absolutely summable... but it is square-summable (i.e. it has finite energy). Therefore, the Fourier transform will converge but in the mean-square sense. In these cases, the transform converges to a discontinuous function (in the case of this specific signal, to a rectangular pulse).

Similarly, sequences that involve sinusoids like $x(n)=\sin(\omega_0n)$ are neither absolutely nor square summable, but a Fourier transform can be defined using distributions (Dirac delta functions).

In both cases the Fourier transforms are not continuous. They are not infinitely differentiable either.

As a general rule we can state that if some impulse response $h(n)\in l^1(\mathbb{Z})$, then its Fourier transform will be continuous. If the Fourier transform is not continuous, then you know that its inverse is not in the sequence space $l^1$.

You mentioned BIBO stability, so I think this may be of interest. There is a theorem that states:

  An LTI system with impulse response $h(n)$ is BIBO stable if and only if $h(n)\in l^1(\mathbb{Z})$.

So one can judge the stability of a system by looking at the continuity of its transfer function.


Similarity in absolute value of nearby DFT coefficients Is there any approximate relation in FFT coefficients?

For example, consider time domain signal $x_n$. Its $N$ point (assume $N$ to be even) DFT is given as $X_0,X_1,...,X_{N-1}$

Can I say something about the nature of $abs(X_1)/abs(X_0)$? I am trying to see the variation in DFT coeffcients and for that I thught one can sonsider the ratio of successive coefficients. My motivation is to see how similar the magnitude spectrogram at frequency $f1$ is to spectrogram at frequency $f1+1$ and $f1-1$

I tried to simplify the expanded equation but got stuck after a step -

$abs(X_0) = \Sigma_0^{N-1} x_n$

$abs(X_1) = abs(x_0+x_{N/2}+2(x_1.e^{-j.2\pi/N} + x_2.e^{-j.2\pi.2/N}+..+x_{N/2-1}.e^{-j.2\pi.(N/2-1)/N}))$

How do I (approximately) simplify $abs(X_1)/abs(X_0)$?

Because the basis vectors for adjacent DFT bins are absolutely orthogonal, the values in the adjacent DFT result bins can be completely unrelated, and thus their ratio can be anything, depending on the data.

If you use a non-rectangular window function, then there will be some leakage of spectral peaks into adjacent bins due to convolution with the transform of the window.  But this depends on the window chosen.


Causal Signal - Fourier Transform or Laplace Transform I am dealing with a physics problem which is related to signal processing. The problem requires me to calculate the instantaneous force acting on a body which depends on some physical parameter $x$. Assume that $x(t)$ is periodic in time for the moment. Since $x(t)$ is periodic, then it can be expanded as a Fourier series with different frequency components (and it doesn't really matter if $x(t)$ is causal). The calculation for the instantaneous force involves adding a complex phase shift (which may depend on the frequency) to each of the frequency component. To do that, I can use the convolution theorem and take the convolution of $x(t)$ with some kernel $\kappa(t)$ whose Fourier transform gives me the required phase shifts, i.e. $\tilde{\kappa}(\omega) \propto e^{i\delta(\omega)}$ where $\delta(\omega)$ is the phase shift.

Now if in reality $x(t)$ is not periodic and is causal since I only know its values in the past, can I still apply the same kernel to get the instantaneous force? I have been told that I should use Laplace transform instead of Fourier transform. I see the point of it being bilateral by definition, but I am not sure how it is actually different to Fourier transform. Does applying the convolution theorem to a causal signal still give me the desired phase shifts?

What you want is an all-pass filter with frequency response

$$H(\omega)=e^{j\phi(\omega)}\tag{1}$$

where $\phi(\omega)$ is the desired phase shift (and $j$ is how we denote the imaginary unit over here). This system is called an all-pass filter because clearly $|H(\omega)|=1$ holds.

The type of input signal is irrelevant, it can be periodic, non-periodic, causal, or non-causal; if you filter it with a linear time-invariant (LTI) filter with a frequency response given by $(1)$ then the desired phase shift will be achieved.

Your problem is the (causal and stable) realization of such a filter. In general, for a given phase shift $\phi(\omega)$ the frequency response given by $(1)$ cannot be implemented exactly; it can only be approximated.


Filters performance analysis I am working on some experimental data which, at some point, need to be time-integrated and then high-pass filtered (to remove low frequency disturbancies introduced by integration and unwanted DC component).

The aim of my work is not related to filtering, but still I would like to analyze more in detail the filters I am using to give some justification (for example to motivate why I chosed to use a 4th order filter instead of a higher/lower one).

This is the filter I am using:

delta_t = 1.53846e-04;
Fs = 1/delta_t;
cut_F = 8;
Wn = cut_F/(Fs/2);
ftype = 'high';
[b,a] = butter(4,Wn,ftype);
filtered_signal = filtfilt(b,a,signal);

I already had a look here: https://stackoverflow.com/questions/5591278/high-pass-filtering-in-matlab to learn something about filters (I never had a course on signal processing) and I used

fvtool(b,a)

to see the impulse response, step response ecc. of the filter I have used.

The problem is that I do not know how to "read" these plots.

What do I have to look for?

How can I understand if a filter is good or not? (I do not have any specification about filter performances, I just know that the lowest frequency I can admit is 5 Hz)

What features of different filters are useful to be compared to motivate the choice?

For example, what info can I see plotting the round-off noise power spectrum of 2 different filters?(high pass Butterworth of different orders: 4 and 22)

EDIT

I tried decimating the signal to reduce the sampling frequency as suggested in the comments and this is what I get

The green one is the new filter's round off noise spectrum having reduced the sampling frequency (now = 325 Hz), the blue one is the older with sampling frequency = 6.5 kHz.

Things look a bit better but still not good enough, am I right?

OK, doing a design of a Butterworth high pass filter with a cutoff frequency of 8Hz and a sampling frequency of 325Hz using scilab by:

 bw = iir(4,'hp','butt',[8/325 0.5],[0.1,0.1]);

then the filter is:

                                  2
0.8168388 - 3.2673552z + 4.9010328z
               3            4
   - 3.2673552z + 0.8168388z
------------------------------------
                                   2
0.6672256 - 2.9384352z + 4.8677204z
               3   4
   - 3.5960396z + z

and the frequency response looks like:

So I'm not sure what you're doing for the design, but it doesn't appear to be high pass.

One thing that may help is if you could share what your data looks like (say 256 or so points?).


Python FFT (average) Hanning spectrum size 4096 and 50% overlap I want compute a FFT (average) with these parameters:

Spectrum size: 4096
Windows function: Hanning
Overlap: 50%

This is my python code:

data = np.genfromtxt(file_path)
Freq = data[:,0]
Pdyn = data[:,1]
ft = np.fft.rfft(PdynA*np.hanning(len(Pdyn)))
mgft = np.abs(ft)

fig = plt.figure()
xVals = np.abs(np.fft.fftfreq(len(PdynA),d=1/12000))
plt.plot(xVals[:len(mgft)], mgft)
plt.show()

So I have my Hanning windows. But now, I don't know how I can set the Overlap to 50% and the spectrum size to 4096 ?

I can tell you what to do but I can't tell you how to do it in Python code. To average four spectra, do the following:

1) Multiply input samples x[0] –thru- x[4095] by a 4096-point Hanning sequence.

2) Compute a 4096-pt FFT of the above product sequence producing an X1[m] spectrum.

3) Multiply input samples x[2048] –thru- x[6143] by a 4096-point Hanning sequence.

4) Compute a 4096-pt FFT of the above product sequence producing an X2[m] spectrum.

5) Multiply input samples x[4096] –thru- x[8191] by a 4096-point Hanning sequence.

6) Compute a 4096-pt FFT of the above product sequence producing an X3[m] spectrum.

7) Multiply input samples x[6144] –thru- x[10239] by a 4096-point Hanning sequence.

8) Compute a 4096-pt FFT of the above product sequence producing an X4[m] spectrum.

9) Finally, average the magnitude-squared samples of the four spectra to produce a final |Xf[m]| magnitude spectrum.  That is,

|Xf[0]| = (|X1[0]| + |X2[0]| + |X3[0]| + |X4[0]|)/4.

|Xf[1]| = (|X1[1]| + |X2[1]| + |X3[1]| + |X4[1]|)/4.

|Xf[2]| = (|X1[2]| + |X2[2]| + |X3[2]| + |X4[2]|)/4.

...

...

|Xf[4095]| = (|X1[4095]| + |X2[4095]| + |X3[4095]| + |X4[4095]|)/4.

If you desire a power spectrum, then merely square each of the |Xf[m]| magnitude samples. user25139, you're implementing Welch's Method of spectrum analysis. You can search the web for "Welch's Method" for more information.


Beginner BPSK demodulation question I'm writing a simple application in Python that sends and receives data over binary PSK using audio signal. I used cos() for 0 and -cos for 1. On the receiver side I figured I can demodulate by multiplying the signal by 2*cos so this will essentially shift the signal below and above x axis for corresponding bits. It works great. Now, I'm thinking about using low pass filter or maybe averaging the data in the time frames to get the output "square wave"? What would be the best/most efficient approach here? I'm might be a pro software engineer but tinkering with this as a hobbyist, thus probably very basic question :)

An "integrate and dump" filter is a common approach for this application, which is essentially an averaging filter averaged optimally over your symbol period. To implement this, you simply reset your accumulator at the symbol boundaries, accumulate your output over the symbol period and then "dump" your result by taking the accumulated output as the result before resetting it again. For a rectangular pulse shape BPSK in the presence of additive white Guassian noise, this is the optimum matched filter.

Note, unless you have a shared clock between your transmitter and receiver, you will likely also need to implement some form of carrier recovery to track clock variations between your transmitter and receiver (otherwise once you drift 90 degrees, your cosine will be a sine and you will get no output in that condition!). And as described above, you will need some form of timing recovery to establish when your start of the symbols are referenced to your signal as received.

I have detailed these in the following posts:

Phase synchronization in BPSK

Recovering signal for psk

Gardner Timing Recovery for Repeated Sybmols


Frequency at which the maximum amplitude occurs of a signal How can i  find the frequency corresponding to maximum amplitude of each frame after uniform segmentation.

What you want in this case is to find the argmax of the FFT, the index (function 'argument'/arg) which yields the maximum (max), which is obtained via the findpeaks() function in matlab. If you're working in a different language, just search for argmax algorithm. It's not so easy to find the argmax (although the derivative of the signal will tell you a great deal).

The argmax given via findpeaks() will give you the index into the FFT at which the peak amplitude occurs.

Note: Using this technique will only give you a peak amplitude for a frequency $$\pm \, \frac{Fs}{2 \cdot \text{NFFT}}$$

But first off, you need to take the real part of your signal, just look at the second half of it. It's symmetrical.

When you get the argmax, you can compute the frequency of the peak amplitude by taking the argmax and multiplying by the width of your frequency bins, which is $$\frac{Fs}{\text{NFFT}}$$

For example:

If $$argmax(FFT(x[n]) = argmax(X[n])$$

then $$f_{peak} = \frac{argmax(X[n])*Fs}{\text{NFFT}}$$

where $$Fs$$ is the sampling frequency and $$\text{NFFT}$$ is the size of the FFT.


How the zero - phase filter without filtered signal truncation at the end can be implemented? I mean what for the window FIR filter the filtered signal is truncated at the end because impulse response of the filter is symmetric. For example the code from here (in python) (lowpass FIR filter) gives the next results:

color meanings:
blue - input signal
green - shifted filtered signal
red arrow - the end of the filtered signal

The red arrow shows where the filtered signal is ended. How can I implement the zero - phase or nearly zero - phase filter so what the filtered signal would end at the same point in time as the original one? The end of the filtered signal is of most significance.

Thanks.

You can either gather more data at the end than you need for output, or zero pad past the end of your data with at least half the length of your zero-phase filter.

If you can’t get data past the “ends”, and don’t like the artifacts from zero-padding, then you might try other reasonable assumptions (circularity, continuity, etc.) depending on your model.


unexplainable aliases in spectrogram My scope has 200 MHz bandwidth and can sample at 4GS/s

I'm sampling a 138 MHz RF signal at 500 M Sa/s storing 14 M points.

When I plot the spectogram in matlab I get the following
spectrogram(Data,12000,8000,[],fs,'yaxis')

Finally here is the frequency spectrum

While the signal does show up at 138 MHz, there appear to be aliases at 50 MHz, 190 Mhz and complex conjugate aliases at 86 MHz, 172MHz and 224 MHz. Does anybody know why this could happen?

edit:

Looking at the spectrogram, the prominent artifacts go up or down in frequency synchronously to the 138 MHz signal but have larger bandwidths. That is an indication that they are its harmonics, due to a nonlinearity somewhere in the system. With the sampling frequency of 500 MHz some of the positive and negative harmonics alias to the following (bold) frequencies that match with the prominent artifact bands:

-138 MHz * 2 = -276 MHz. -276 MHz + 500 MHz = 224 MHz
-138 MHz * 3 = -414 MHz. -414 MHz + 500 MHz = 86 MHz
138 MHz * 4 = 552 MHz. 552 MHz - 500 MHz = 52 MHz

The apparent bandwidths and complex conjugacy seem to match those calculations.

You did not mention what equipment you are using, but it looks like there is no proper antialias filter.

The less prominent artifact bands can be explained similarly:

138 MHz * 5 = 690 MHz. 690 MHz - 500 MHz = 190 MHz
-138 MHz * 6 = -828 MHz. -828 MHz + 2 * 500 MHz = 172 MHz


Gibbs phenomenon in Hilbert transform I have just read that the Gibbs phenomenon can cause artifacts in the use of Hilbert transforms for audio processing applications, but the reference contained no more specific information.  Can anyone expand on this remark?  When does the Gibbs phenomenon cause problems for the Hilbert transform, and how does one solve this problem?  Specific references to books or articles, and specific examples of signals with difficulties with Gibbs phenomenon for the Hilbert transform, would be appreciated.  Thanks.

Bruce

The Gibbs phenomenon occurs wherever a discontinuous function is approximated by a truncated (Fourier) series. For the design of systems, the Gibbs phenomenon occurs when approximating ideal brick wall filter responses (like low-pass, band-pass, etc.), but also for all-pass systems such as the Hilbert transformer because of its discontinuity at $\omega=0$:

$$H(e^{j\omega})=-j\textrm{ sign}(\omega),\quad |\omega|<\pi\tag{1}$$

The impulse response of the ideal discrete-time Hilbert transformer is

$$h(n)=\left\{\begin{array}{lc}\frac{2\sin^2(\pi n/2)}{\pi n},&n\neq 0\\
0,& n=0\end{array}\right.\tag{2}$$

A naive design approach would be to simply truncate the impulse response given by (2) to obtain an FIR approximation of a Hilbert transformer. The magnitude of the resulting frequency response for an FIR filter with 23 taps is shown by the blue curve in the following figure:

The point '1' on the x-axis is the Nyquist frequency. One can clearly see the Gibbs phenomenon with large ripples near the discontinuity. Note that with an even length filter the approximation would be better around the Nyquist frequency (because such a filter doesn't have a zero at Nyquist), but the Gibbs phenomenon around 0 would remain. The green curve in the figure is the magnitude of the frequency response when the truncated impulse response is windowed by a Hamming window. Clearly, the effect of the Gibbs phenomenon is greatly reduced, but this comes at the expense of a wider transition band.

Using more sophisticated methods for designing FIR Hilbert transformers, such as the Parks McClellan algorithm, weighted least squares, or constrained least squares, the Gibbs phenomenon can be avoided. Of course, there will always remain an approximation error, but it can be made sufficiently small by choosing an appropriate filter length. Note that this is not the case when simply truncating the ideal impulse response. In this case, the maximum ripple size near the discontinuity remains the same, regardless of the filter length. Also note that with FIR filters the phase approximation error can be made exactly zero by choosing an anti-symmetric impulse response. The only error is the error in the magnitude of the response, as shown in above figure. This is different from IIR approximations, which can achieve an exact magnitude of 1, but there will always be a phase approximation error. You can find some basics about the design of Hilbert transformers in the textbook Discrete-time Signal Processing by Oppenheim and Schafer.


SURF feature extractor and descriptor I am trying to understand how SURF extracts features using a Hessian matrix.
What I am a bit confused about, is why a second derivative Gaussian filter is used, rather than a standard Gaussian filter? What effect does this have on the system?

First of all, SURF only uses it as a blob detector to find interest points.
The first order derivatives only give you directional information such as edge-ness. However, the Hessian (as the matrix of second derivatives) contains local structure information such as curvature / concavity etc. For instance, its eigenvectors can be used for cornerness measures. It is just more informative. Also, its determinant is used as an automatic scale selector.

Another common alternative would be a Laplacian detector, but Hessian shows better scale selection properties.


When i have already found beamforming,do i still need to use MRC?i have some problems for them I have some problem when i learning MIMO,first thing is that i know designing beamforming and MRC can ensure better receive signals,but they are not the same,so i have some question about them.

1.I know the MRC is the one of method to ensure better received signal,so if I have already  found the beamforming direction now,do i still do MRC when calculating the SNR?,like the example below

2.If the transmitter use $N_T$ antenna to transmit one signal to the receiver which is with $N_R$ antenna,now  i do SVD to the MIMO channel $\mathbf H$ ($N_R \times N_T $ matrix ),and find the best beamfoming direction $f_A,(N_T \times 1$ column vector), from it.Now i have beamforming now ,so i don't need to do MRC?

3.Now the received signal is $y_t=\sqrt{P}\mathbf H \vec f_Ax+\vec n$,so now the $SNR=\frac{P||H\vec f_A||^2}{\sigma^2 _n}$?

Because i heard my classmate said that $SNR$ is not $\frac{P||H\vec f_A||^2}{\sigma ^2 _n}$,but $\frac{P||\vec u^HH\vec f_A||^2}{\sigma ^2_n}$,and this $\vec u^H$ is calculate the SVD too.The first column of $U$,SVD=$U\Sigma V^H$.Because $H\vec f_Ax$ is still a vetor,and we can't not calculate SNR from a vector,it must be calculate from a value,and $\vec u^HH\vec f_A$ is a value

Is my SNR formula right or my classmate's is right?Does anyone know the answer about my question?It confused me for lots of months

You are doing transmit beamforming and choosing your vector $\mathbf{f}_A$ so that you beamforming in the "best" direction. So you transmit the signal: $\mathbf{z}=\sqrt{P} \mathbf{f}_A x$, where $\mathbf{z}$ is a length $N_T$ vector. Now you receive the signal: $\mathbf{y} = \sqrt{P}\mathbf{H}\mathbf{f}_A x + \mathbf{n}$, which is a length $N_R$ vector.

The post-MRC SNR should still be a vector. You will have a SNR for each transmit stream, so $N_T$ SNR values. To compute the SNR of the $k^{\text{th}}$ transmit stream after MRC you take the channel gains from transmit antenna $k$ to each of the $N_R$ receive antennas, $\mathbf{h}_k$, and compute its output power: $\gamma_k = P ||\mathbf{h}_k^H \mathbf{H} \mathbf{f}_A||^2$. Or, you can also equivalently do it all at once by: $||\mathbf{H}^H \mathbf{H} \mathbf{f}_A||^2$.

Your formula is actually computing the pre-MRC SNR, whereas your friends was showing how to get the post-MRC SNR. Hope this helps!


Why is $x[n]=\sin(\frac{12\pi n}{5})-\sin(\frac{2\pi n}{5})=0$ I was trying to find the period of $\sin(\frac{12\pi n}{5})-\sin(\frac{2\pi n}{5})$, each of the sinusoids has a period of 5 however their difference has a period of 1. It turned out that they're identical in the discrete-time domain.
besides actually solving the equation for reals and concluding integer roots, is there anything that explains this in terms of periodicity?
Thank you.

HINT:
$$\sin\left(\frac{12\pi n}{5}\right)=\sin\left(\frac{2\pi n}{5}+\frac{10\pi n}{5}\right)=\sin\left(\frac{2\pi n}{5}+2\pi n\right)$$
Taking into account that $\sin(x)$ is $2\pi$-periodic should make the result obvious.


Is there a way to double all frequencies? Let's imagine I have an input signal that looks like this:

$$
f(t)=\sin\Big(2\pi(2)t\Big)+\sin\Big(2\pi(3)t\Big)+\sin\Big(2\pi(5)t\Big)
$$

And I want an output signal that scales each frequency proportionally by some factor. Here, I scale everything by a factor of '2':

Is there a way to do this?  The frequency domains look like this:

That's typically referred to as "pitch shifting" and it's typically a combination of re-sampling and time stretch/compression. Re-sampling is straight forward but time stretch/compression is a bit more tricky and the best methods depends on the application requirements, the characteristics of the signal and the amount of stretch/compression you'll need. See http://blogs.zynaptiq.com/bernsee/time-pitch-overview/ and https://en.wikipedia.org/wiki/Audio_time-scale/pitch_modification to get started


Finding Frequency Differences In Differing Amplitude Signals I have two audio signals, both with supposedly matching frequencies at each time window.  What I would like to do is to verify that the frequencies are identical.  However, one was created at a higher power than the other, so simply checking for a match at each sample won't work.  I've already cross-correlated the signals, so there's no time lag between them.

Is there any way to compare the "shapes" of two waveforms of differing amplitudes?  If there's some way to drop the amplitude domain from the signal, that would work - I don't know nearly enough about DSP to know if it's easy or not.

Sounds to me like you need to look at the signals in the frequency domain.  Take the Fourier transform of each signal (use an FFT), and look at the spectrum.  I'm not sure if you've got signals with just a single frequency (a tone) or with many frequencies, but in either case if they have the same frequencies in them then the spectra should look the same (except for one being larger than the other).

If you're unfamiliar with frequency domain concepts, I'd suggest you read up on them first.  The wikipedia page on Fourier transforms is pretty large, so I'm not sure how good that would be as an introduction, but there are a large number of books on the topic and many websites too.


How to demonstrate the orthogonality of OFDM subcarriers? Assume if we have $N$ OFDM subcarriers represented  by results of the inverse FFT of $N$ data symbols $\mathbf x$. As I know, the subcarriers of OFDM should be orthogonal. It means that
$$X(n)X(n+1) = 0\quad\text{where}\quad n = 0,\ldots,N-1$$
My question is that I can't get that right when doing it in MATLAB. For example suppose that $N=4$ and
$$\mathbf x = \begin{bmatrix}
0.7+0.7i\\
0.7-0.7i\\
-0.7+0.7i\\
0.7+0.7i\end{bmatrix}$$
is a data symbols gotten after QAM modulation.  The OFDM $N\times 1$ subcarriers are gotten by iFFT of the data symbol as below:
X = ifft(x)*sqrt(4);

It means that
$$X = \begin{bmatrix}
0.7+ 0.7i\\
1.4 + 0.0i\\
-0.7 + 0.7i\\
0.0 + 0.0i\end{bmatrix}$$
Then,
$$\begin{align}
X(1)\cdot X(2) &= (0.7000 + 0.7000i)\cdot(1.4000 + 0.0000i)\\
& = 0.9800 + 0.9800i
\end{align}$$
If subcarriers are orthogonal, that wil be zero.
Thus, my question, how can I demonstrate the otrhognality of OFDM subcarriers ?

Orthogonality is defined as "the inner product of two vectors equals zero".
Now, in OFDM, the transmit vector for a single subcarrier is exactly one row vector $\mathbf D_k$ of the DFT Matrix $\mathbf D$, multiplied by the complex value of a symbol $c_k$, i.e. $c_K \mathbf D_K$.
Two different subcarriers $k, l, k\ne l$ hence have the inner product $\langle c_k\mathbf D_k,c_l\mathbf D_l\rangle$; inner products are linear things, hence that's
\begin{align}
\langle c_k\mathbf D_k,c_l\mathbf D_l\rangle &= c_kc_l \langle \mathbf D_k,\mathbf D_l\rangle\\
&= c_kc_l \begin{cases}0&k\ne l\\\|\mathbf D_k\| & k = l\end{cases} &\text{q.e.d.},
\end{align}
because the DFT matrix is unitary.


Bandwidth range for Fast Fourier vs principal component analysis? I've read somewhere that the Fast Fourier is only applicable to those processes exhibiting bandwidth. Where as principal component analysis can be applied to a process exhibiting any finite  bandwidth. Why is this?

The band width is simply the difference in the upper and lower frequencies in a contentious set of frequencies? So, for the Fast Fourier only a set of basis functions of discrete frequencies (i.e. $\sin(2\pi f x)$) are considered, but why does the bandwidth ($f_\text{max}-f_\text{min}$) need to be infinitesimal? Why can it not be finite? Is is due to the fact that the duration of the sine wave is infinite and therefore the bandwidth infinitesimal? Or have I miss understood? I cannot find this information anywhere.

From Karhunen–Loève theorem, when talking about stochastic processes:

  In the theory of stochastic processes, the Karhunen–Loève theorem
  (named after Kari Karhunen and Michel Loève), also known as the
  Kosambi–Karhunen–Loève theorem is a representation of a
  stochastic process as an infinite linear combination of orthogonal
  functions, analogous to a Fourier series representation of a function
  on a bounded interval.

Basically, K.-L. yields an adaptive representation, while Fourier provides a fixed representation with  sinusoidal functions.

There are many flavors of Fourier tools, so since you are talking about "fast" Fourier and PCA, I assume for now you are dealing with discrete data, thus the empirical version of the  Karhunen–Loève transform (PCA) and the discrete Fourier transform (DFT). The idea of a continuous set of frequencies does not fit gracefully here.

Both PCA et discrete Fourier can be cast into linear and orthogonal transforms, so they are applicable to any data, and may preserve all information.

From Principal Component Analysis (p. 44 sq.), you can find that you can derive DFT bases from PCA, when you are studying a process that follows a "correlated" Markov model.

Finally, when you talk about infinite sines with infinitesimal bandwidth, you are more in the context of continuous functions (or distributions) and standard Fourier analysis, which might be a cause for confusion. A solid book on the topic could help, like Mathematical principles of signal processing: Fourier and wavelet analysis, by P. Brémaud.


What should be my sampling frequency? I have an accelerometer that I am reading the x, y, z values for m/s^2. I've written up my code in Python as follows.

while True:
  x, y, z = sensor.accelerometer
  with open('data.txt', 'a+') as f:
    f.write('{},{},{}\n'.format(x, y, z))
    f.flush()
    f.close()
  time.sleep(0.1)

The time.sleep(0.1) says to sleep for 1/10 of a second, so theoretically, in 1 second, I should have 10 samples. Thus, I would assume my sampling rate is 10 samples per second. However, what I have noticed is that I do not get 10 samples per second; there is overhead with the file writing (in actuality, I am doing more stuff than just reading and writing the accelerometer data). In fact, in one second, I get a distribution over the number of samples per second.

Now that I have recorded the data, I need to analyze these signals, but with a lot of the signal processing examples, I see that the sampling frequency is a required parameter.

My question is, what should I set as my sampling frequency?

Should I put in 10? If so, empirically, I am not consistently (if at all) getting 10 samples per second.
If I do not use 10, I have a distribution over the number of samples per second, should I just use the mean? If I use the mean number of samples per second as the sampling frequency, do I need to remove samples within a second interval?
What do I do about a second of interval with less than the mean number of samples?

Any help into these questions is appreciated.

Much depends on what you want to deduce from your data.

In general,  if your sample rate is not uniform, your measurements should be accurately time stamped.

The nearer to an average of 1/10 a second your sample intervals are the better.

A typical hueristic used in Engineering is the 1/10 rule, so if your samples are within 1/100 of a second of 1/10 second sampling, you probably can just assume they are uniformly sampled, but again, it really depends on what your measurement's purpose is.

Again if you have timestamps, there are methods to interpolate your measurements to uniform time increments. There is the timeseries class in Matlab which has a method for this purpose, which is an example of how to treat nonuniform samples.

If you don't have the ability to accurately time stamp data, you could use the average of sample intervals.  Like all practical solutions to problems you need to use your own judgement to determine if what you do is adequate.

We really can't make that judgement for you.


Is there a rule of thumb for selecting the variance of the input? Considering an estimation problem, where I want to estimate the unknown input $x$ using the known channel parameters. This estimation problem can be solved by Least Squares.

Thus, the model is $$\hat{x}[n] = \mathbf{A}^T\mathbf{x}_{n} + w[n]$$ where $\mathbf{A} = [a_0,a_1,\ldots,a_{L-1}]^T \in \mathbb{R}^{1 \times L}$ represents the channel's impulse response, and $\mathbf{x}_{n} = [x[n], x[n-1], x[n-2],\ldots,x[n-L+1]^T$  and $L$ is the order.
The model is an FIR filter and $w[n]$ is the measurement noise which is $w \sim N(0,\sigma^2_w)$.

After the estimation, I want to see the performance using Mean square error defined as $MSE = \frac{1}{NT}\sum_{n=1}^Ne[n]^2$ where $T$ is the number of independent runs, that is for every SNR, I generate $T = 10$ different channel parameters at random and for each  $T$ I perform estimation.

The error is $e[n] = x[n] - \hat{x}[n]$.

The graph to plot would be SNR on $X$ axis vs. $MSE$ on $Y$ axis. SNR is defined as $\frac{\sigma^2_x}{\sigma^2_w}$.

But I don't know what is the value of the numerator, $\sigma^2_x$. I cannot understand how to plot the curve. Can somebody please explain (with an illustration / graph would be very useful) what is the procedure to plot the graph, what is the value of variance of the input, $\sigma^2_x$?

Please correct me where ever I have done mistake (especially, in the formula ).

You can define your model as $$\hat{x}[n] = \frac{1}{||\mathbf{x}_{n}||}\mathbf{A}^T\mathbf{x}_{n} +\frac{1}{\sqrt{\mathsf{SNR}}} w[n]$$
In this way, you define a range for $\mathsf{SNR}$ in dB (e.g. between 0 to 10 dB), then translate it to natural (non-dB) scale, and then $\frac{1}{\sqrt{\mathsf{SNR}}}$ gives you the noise power. The signal power is normalized to one which makes it easy.

Here, $w\sim \mathcal{N}(0,1)$ is a standard Gaussian RV and you can simulate it as w = randn;. In fact, $\frac{1}{\sqrt{\mathsf{SNR}}}w=w'\sim \mathcal{N}(0,\sigma^2_{w'})$, where $\sigma^2_{w'}$ has the right proportion wrt the SNR.


Impulse Response for an Anti Causal Linear Shift Invariant System I have a system given by
$$y[n] - \frac{1}{4} y[n-1] - \frac{1}{8} y[n-2] =3x[n] $$
I'm asked to determine the impulse response of both an anti-causal and a causal Linear Shift Invariant System characterized by this difference equation.
The impulse response to the causal system is $$2\left(\frac{1}{2}\right)^{n}u[n]+\left(-\frac{1}{4}\right)^{n}u[n]$$
I found this without using Z transform or Fourier Transform but by solving the difference equation using the initial condition $h[n]=0$ for $n<0$.
How can I find the anti-causal impulse response without using Z-transform or Fourier Transform.
What boundary conditions do I use now?

You've already found the solutions of the characteristic equation as $\lambda_1=\frac12$ and $\lambda_2=-\frac14$, so you know that the solution must have the form
$$y[n]=c_1\left(\frac12\right)^n+c_2\left(-\frac14\right)^n\tag{1}$$
For the anti-causal solution we know that $y[n]=0$ for  $n>0$. The constants $c_1$ and $c_2$ are simply determined by making sure that the difference equation is satisfied for $x[n]=\delta[n]$, given that we require $y[n]=0$ for $n>0$:
$$\begin{align}n=2:\quad&&&-\frac18y[0]&=0\\n=1:\quad&&-\frac14y[0]&-\frac18y[-1]&=0\\n=0:\quad&y[0]&-\frac14 y[-1]&-\frac18y[-2]&=3\end{align}\tag{2}$$
From the first equation in $(2)$ we see that the impulse response actually starts at $n=-1$, not at $n=0$. The constants $c_1$ and $c_2$ can be determined from the second and third equation, which simplify to
$$y[-1]=0\quad\textrm{and}\quad -\frac18 y[-2]=3\tag{3}$$
Plugging $(3)$ into $(1)$ finally gives the actual values of $c_1$ and $c_2$, which I'm sure you can figure out yourself. The final anti-causal solution has the form
$$y[n]=\left[c_1\left(\frac12\right)^n+c_2\left(-\frac14\right)^n\right]u[-n-1]\tag{4}$$


Real time background noise reduction or removal I am trying to build a real time simple audio recognition that can recognize few keywords. Following this tensor flow tutorial i was able to train a simple neural network.
Currently the signal flow looks like this

  microphone --> Pulse Code Modulation(PCM - range [-1, 1])  --> stfts
  --> spectrograms--> mel_spectrograms --> log_mel_spectrograms --> mfcc's --> Neural Network.

So far, I just grab the data from the microphone input and preform some steps and input to Neural network. Not performing any kind of back ground noise reduction or some kind of signal filtering techniques. So far the result are very much influenced by the back-ground noises( single channel sampling frequency 16000 Hz).

I am new to this subject and wondering if you can share some valuable insight and techniques/python packages to enrich the speaker voice before feeding it to neural network.

A typical first approach to solve your problem would probably be to use one of the many existing noise reduction algorithms that target speech signal applications. These algorithms typically consist in

noise PSD estimation and
noise reduction filtering.

For each task, there are a variety of algorithms to choose from. Example noise PSD estimation algorithms are those based on minimum statistics of the input signal [1] or based on the minimum mean-square error [2]. The actual noise reduction filtering can, e.g., be performed by using spectral subtraction or a Wiener filter.

References

[1] Rainer Martin, “Noise Power Spectral Density Estimation Based on Optimal Smoothing and Minimum Statistics,” IEEE Transactions on Speech and Audio Processing, vol. 9, no. 5, pp. 504–512, Jul. 2001. PDF MATLAB implementation

[2] T. Gerkmann and R. C. Hendriks, “Unbiased MMSE-Based Noise Power Estimation With Low Complexity and Low Tracking Delay,” Audio, Speech, and Language Processing, IEEE Transactions on, vol. 20, no. 4, pp. 1383–1393, May 2012. PDF MATLAB implementation


What's the advantage/disadvantage of oversampling followed by decimation, verses sampling at the correct rate to begin with? I read somewhere that oversampling somehow reduces sampling noise?  How to quantify the reduction in noise?  and why would it reduce noise if you are decimating the oversampled signal.  What's the advantage of oversampling followed by decimation, verses just sampling the signal at the decimation rate to begin with?  Or is the advantage moot, and its generally better to just sample at the decimation rate to begin with and forgo the oversampling.  would there be a case where one is better than the other, and vica versa?

I may have originally misinterpreted your question so I'll change my answer here.

As to why we don't typically just sample at our a decimation rate (i.e. the Nyquist rate for the bandwidth of interest). It really comes down to the practicality:

If you had a 200hz signal of interest, yes you could sample at 400hz assuming there was absolute no noise in your system. However, there will be noise and you would need an anti-alias filter to block all frequencies greater than 200hz. The issue with this is that many anti-aliasing filters are slow roll-off because they're cheaper and easier to implement. Higher order analog filters are simply more expensive and difficult to implement (they're more sensitive to component tolerances).

Instead, it would make more sense for us to over-sample our signal such that we could use a lower order anti-alias filter.

Let's say your signal of interest is still 200hz and you decide to sample at 10khz this time. You put in an anti-alias filter at say 5kHz with a slow roll off. This would prevent aliasing frequencies of getting into your digital signal while being able to use a cheap/easy anti-alias filter. You would then perform your decimation operation. Note that the decimation operation includes a digital low-pass filter such that aliasing doesn't occur in the down-sampling process.

Summary: the further we separate our bandwidth of interest from the Nyquist rate, we can use lower-order anti-alias filters which are cheaper to implement.


Does this lowpass filter have a noncausal impulse response? I'm doing some measurements with my sound card, I have one output channel connected to one input channel.

I then send the sound card a unit impulse, i.e. a signal 1 value followed by zeros.

The recorded response is shown below.

I know that the Fourier transform of a Rect function, or LPF, will result in a Sinc function - which is what the response below seems to be.  But I'm not quite sure as to the reason for the pre-ringing in the response.

Somewhere in the chain between the sound card analog input and the samples that you're plotting, there is certainly a lowpass filter. There is likely an analog anti-aliasing filter before the ADC; in addition, there are likely one or more lowpass filters applied during resampling processes either on the card or in the operating system's audio driver stack (resulting in a stream of samples at the sample rate that you requested).

As you noted, lowpass filters often have an impulse response that resembles a $\text{sinc}$ function. When you excite a linear system with an impulse, the resulting output is just a copy of the system's impulse response (due to the convolution theorem). So, the fact that the observed signal looks like a $\text{sinc}$ function is to be expected. It doesn't sound like this surprised you, either.

However: you did not observe noncausal behavior. If you had, then you should have run immediately to your local patent office. Noncausal systems are not realizable in the real world. Remember the definition: in order for a system's response to be noncausal, its output must lead the input in time. Stated differently, the filter would start outputting its response to the input before you put the input in. Obviously, that's not going to happen.

So what's with the pre-ringing? Simple: the large spike in the middle corresponds to the center of the filter's impulse response. However, that spike does not correspond to $t=0$, the time at which the impulse was inserted into the filter. Consider the following MATLAB example:

% generate a 250th order lowpass filter
b = fir1(250, 0.5);
% plot its impulse response
plot(0:250, b); grid on;

The resulting plot looks like this:

As you can see, the peak in the plot is not at sample index zero. The filter has an overall delay of 125 samples (indeed, all linear-phase FIR filters have a bulk delay of $\frac{N}{2}$ samples, where $N$ is the filter order), and the pre- and post-ringing are centered about that delay. So, when you insert the impulse, you initially see very little response. As the impulse makes its way through the filter taps, the ringing ramps up to a peak at the impulse response's center, then ramps back down to zero.

The takeaway: There is no noncausal behavior in the example you gave. It is possible to simulate noncausal filters in practice by adding enough delay, similar to the delay shown in the lowpass filter's impulse response above.


How to remove background noise from sound file for analysis I'm working on a project at my university and it involves background noise recognition/classification.  I'm very new to this so I'm unsure what to really search/google or where to really begin any research.

I'd like to remove background noise from any audio file and analyze THAT and not the primary signal (for example, if you had a recording of someone talking on a phone, remove their voice and keep the remaining background noise).  I then wish to train a classifier to cluster or group possible sound samples.

How does one do this (extraction and preseervation of background noise)?  I'm guessing I need to do a FFT on the raw data?  Are there proper file formats one needs to have an audio file in to do this?  Additionally, what are proper learning or classification techniques for audio recognition?

So you basically have a noisy speech signal and want to analyze only the noise without speech? This seems like a task for a voice activity detector (VAD). Assuming it works properly it should give you the parts of the signal containing speech and the parts of the signal containing noise. This approach is often used in speech enhancement where an estimate of the noise is required in order to suppress it from the noisy speech signal (usually the noise power spectral density or PSD is calculated). This seems to be pretty close to what you want to do apart from the speech enhancement part which you don't need. This is the way to go instead of "removing the voice" as this is more difficult. If you have access to the Digital Speech transmission book by Peter Vary, it has a nice section on VADs. You can also look at http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html - very useful toolbox for any type of audio processing.

So let's say you have the parts of the signal containing only noise. So you then want to classify it in how many classes? The most important thing for the classification is your choice of classification features - basically, how will the different noise types differ. You can have external/ambient noise as a property of the environment in which the communication is taking place but also internal noise as a property of the communication system itself. Keep in mind that competing speakers can also be considered noise, based on the application. In that case the VAD will have more problems.

So, you need to define your features and then you can extract them. Yes, you can do the FFT on raw data and calculate the PSD for example and look at the energy in certain frequency bands. This is a nice (and free) book on spectral analysis you can use as a reference: http://user.it.uu.se/~ps/SAS-new.pdf.

Anyhow, try to get started on it (people have already done similar things -google noise classification) and do it step by step. So there's no use worrying about classification when you don't know what features you'll use to distinguish between the signals.

Hope this helps


Why does QAM use two cycles instead of one per byte? If you look at this image, you can see that you have 256 possible binary values by modulating the amplitude (the radius of the trigonometric circle) and the phase (the degrees) of the carrier.

However, multiple diagrams like this one show that each value is actually two cycles long.

My question is why is there two cycles per value instead of one cycle? The constellation only goes to 360 degrees, so why is two cycles used? Wouldn't you be able to transfer data 2x faster if you used one cycle per value?

It doesn't matter how long you stay on one constellation point. The constellation diagram has NO ties to time whatsoever. The two base vectors $e_x, e_y$ (x-axis and y-axis) represent your real and imaginary part or "basic functions" your signal is made of. For example, if your signal consists of two functions (in this case, QAM, it consists of a real and an imaginary part including the amplitude)

$$
e_x = A\cos(2\pi f_{carrier}t) \\
e_y = A\sin(2\pi f_{carrier}t) \\
$$

If you stay on one constellation point for an hour, that's fine. The constellation diagram only shows you which "constellation" of amplitude and phase you should choose to transmit a specific "constellation" of bits (001 010 110, etc ..).

The question you ask is about the symbol rate. If the symbol rate is half your carrier frequency, then you will always have two cycles before the next symbol will be sent. If your symbol frequency is a quarter of your carrier frequency, you will have four cycles to the next symbol.

I think there are always at least one and a half cycle shown to make it easier on the eye and for the brain to understand because we can follow the sine easily. This wouldn't be the case if you'd change the phase and amplitude every cycle. That being said, there's also the problem of having enough time for the demodulator to detect the phase shifts. This is why the symbol frequency is often a lot lower than the carrier frequency.


Emotional Tone detection using audio on mobile phones Which audio features should be used to detect tone from audio? Are there any open source frameworks or implementations which can be ported to mobile phones to process audio?
Alternatively, would it be possible to collect audio features from the phone and then process this on the server so as to detect tone?

I dont know how bound your problem is... Is this dealing only with words that are already detected?

This is a very very hard problem, as @JimClay mentions. For such problems, one cannot simply say "pick highest "x"", because the features are not obvious to us. (Even though they might be obvious to non-observable states within our brains, whose signal processing paradigms are unfortunately inaccessible to us via introspection).

So in problems such as those, I would first collect a large sample of "angry" tones, and a large sample of "happy" tones to contrast. You now have two options:

1) I would 'brute force' it and start looking for features to extract that separate them amicably in feature space. It is honestly hard to say where to begin without looking at real data. VERY coarsely speaking, I imagine angry tones have more low-frequency information, VS higher-pitched happy ones. This then lends itself to including the DFT of the data as 1 feature in DFT_length dimensional space. Of course, not all 'low pitchy' sounds are angry, and not all high-pitchy sounds are happy, so you need another feature vector(s) on which to discriminate, lest you inadvertently create a male-female detector...

2) I would use a different approach and simply label my data vectors, (here, each data vector is the time domain snippet of detected words, all of which are re-sampled so that they have the same pre-determined length). I would then start off with a simple linear perceptron to try and create a hyper-plane to best separate the two classes. (You can add an extra dimension to your classifier to make linearly un-separable clusters separable by taking your hyper-plane up one more dimension). Once adequately trained, you can use that weight vector as your classifier and project all data vectors unto it to determine which class it falls under. Bear in mind you need A LOT of training data for this to work properly, because there are many situations where emotions can manifest itself in words and phrases. This is why I asked you how bound your problem is. (Words? Phrases? Sentences?, etc).


Difference between a digital lock-in amplifier and a FFT when extracting phase of a signal? I am trying to measure the relative phase of a sine wave fixed at a particular frequency in a noisy environment. My initial approach is to simply collect $N$ samples, take an FFT, and then just extract the phase at the operating frequency (which I know apriori).

When reading about methods to extract signals from noisy environments, I came across the lock-in amplifier. However, I am confused about whether I should have any reason to expect it to perform better than the FFT (ignoring run time, just looking at ability to extract my phase).

Specifically, imagine that I implement the lock-in amplifier digitally. I would then numerically calculate the following two integrals (source):

\begin{align}
X_{LI} &= \frac{1}{T}\int\limits^{T}_{0} U_{\rm in}(s) \cos(2\pi f s) \,ds\\
Y_{LI} &= \frac{1}{T}\int\limits^{T}_{0} U_{\rm in}(s) \sin(2\pi f s) \,ds
\end{align}

and then I would get my phase using:

$$
\theta_{LI} = \tan^{-1}\left(\frac{Y_{LI}}{X_{LI}}\right)
$$

However, if I were to use an FFT, the FFT would evaluate

\begin{align}
X_{FFT}(f) + iY_{FFT}(f) &= \int\limits_0^T U_{\rm in}(s)e^{-2\pi i s f}\,ds\\
&=\int\limits_0^T U_{\rm in}(s)\cos(2\pi f s)\,ds + i \int\limits_0^T U_{\rm in}(s)\sin(2\pi f s)\,ds
\end{align}

Which seems to imply that just as with the lock-in amplifier, I would get

\begin{align}
X_{FFT} &= \int^{T}_{0} U_{\rm in}(s) \cos(2\pi f s)\,ds = \frac{1}{T}X_{LI}\\
Y_{FFT} &= \int^{T}_{0} U_{\rm in}(s) \sin(2\pi f s)\,ds = \frac{1}{T}Y_{LI}\\
\theta_{FFT} &= \tan^{-1}\left(\frac{Y_{FFT}}{X_{FFT}}\right) = \tan^{-1}\left(\frac{Y_{LI}}{X_{LI}}\right) = \theta_{LI}
\end{align}

So does this mean that I would always get the same result whether I use a digital lock-in amplifier, or a digital FFT  on the same data set?
What is the benefit of using a lock-in amplifier?
Is there no advantage to the lock-in amp in my application?
What is then an example of an application that the lock-in amplifier is best suited for?

Is there a difference between the FFT and a lock-in amp?

Yes, two of them:

The FFT assumes the signal at its input to be periodic.

What this means about your FFT integrals is that they are missing a phase variable (let's call it $\phi_r$) which will be random because it depends on the ratio of your signals period (or its frequency), the window of the FFT and the initial phase that the FFT will "catch" the signal on the first frame (and any subsequent phase changes whether desired or not). By the way, Frequency is integrated Phase.
Consequently, your $\theta_{LI}$ will depend on that $\phi_r$ (which is all over the place).
Intuitively now, the FFT assumes the signal to be periodic. That is, it repeats in the same way both to the left and to the right of the observation window. This is "alright" if the window happens to be an integer multiple of the signal's period but this is extremely unlikely because of noise and possible changes in the phase of the signal. It will therefore sound like a whistle with regular "pops" because of the discontinuities and at each "pop" the phase estimation will be disturbed. There are workarounds to this (and this one) and this brings us to the second difference.

There is a concept in the FFT called Spectral Resolution.

Spectral resolution relates the physical frequencies (in Hz) with each of the distinct harmonics that the discrete FT evaluates its integrals at. Therefore, you would have to accurately calculate the length of the window of the FFT with respect to the sampling frequency and the frequency of the input signal so that it lands exactly on one of the bins. And after doing this of course, it would seem a waste to be evaluating the FT for all harmonics when all you are interested in is just one of them.

Having said this, let's try to tackle the rest of the questions:

  So does this mean that I would always get the same result whether I use a digital lock-in amplifier, or a digital FFT on the same data set?

Provided that the above details are taken care of, the result will be almost identical. I say almost because there is another consequence from #1 above that is not exactly obvious. Even if you were to use overlap-add or overlap-save, you would still have no control over the initial phase at which the FFT would "catch" the incoming signal. And once the FFT "starts", it's own local oscillators are not going to adapt to the incoming signal. Therefore, both techniques would provide information about phase but the FFT's estimate would suffer by a systematic error proportional to the initial phase of the incoming signal. The lock-in amp takes care of this with its Phase-Locked-Loop (PLL). For (much) more information, please see this link.

  What is the benefit of using a lock-in amplifier?

Should be becoming clearer by now. A lock-in amp can return amplitude, power and phase information at one specific frequency with the accuracy of its estimate depending on the clarity of the reference signal (the local oscillator).

  Is there no advantage to the lock-in amp in my application?

On the contrary, there is a massive advantage. And this advantage is coming from the use of the PLL. The PLL is an automatic control system that adjusts the frequency (and therefore phase) of a local oscillator to match the frequency and phase of the incoming signal. It therefore can "adjust its frame of reference" and it returns much more accurate phase information.

  What is then an example of an application that the lock-in amplifier is best suited for?

Please see above PDF.

Hope this helps.


Deconvolution - Richardson Lucy vs. Wiener Filter I am studying some deconvolution techniques, In order to remove motion blur, like:

Richardson-Lucy
Wiener

Are there any pros / cons of using one versus another?
For example which are the pros / cons of Richardson-Lucy technique?

Both are the MMSE estimators.

The main difference is Wiener is the optimal for Gaussian Noise while Richardson Lucy assumes Poisson Noise.

Poisson Noise is a better model for noise in photos captured by a Photo Diode.

Computationally, in the case of Gaussian Noise and Linear Convolution the solution has a closed form solution in the Maximum Likelihood / LS sense.

The Lucy Richardson method, which is a Maximum Likelihood in the Poisson case, has no closed form solution and requires iterative approach.

See Noise, Image Reconstruction with Noise (EE367/CS448I: Computational Imaging and Display, Class 10, Gordon Wetzstein, Stanford University).


Find autocorrelation of exponential signal $a^nu[n]$ I need to find the autocorrelation of the following discrete signal $$x[n]=a^nu[n] $$ So I tried finding the convolution of $x[n]$ and $x[-n]$.
\begin{align}
\phi_{xx}[n]&=\sum_{m=-\infty}^\infty x[m]x[m-n]\\
&=\sum_{m=-\infty}^\infty a^m u[m]a^{m-n}u[m-n]\\
&=\sum_{m=n}^\infty a^m a^{m-n}\\
&=a^{-n}\sum_{m=n}^\infty  a^{2m}\\
&=a^{-n}\sum_{l=0}^\infty  a^{2(l+n)}\quad\tag{with $m-n=l$}\\
&=a^{-n}a^{2n}\sum_{l=0}^\infty  a^{2l}\\
&=a^{n}\sum_{l=0}^\infty  \left(a^2\right)^l\\
&=a^n\frac{1}{1-a^2}
\end{align}

However the result should be an even function of $n$ and instead of $$\frac{a^n}{1-a^2}$$ I should have found $$\frac{a^{\lvert n\rvert}}{1-a^2}$$

Let $x[n]=a^nu[n], |a|<1$. Autocorrelation is

$$\phi_{xx}[n]=\sum_{m=-\infty}^{\infty}x[m]x[m-n]=\sum_{m=-\infty}^{\infty}a^mu[m]a^{m-n}u[m-n]$$
First assume that $n>0$. In this case, we have

$$u[m]u[m-n]=\begin{cases}0,& \forall m<n\\
1,& \forall m\ge n\end{cases}$$
Therefore,
$$\begin{align}
\phi_{xx}[n]&=\sum_{\color{red}{m=n}}^{\infty}x[m]x[m-n]\\
&=\sum_{m=n}^{\infty}a^ma^{m-n}\\
&=a^n(1 + a^2 + a^4 +\cdots )\\
&=\frac{a^n}{1-a^2}\end{align}$$

For $n<0$:

$$u[m]u[m-n]=\begin{cases}0,& \forall m<0\\
1,& \forall m\ge 0\end{cases}$$

$$\begin{align}
\phi_{xx}[n]&=\sum_{\color{red}{m=0}}^{\infty}x[m]x[m-n]\\
&=\sum_{m=0}^{\infty}a^ma^{m-n}\\
&=a^{-n}(1 + a^2 + a^4 +\cdots )\\
&=\frac{a^{-n}}{1-a^2}\end{align}$$

and since $\phi_{xx}[n]=\phi_{xx}[-n]$, we can write it for all $n$ as
$$\phi_{xx}[n]=\frac{a^{|n|}}{1-a^2},\ |a|<1$$


How to determine which Fouriers Series terms to use to approximate a signal? I have a signal (a time-series of air temperature values) that I can approximate quite well with a Fourier series. However, the number of terms in the series grows rapidly, to the point that 30 - 40 terms are needed for a good fitting.

So I was trying to understand how can I select only the Fourier series coefficients that convey the most information about the signal, and discard the others.

I could simply choose values of the series parameters a and b above a certain threshold value, but I do think this is not the right way to go.

I need to limit the number of terms because this approximate function will later be used for further elaborations, and the function must be input manually in the code (suggestions on how this task can be automated would be appreciated too).

I am not trying to predict future temperature or anything like that, I just need a function that can reproduce temperature oscillations with minimum period of ~12 hours.

Your inclination to choose the largest magnitude coeffecients is the correct one.  There is a Linear Algebra explanation for this having to do with orthogonaility.

However, your second question is more pertinent, I think.  If your sole motivation for reducing the coefficient count is to reduce data entry effort then you are much better off formatting the output of your DFT analysis as syntax for inputting into your code.  Then a simple copy and paste will do the job for you.  You also greatly reduce your exposure to data entry errors.  Depending on the volume, a custom programming solution which automates the copy and paste process might be economical.

Hope this helps.

Ced


Estimating effect of filter on headroom I am working with beamforming algorithms for speaker arrays. One of the problems we face is that we need to estimate the maximum amplitude of our test signals after we applied FIR/IIR filters to them. I am looking for ways to express the reduction in headroom that a given set of filters would give us. Ultimately I want to know by how much I have to attenuate my input signal to definitely avoid clipping in the DAC, with respect to the filter.

For sine sweeps, this is easy: The maximum of the filters magnitude response is also the maximum amplitude gain we can expect. However, for broadband signals, this is not as simple.

I considered adding up $N$ sine waves of frequencies $\omega_n$ covering the range of 0Hz - Nyquist. If my filter input signal consisted of sines with amplitude 1, the maximum amplitude of this "stacked sine" signal would be $y_{max} = N$.

Looking at an output signal of the filter, the same method can be used to estimate the maximum amplitude: This time, the amplitudes would not be 1 but $|H(\omega_n)|$ and the maximum output amplitude would be $\sum_{n=0}^{N-1}|H(\omega_n)|$.

I could now compare the maximum input amplitude with the maximum output amplitude:

$AmplitudeGain = \frac{\sum_{n=0}^{N-1}|H(\omega_n)|}{N}$

My hope was that I may see by what amount I'd have to attenuate my broadband input signal to avoid clipping. The truth is that I'm really calculating the average of my magnitude response and that there are signals for which the calculated gain would produce clipping. Consider the following magnitude response, which is the result of a superdirective beamforming optimization:

As you can see, the maximum gain is about 13dB, but my formula would return the average which is in this case -12.8dB. This tells me that there is an error in my thinking. Adding up the respective maximum values of all sine-components of a broadband signal should give me the highest amplitude this signal could ever reach, right? Apparently my approach in itself is wrong.

How could I better estimate the maximum amplitude gain of my filter for broadband signals?

For an input with limited amplitude, the maximum amplitude gain of a filter can be calculated as the absolute sum of the impulse response
Scaling with the inverse of the absolute sum, will guarantee that your max output amplitude will never exceed the max input amplitude
In practice, this does NOT work well: It's a very conservative scaling that often results in bad signal to noise ratio and the input signal that would actually trigger the max amplitude gain is extremely unlikely to occur. Scaling this way leaves a lot of headroom on the table.
A better approach is to look at your worst case sine wave gain. Sine waves (or similar signal) do happen in music and audio a lot, so that's a much more likely scenario. Especially in your case, where the max is around 200 Hz.
That still may not result in acceptable signal to noise ratio: in this case you can "eye ball" the gain at the most important frequencies and use a limiter or soft-clipper to manage the occasional overshoot.

In your particular case, I'd probably go with -10dB scaling with a limiter (provided, it's implemented in floating point). At the end of the day this is a trade off between managing clipping/overshoots, efficiency, signal to noise ratio and max output capability. It needs to be optimized within the specific requirements of your application


neural network for 3D pictures I'm wondering if they are some state-of-art algorithms to apply the convolutional neural network approach to 3D pictures, eg. input is no longer a grid of pixels, but voxels. My objective is to extract automatically features from 3D pictures and, cherry on the cake, perform supervised learning using 3D pictures as input data. Maybe others approach exist to extract features, but I would like to stay in the scope of machine learning approach (eg. don't want to use methods such as Fourier transform). Thanks in advance for your insights!

3D Shape Nets is published this year at CVPR (2015). You could access it here. This is a good way to use voxels as input to DNNs. There are also works applying DNNs to medical imaging such as this one and this one. I guess these would help you a lot.


How to tune PID parameters to track a family of reference signals? I have a system that is controlled by a PID controller. I tune the controller to track a reference signal that is a unit step function.
Is there a way to extend the obtained results to be able to track any step function with arbitrary start end end value that occurs at an arbitrary time instant?

As already mentioned, a linear system that tracks a step reference, will track it for arbitrary values and time of application. Due to the internal model principle, in a feedback control system, if the denominator of the product
$$
C(s) P(s)
$$
contains a pole at $s = 0$, then it will robustly track a step function.
The control law
$$
C(s) = \frac{k_{d} s^{2} + k_{p} s + k_{i}}{s}
$$
has a pole at $s = 0$. If the plant  $P(s)$  does not have a zero at $s = 0$, then your system will track an arbitrary step function $R(s) = \frac{A}{s}$, if it is closed-loop stable.

Since the complementary transfer-function, from reference to output, is $Y(s)/R(s) = T(s)$, and
$$
T(s) = \frac{C(s) P(s)}{1 + C(s) P(s)} = \frac{\frac{1}{s} \bar{C}(s) P(s)}{1 + \frac{1}{s} \bar{C}(s) P(s)} = \frac{\bar{C}(s) P(s)}{s + \bar{C}(s) P(s)} ,
$$
this can be seen by using the final value theorem:
$$
\lim_{t \rightarrow \infty} y(t) = \lim_{s \rightarrow 0} s \frac{A}{s}\frac{\bar{C}(s) P(s)}{s + \bar{C}(s) P(s)} = A
$$