bandpass waveform detector? I am hoping that somebody might have a good idea for detecting low-frequency pulses in some rather noisy data (I'm working on improving the snr). I have repurposed a "Raspberry Shake 1D" to record the vibrations made by an adult filaria just below the skin and using a lowpass filter I was able to collect the following as a proof of concept to see if this was even possible. It seems to work.

https://raspberryshake.org/products/raspberry-shake-1d/

The issue here is I don't get any feedback from the device when it is attached and I need to download the stored data later as a "miniseed" file. The device is, however, able to send UDP packets in real-time and I would like to attach a display for a live picture of the noise filtered data. If possible I would like to add a pulse detection algorithm of some kind to make it more obvious when a distinct pulse of the appropriate frequency range, attack, and decay of this pulse is within some basic heuristic bounds.

I could take a deep learning approach but for that, I would need lots of tagged data that has not been collected yet. I am hoping for some form of fft like algorithm just to give me enough real-time feedback while placing the device to better collect my data.

Any thoughts?

Edit1:

What is being detected (shown in the graph) is a small subsonic vibration where the data leading into that waveform is the background ambient noise level, then the oscillations ramp up and back down, then the background level noise again as the graph ends. I think the frequency here is about 4Hz but that will likely vary, as well as the amplitude, so from the detection level, I may need something like what a PLL lock indicator circuit would do if this were a hardware problem, but then this is not a nice sine wave either. I have not done anything like this in software so I am not familiar with what software techniques might even be possible.

What I am envisioning is a python program receiving UDP packets, passing that data through a low-pass filter, then through some kind of FFT processing to find something near the center (4Hz?) frequency range, and the output from that to produce both my real-time graph and some indicator (noise or color change on the display) that the waveform's desired frequency and amplitude are present. It would need to ignore a signal over some specific amplitude, but I'll cross that bridge when I get something that mostly works. I have done this with the static recorded data, and visually identified the event, but I have no clue how to do this from live data in realtime.

btw - The note on "filaria" was only for context, so you could see the problem you are helping to solve, but it is not really important for this discussion. But since someone asked, it's an investigation of parasitic worm one acquires through being bitten by an infected insect, for which there is no diagnostic test in humans for the disease, no disease-specific symptoms, and no known cure for the disease. I'm trying to research a total enigma and these very minute vibrations are just one general property of the disease being investigated.

Very interesting problem.

Your pulse locations can be found by using a sliding window RMS measurement.  Once found, I would recommend centering a DFT frame on the pulse and applying a VonHann window (broadens narrow peaks and narrows broad peaks), then take the FFT.  Looking at the data, you should get a pretty distinct and clear profile.

Your backgound noise may contain some small peaks too, so you may want to analyze a few "quiet periods" to determine what they are so you can discount them in your evaluation.

Precise frequency determination shouldn't be necessary (or even possible) so a simple estimator, like Jacobsen's or MacLeod's, should suffice.

Like wise, a rough amplitude estimate should also be adequate for each peak.  (See my answer here for phase and amplitude estimation: FFT Phase interpretation of input signal with non-integer number of cycles in FFT window)  Don't use the window if you are going to calculate these.

If you know the frequencies ahead of time (or want to do it twice), selecting a DFT bin that has your frequencies near bins will improve these estimation methods considerably.


How to Quantify a Difference in Spectral Information of two Signals? As the title states, I am looking for established methods to quantify the difference in two FFTs. I have come across methods such as using cross correlation in the time domain or the coherence estimate of the power spectral density of two signals, but these methods seem more geared toward finding similarities in signals rather than differences.
I am running an experiment with a treatment and control group and I think the treatment is causing a change in some periodic motion which I am measuring optically and ultimately I want to quantify that change into a single number or metric.
I have thought of perhaps normalizing the FFTs of the two signals I am comparing and finding what frequency value is 50% of the area under the curve for the FFT. The idea being that if the peaks of the FFTs in my treatment group shift then so should the point of 50% of the area under the curve. But I'm not sure if this is a sensible approach.
Thanks for any assistance!

Quantifying a difference in the frequency domain can be useful if the information on how your signals differ is better expressed in the Fourier space. If this assertion is valid in your case, hereafter are hints on methodology, and basically the sort of you can find in Shazam to distinguish music from a short listening period.
By the way, I think that I have been trying to find an answer to a related question  in 50% of my research hours pertaining to spectral or harmonic  analysis (the rest was coding). I am not fully successful yet.  I would however rephrase it as:

How could I easily find quantitative or selective (Fourier) features for
relatively stationary signals?

By quantitative I mean that the feature set should be able to measure scalar values related (in a monotonic way) to a quantity I seek in the original signal. By separative I mean that features that could tell whether signals are "apart". Those two concepts could be related by proper distances or metrics between the feature sets of two signals.
In your proposal, there is a notion of location in the frequency axis, a sort of median of spectrum amplitudes : 50% of the area under the curve for the FFT. A notion of amplitude is skipped, as you may normalize spectra. There is no notion of uncertainty or spread either. I am often trying to consider, as first, feature  triplets of Location-Amplitude-Spread ($f^k_0,f^k_1,f^k_2$) to describe the morphology of signals, spectra, or other suitable sparsifying transformations on signals.

Those three parameters could be compared to higher-order cumulants or  moments of a distributions. They are related to integrals of the shape:
$$M_d=∫^∞_{−∞}(t−c)df(t)dt$$
or their ratios. A related question is Finding the right measure to compare sound signals in the frequency domain.


Connection with system analysis and laplace&Z transform Laplace&Z transform is just frequency analysis of a system with a multiplication of decaying exponential.  We analyse frequency with varying Laplace Exponential, which can be seen in the formula below for laplace transfrom.

From my simple aspect I can't make the theoretical connection, how inspecting frequency  with decaying exponential helps with understanding how system will behave? And how it helps with coping with differential functions for example "Analysis of Electric Circuits". And even with IIR filters.

What further surprises me the things laplace&Z transforms used. They seem to me completely different .But somehow multiplying the system with decaying exponential works for them all.
The only similarity I see between; for example RLC curcuits and IIR filters is; these systems are not memoryless(output depends on prev. outputs).

When I search the net I see only laplace-table's, questions etc.. I don't want to solve a problem, I want to understand the problem. If somebody can explain or direct me to a link I will be thankful.

You are right that the (bilateral) Laplace transform can be interpreted as the Fourier transform of $e^{-\sigma t}f(t)$. However, I think that the significance of the Laplace transform only becomes clear when $s=\sigma+j\omega$ is viewed as a complex variable because then we can study the analytic properties of the system function. E.g., electrical networks are characterized by the poles and zeros (in the complex $s$-plane) of their system function.

The practical importance of the Laplace transform for the analysis of continuous systems lies in the fact that it transforms linear differential equations with constant coefficients into algebraic equations, which are obviously much easier to solve. This is achieved by the fact that in the Laplace transform domain, integration is equivalent to multiplication by $1/s$, and differentiation becomes multiplication with $s$.

Any linear time-invariant (LTI) system, such as the RLC-circuit you mentioned, will result in linear differential equations with constant coefficients. This is why any LTI system can be analyzed by the Laplace transform. Note that in practice usually the unilateral transform is used (where the integral starts at $t=0$ instead of $t=-\infty$), because practical systems are causal, and because in this way initial conditions (such as charged capacitors at $t=0$) can be easily taken into account. The unilateral Laplace transform is thus an appropriate tool to solve initial value problems.

Everything said so far carries over to the $\mathcal{Z}$-transform if you replace 'continuous system' by 'discrete system', and 'differential equation' by 'difference equation'. Linear time-invariant discrete systems consisting of adders, multipliers and delay elements are treated by the $\mathcal{Z}$-transform in the same way as LTI electrical networks are treated by the Laplace transform. And indeed, IIR filters are analogous to RLC circuits in the sense that both have an infinitely long impulse response and that both can be characterized by the poles and zeros of their system functions.


Approximating Lorentzian Fourier transform with FFT I have a Lorentzian frequency distribution

$F(w) = \frac{1+iz}{1+z^2}$

Where

$z = \frac{w-\Omega}{R}$

With $\Omega$ being the peak frequency and R the decay constant. I know that analytically the Fourier transform should be

$F(t) = exp(i\Omega 2\pi t)exp(-Rt)$

When i take the FFT of this expression, it doesn't return the original frequency plot. I understand that there may be a scaling factor (1/n) that should be in there somewhere, but even when i scale for one frequency, if i then change $\Omega$ or R, the amplitude is no longer scaled properly, suggesting that the scaling factor is a function of $\Omega$ and/or R. The FFT also seems to be mirrored along the intensity axis.

I'm fairly new to DSP, but i do understand that the continuous fourier transform is not the discrete fourier transform. I've read this (https://dspillustrations.com/pages/posts/misc/approximating-the-fourier-transform-with-dft.html) but that approach makes the approximation worse.

I'd like to return the original frequency distribution when i take the FFT of my time signal. Am i missing something fundamental or is it a fairly simple scaling error? I've attached my code below.

Cheers.

# R script to compare FFT and Analytical fourier transform
library(SynchWave)

#-------------------------------------------------
# Frequency and time axes
n <- 100
f <- seq(0, 1, length.out = n)
t <- seq(0, n, length.out = n)

# peak paramaters
O <- 0.3 # Frequency values from 0->1
R <- 0.04 # Decay in arbritrary units

z <- (f-O)/R

# The original lorentzian frequency
ff <-complex(re = 1, im = z)/(1 + z^2)

# creating the time domain signal
ftideal <- exp(-R*t)*exp(complex(i = (O)*2*pi*t))

unscaled <- (fft(ftideal))
scaled <- unscaled - min(Re(unscaled))

plot(f, Re(ff), type = 'l')
lines(f, Re(scaled), type = "l", col = 'red')

```

First of all, take a look at this very related question and its answer(s). Second, your analytic solution is wrong, so it's no surprise that you don't see the correct result.

The time domain function

$$x(t)=e^{j2\pi f_0t}e^{-Rt}\tag{1}$$

has the following Fourier transform:

$$X(f)=\frac{1}{R}\frac{1-j\frac{2\pi(f-f_0)}{R}}{1+\left(\frac{2\pi(f-f_0)}{R}\right)^2}\tag{2}$$

Also, scaling should be multiplicative and not additive (unless you're in the logarithmic domain), and you should also not only look at the real parts of complex-valued functions.

As explained in this answer, approximating the CTFT by the DFT usually results in two types of errors: the truncation error (due to truncation of the time domain function), and the aliasing error (due to sampling the time domain function). These errors can be made small by choosing a large sampling frequency and a large DFT length.

Below is an example code in Octave/Matlab showing how the CTFT of the given function can be approximated by the DFT:

F0 = 0.3;
R = 0.04;

Fs = 30;            % sampling frequency
Ts = 1/Fs;
Tmax = 100;         % length of time domain signal
N = round(Tmax/Ts); % DFT length

t = (0:N-1) * Tmax / (N-1);
N2 = round(N/2);
f = (0:N2-1)/N*Fs;

fun = exp(1i*2*pi*F0*t) .* exp(-R*t);

% analytic expression for CTFT
z = 2*pi*(f - F0) / R;
FTfun = (1 - 1i*z) ./ (1 + z.^2) / R;

% DFT approximation of CTFT
FTfun2 = Ts * fft(fun,N);
FTfun2 = FTfun2(1:N2);

plot(f,20*log10(abs(FTfun)),f,20*log10(abs(FTfun2)))
    title('magnitudes (dB)'), xlabel('f'), legend('CTFT','DFT'), grid on
    axis([0,Fs/2,-40,30])


Simulation of FMCW radar in frequency domain I would like to simulate a FMCW radar in the frequency domain. I used the code (Matlab) I attached below. Specially I am interested in the phase spectrum of the beat signal. In a further research I would like to calculate the response of the FMCW signal when it is hitting a dielectric slab. But for now I want to start with a simple example, a single target in a certain distance.

I believe the code is working, but I have the following problems:

The phase spectrum over the range is depending on the sample frequency. Why is that and how can I avoid that?
I am not quite sure if the phase response over the range is correct, does anyone have an idea what to expect?
This is just an 1-D simulaion, but in the future I want to extend it to 2-D. Is there a way to run the code faster?
Is there another way of attacking this simulation problem in Matlab?

Thanks Pavel

%%
function ProperShiftChirp

clear all; clc;
%generate FMCW Signal
c=299792458;

f_start=65e9;
f_end=90e9;
fc=80e9;

lambda=c/fc;

B=f_end-f_start;
T=4e-6; %Duration for one sweep;
PRI=T;

R=0.5; %distance target [m]
td=R/c; %roundtrip simple [s]

fs=2*B;
slope=B/T;

t=-T/2:1/fs:T/2;

%Generate ChirpSignal
sig_transmit=exp(1i*(pi.*slope.*t.^2));

%frequency domain
a=nextpow2(length(sig_transmit));

xfft=fft(real(sig_transmit),2.^a);

%scaling frequency domain
df=fs/length(xfft);

freq=(0:df:fs-df);

% transfer function
absxfft=abs(xfft);
phasexfft=angle(xfft);

S=absxfft.*exp(1i.*phasexfft).*(exp(-1i*2*td*2*pi.*freq));

%mixing transmitted and received signal
erg=fftshift(fft(ifft(xfft).*ifft(S)));

%single band
idx=(length(erg)/2);
abserg=abs(erg(idx:end));
phaserg=angle(erg(idx:end));

meterabs1=(c*freq(1:idx+1))./(2*slope);

figure;plot(meterabs1,abserg);
xlim([0 R+0.2*R]);
xlabel('metre');ylabel('|S|')
figure;plot(meterabs1,(unwrap(phaserg)));
xlabel('metre');ylabel('rad')

The function of the FMCW radar is equal to the scanning of the environment with a si-pulse. Thats why the following approach can be made:

Formulate reflection transfer function, with bandwidth $\Delta\omega$, center frequency $\omega_c$, reflection coefficient $r$, distance to target $l$, phase constant $\beta=\frac{\omega}{c_0}n_0=kn_0$ i.e.

$$R(\omega)=re^{-j2\frac{\omega}{c_0}n_0 l}= re^{-j\omega \tau n_0 }$$

Calculate the $R(\omega)$ for a wished center frequency and bandwidth in MatLab
Calculate the inverse fourier transformation of $R(\omega)$ to get $R(t)$, zero padding is recommended

You finally get the time domain response for a single target, which is a si pulse at the time $\tau=\frac{2l}{c_0}$ which can be transformed into range by the previous in line equation

Hope I could help some of those who are looking for a start.


How to verify the obtained SNR of the noisy signal after adding noise? I defined the SNR as follows:
$$\rm SNR = \frac{\text{Signal power density}}{\text{Noise power density}}$$
Given noise-free signal, I added the noise to this signal by using the formula above. Now, I am looking for a way to verify if the obtained SNR of the noisy signal is equal to the SNR value that I set when I generated the noise.
Thanks in advance!

Here is a common/straight-forward way using the discrete time domain samples. If you have the noise free signal, $x[n]$, and you created a noise signal, $w[n]$, then you can calculate the SNR by using the formula:
$$\hat{\text{SNR}}=\frac{\sum_n |x[n]|^2}{\sum_n |w[n]|^2}$$
Edit
An unstated assumption above is that $\mathbb{E}\big[x[n]\big]=\mathbb{E}\big[w[n] \big]=0$. In general, you'd use the variance not the sum of squared values, as to take care the possible non-zero means:
$$\hat{\text{SNR}}=\frac{\text{var}(x[n])}{\text{var}(w[n])}$$


A question about the meaning of pole in time domain Lets say I have a transfer function $H(s)$ of a system defined in $s$-domain as:
$$H(s) = \frac{1}{s - (-1-j)}$$

So I conclude that the pole on the $s$-plane is where $s = 1+j$. So far so good.

Now does that mean if the Laplace transform of the input to the system is $s = 1+j$ the system goes crazy/infinity/oscillate?
Does that mean if I take the inverse Laplace transform of $1+j$ I can find what input in time domain makes the system unstable?

Let $H(s)$ be a transfer function of the form
$$H(s) = \frac{1}{s-p}$$
where $p$, which is a pole of $H(s)$, can be written as a complex number $a+jb$. Taking the inverse Laplace transform of $H(s)$ gives the corresponding impulse response $h(t)$ (that is, the output of your system when given $\delta(t)$ as input). Noting $\mathcal{L}^{-1}$ the inverse Laplace transform, we have
$$h(t) = \mathcal{L}^{-1}\{H(s)\} = e^{pt} = e^{at}e^{jbt}.$$
Now let's look at what this impulse response looks like. The term $e^{at}$ is a simple exponential which will be either decaying (if $a < 0$) or growing (if $a$ > 0) with time. The term $e^{jbt}$ will be responsible for oscillations in the output of your system (remember that $e^{jbt} = \cos(bt) + j\sin(bt)$). From this, you can infer the stability of your system and understand why we need poles in the left-hand side of the $s$-plane (i.e. we need $a < 0$) for the system to be stable.

Often, the numerator and the denominator of your transfer function have real coefficients, and in this case poles appear in complex conjugate pairs. You could for example have
$$h(t) = e^{at}(e^{jbt} + e^{-jbt}) = 2e^{at}\cos{bt}.$$

I like to keep this picture in mind (taken from here) which greatly summarizes this.

For more complex transfer function, partial fraction decomposition can be used to go back to simple cases as presented here.


Image Cross-Correlation interpretation I'm working on this idea of choosing best cover image from given db for image steganography(hiding secret image inside cover).
I've started with calculating cross-correlation matrix between images in db
 with given secret image(there are images with different content/texture, like secret could be a building while cover image could be a face and so...).
So first of all i'm not sure if cross-correlation is a proper choice!? but anyway how can i interpret cross-correlation matrix for image similarity? should i only rely on maximum of the result matrix or mean would do better comparison?
What other comparison methods would you recommend for this problem?
I was thinking of binary similarity measures but i'm not sure...(i guess i need proper features first).
Thanks.

I'm using scipy.signal.correlate2d for calculating cross-correlation matrix.

In the case of 2 images $m_{1}(x,y)$ and $m_{2}(x,y)$ of size $m×n$ , where $(x,y)$ are the coordinates of the image, the correlation constructs a matrix $A$ where each $A(i,j)$ is the degree of similarity between the image $m_{1}(x,y)$ and $m_{2}(x+i,y+j)$ (note that $m_{2}(x+i,y+i)$ is the image $m_{2}(x,y)$ translated $i$, $j$ pixels). This degree of similarity is basically a measure of the angle between the vectors $v_{1}$ and $v_{2}$ of $R^{m×n}$ that could be constructed by vectoring the images $m_{1}(x,y)$ and $m_{2}(x+i,y+j)$, this is because for each $i, j$ What the correlation really does is the normalized dot product $\frac{v_{1}.v_{2}}{|v_{1}||v_{2}|}$. Remember that $v_{1}.v_{2}=|v_{1}||v_{2}|cos(θ)$, where θ is the angle between the vectors $v_{1}$ and $v_{2}$.

In conclusion, the maximum $A(i,j)$ tells you that for the degree of similarity (angle θ) between $m_{1}(x,y)$ and $m_{2}(x,y)$ to be maximized, $m_{2}(x,y)$ must be translated $(i, j)$ pixels.

This measure of similarity is not very good, since it is not robust to changes in tones, lighting, noise, among others. There are other comparison alternatives such as PCA, SRC, SIFT or SURF.


Fourier-series odd vs even square wave? I really need to know the difference when doing a fourier seires between even and odd square waves. I've been trying to understand but I just get the same results and the same spectrums... Is the difference in the formula? If so, I can't find any formulas online that match with that my professor gave me. They're all a little different.

For a trigonometric Fourier series:

$$
\tilde{x}(t) = a_0 + \sum\limits_{n=1}^{\infty} a_n \cos(n \omega_0 t) + \sum\limits_{n=1}^{\infty} b_n \sin(n \omega_0 t)
$$
with
$$
a_0 = \frac{1}{T} \int_{-T/2}^{T/2}  x(t) \ dt
$$

$$
a_n = \frac{2}{T} \int_{-T/2}^{T/2}  x(t) \cos (n \omega_0 t) \ dt, \quad n = 1, 2, \ldots
$$

$$
b_n = \frac{2}{T} \int_{-T/2}^{T/2}  x(t) \sin (n \omega_0 t) \ dt, \quad n = 1, 2, \ldots
$$
and $ \omega_0 = \frac{2\pi}{T} $ .

the Fourier series of an odd square wave will be a sum of only sine terms, and the Fourier series of an even square wave will be a sum of only cosine terms.

This is because the product of an odd and an even function is an odd function, which when integrated over a symmetric interval (i.e. from $-T/2$ to $T/2$), will return zero. For an odd square wave, this means that all the $a_n$ will be zero, and for an even square wave, all the $b_n$ will be zero.

You can save yourself time and potential mistakes by exploiting this fact and only computing the coefficients that are nonzero.


Why multiplying frequency kernel with sinusoid results shifting frequency domain I am having great time while reading "The Scientist and Engineer's Guide to Digital Signal Processing" book.

In "Chapter-14 Introduction to Filters" author indicates we can create high pass filter from low pass filters by multiplying low pass filter's frequency kernel with a sinusodial which has 0.5 frequency. Link to the this chapter

I thought I read every chapter carefully but I can't understand why this will work?
Can somebody please explain the reasoning behind of this process as easy that a computer science graduate could understand.

This can be explained by the Convolution theorem

As you state, in the time domain you have the kernel (f) multiplied by a sinusoid (g).

If we take the Fourier transform of f·g (i.e., F(f·g)), this is equivalent to the convolution between F(f) and F(g) (i.e., F(f)*F(g)). Note that the Fourier transform of the sinusoid is proportional to two delta functions located at the frequency and negative frequency of the sine wave. So when you convolve these deltas with the frequency kernel your frequency kernel "shifts" to be centered around the sinusoid  frequency.


Autocorrelation of a stochastic process which is a sum Let's say that a stochastic process is defined as
$$X(t) = \sum_{n = -\infty}^{+\infty} X_n * f(t- nT)$$
where $X_n$ is the $n$-th symbol among independent random symbols. I have proved that $\mathbb E[X(t)] = 0$.  But, I want to prove that this stochastic process is Wide Sense Stationary(WSS) which means by definition that its autocorrelation function $R_{xx}(t+k,t) = \mathbb E[X(t+k)*X(t)] = R_x(k)$ i.e. the auto-correlation function of $X(t)$ depends only on $k$, which is the time difference between $t+k$ and $k$. What I have proved is that the auto-correlation function of $X(t)$ is zero but I am not sure that this proves that this stochastic process is WSS. I have attached a picture of the problem and what I have done so far in order to be clearer.

Your proof is incorrect.  the step indicated below is not true, i.e.:
$$
E[X_n X_m] \not = E[X_n] E[X_m]
$$
in general.


Computer Vision Libraries for PHP Are there any computer vision libraries in PHP, or any server-side programming languages?

If not, maybe we can make a list of some -already implemented- methods (Viola Jones face detector, landmark detector etc.) to help someone who wants to implement such a library?

Someone's written a php interface to OpenCV - https://github.com/mgdm/OpenCV-for-PHP, though I haven't tested it.

As you've tagged face detection I'll also mention that someone has written a php interface specifically to the face detection functionality of openCV too; http://www.xarg.org/project/php-facedetect/

Alternatively you could try (and this would probably offer more flexibility);

[Note: this is just a short list; searching e.g. "python image processing" will give more options]

Python

OpenCV bindings
scipy.ndimage
Scikit-image

Java

OpenCV bindings
ImageJ api

C++

OpenCV
CImg

The top answer to this question also has a good discussion, though I'd be careful about using Matlab in this instance as you might fall foul of their licence restrictions.


How to find system function, H(z) in the z-domain, given zero-pole plot of the system? I'm given a zero-pole plot (Real/Imaginary Plane) of a system in which poles occur at $\frac{-9}{10}$, $\frac{9}{10}$, and zeros occur at $\frac{10}{9}j$, $\frac{-10}{9}j$. An additional zero occurs at negative infinity.

I'm also given that $H_{F}(f)=1$ for $f=0$, where $H_{F}(f)$ refers to $H(z)$ for when $z=e^{-j 2\pi  f}$. I'm tasked to find $H(z)$.

I know I can express $H(z)$ as $$H(z)=G_0z^{-1}\frac{(z-z_0)(z-z_1)\ldots}{(z-p_0)(z-p_1)\ldots} $$ for zeros $$z_0,z_1,\ldots$$ and poles $$p_0, p_1,\ldots$$

I can find the constant $G_0$ by using the fact  $H_{F}(f)=1$ for $f=0$.

I'm unsure how to express the fact that the system has a zero at negative infinity.

The $z^{-1}$ term in the $H(z)$ expression you provided is the expression of the zero at infinity:

$$\lim_{z\to\infty} \frac{1}{z} = 0$$

And in your specific case:

$$\lim_{z\to\infty} H(z) = G_{0}\frac{1}{z}\frac{(z-\frac{10}{9}j)(z+\frac{10}{9}j)}{(z-\frac{9}{10})(z+\frac{9}{10})} = 0$$

Generally, if the order of the denominator exceeds the order of the numerator, then you have 1 or more zeros at $\infty$


FFT spectre graph measurements y-axis I am very new to this things. Sorry for probably stupid question. I don't understand what units and meaning have the values on Y-axis of Fourier Transform graph? On X-axis it is Frequency (Hz). Pretty clear. Can't find the meaning of Y.

Here is the code in MATLAB I used to produce the graph, and the image of the output:

A = 5;
t = 0:1/1000:1;
f1 = 7;
s = A * sin(2 * pi * fi * t);
F = fft(s);
figure;
plot(abs(F)), (s);

The Y-axis units are the same as the Y-axis units you would have if you plotted the time-domain signal.  For instance, if the samples that you input into the FFT are the signal voltage, then the FFT units are also voltage.  To get energy units you have to square the absolute value of the FFT.


Are there any difference between the FFT magnitude of a signal and average FFT magnitudes of the signal broken into multiple frames? Suppose I calculate the FFT of a 5 seconds sound wave. If I'm only interested in the FFT magnitudes across frequencies, I should obtain the same result if I segmented the sound wave into five 1-second frames, obtain the FFT magnitude for each, and taking the average of the magnitude across the frequencies?

Calculating the FFT across a longer signal would allows for smaller frequency resolution, assuming sampling rate is constant.

Does averaging across the frames offer advantages in reducing white noise in the signal? Or does it not matter in that one would get the same result with a long FFT? When might one uses the average magnitudes across frames?

Thanks

You would consider sectioning the overall data into smaller FFT frames for the possible reasons:

you are interested to know that for certain durations of time (length of ye FFT sections) which frequencies were active
there are not enough hardware cycles or resources to carry out a bigger FFT

If you are only interested in knowing what frquencies were ever present in the signal in it's entire duration and you have computation resources to carry out bigger FFTs then yes, you can carry out the FFT of the entire signal itslef. The bigger the sample size of signal duration in time domain over which the FFT is taken the better the resolution in frequency.

For the case you mention you would need to average out by 5 the result of the bigger FFT


How to Linearise an Input Image before applying Color Correction? I'm doing color correction on Raw images taken from Camera directly. I want to make my Camera reproduce same colors as my Target camera.

Here are the steps I'm following:

Take Target Color values by Capturing a Macbeth Color chart by using Target Camera & Take set of 24 Color Palette's.
Take the Input Color Values by Capturing the same Macbeth Color chart by using Reference Camera & Take set of 24 Color Palette's.

Now if i Compute Color correction directly at this stage, all the White colors are appearing Rose or Not correct.

So I'm applying Gamma correction as follows:

Split the Three channels. For each Color find the Gamma for that color palette by using this formula.

float Gamma_R = log10(Target_R/255.0))/(log10(Input_R/255.0);

Error Cases:1

If the Input_Channel Value == 255 make it as 254

If the Input_Channel Value == 0 make it as 1

If gamma_Channel > 3 or gamma_Channel < 0.2 make gamma_Channel as 1.

Average the Values of Gamma_R for all 24 colors & make this as a
gamma of that channel.
Apply Gamma for Each channel using the Gamma computed for Each channel using the formula.

Corrected_R = 255 * (Input_R/255)^(1/Gamma_R)

My Problem:

How to make sure the Gamma Values that I've computed are correct at this stage? Also correct me if I'm doing it wrong.
What to do when the Following Error case: 1 happens?
After applying Gamma correction & Color Correction using 3x3 Matrix (Assume I'm doing that correctly) still the colors are not reproduced correctly.
If i encode the gamma by using same Gamma Values from Previous stage using this formula, again the Colors are not reproduced correctly.

Corrected_R = 255 * (Input_R/255)^(Gamma_R)
So do i need to compute the Gamma again from the Color corrected output?

Any suggestion or References is greatly appreciated!

Your linearization gamma calculation does not take into account the correction matrix.

You could do a combined numerical optimization of the pre-correction-matrix gamma adjustment, the correction matrix, and the post-correction-matrix gamma adjustment, preferably using a colorimetric error metric. There are 11 variables to optimize (or up to 14 if you wand to add constant terms). That's doable without too much work if you don't need to repeat this very often. I have a C++ Differential Evolution optimization library here that should be able to handle it: http://yehar.com/blog/?p=643


Finding if data in array is periodic I have a very large array (about 500.000 in size). I've plotted the array in chunks and founded at the tail of the array there is a periodic behavior, so I want to know if there is a way to know if it is periodic an if so what is its period?

I don't have any idea about signal processing so I need to know if there is a simple and straight forward method to check it. I'm using Python and Mathlab.

Thanks!

Mathematically, what you probably want is the autocorrelation function. That function is a function on a function over time or a sequence, which gives you info about how similar your signal is to itself when shifted by shift $\tau$.

If your signal is periodic with period $\tau$, then it will be very similar to itself when shifted by $\tau$, $2\tau$, $3\tau$, and so on. So, you'll get $\tau$-periodic peaks in your autocorrelation function

A very common thing to do is to look at the Power Spectral Density (PSD), which is, in fact, the Fourier transform of the autocorrelation function. You can, however, in your sampled case, calculate the PSD directly by Fourier-transforming (Matlab function fft) your sequence and finding the abs² of it. The peaks correspond to periodicities.


Can the Entropy Be Used as a Measure of the Contrast of an Image? I was looking for value that can measure a contrast of an image. And I found a couple of answers like this one https://stackoverflow.com/questions/13397394/how-to-measure-contrast-in-opencv-visual-c/13398126.

But I can't find any explanations WHY is that. No articles or links.
So I wrote this code:

def shannon_entropy(img):

    histogram = img.histogram()
    histogram_length = sum(histogram)

    samples_probability = [float(h) / histogram_length for h in histogram]

    return -sum([p * math.log(p, 2) for p in samples_probability if p != 0])

And found the Shannon entropy for four images:

 - original image "image_1". Entropy is 7.82426730799067

 - more contrast image "contrast". Entropy is 5.682661230050924

 - much more contrast image "contrast2". Entropy is 2.1507540571698587

 - less contrast image "contrastLow". Entropy is 6.565885023104716

I thought, that for the last image entropy will be the highest, but no. And it completely baffled me.

Can somebody explain:

Why entropy is a measure of contrast?
What's wrong with my images or calculations?

Hard facts first.

Why entropy value can be used as measure of contrast of an image?

It can't. Entropy is not a function of contrast.
I'll show why. Time for a bit of basic Information Theory¹:
What entropy is, and what happens to it under transformations
Definitions
Let's define $H(·)$ as the entropy of a source; imagine $X$ is your source picture, $Y$ is the picture after you subjected it to some operation $F:\, F(X)=Y$ (e.g. contrast spreading, or compression, or adding random noise).
Intuition of Entropy
$H(X)$ is a measure for the amount of information in the source image. It's hence a measure for the "uncertainty" of the value of each pixel; if that picture was just one intensity all over, the uncertainty would be low, and so the information, and thus the entropy.
Let's assume your input image has enough contrast to at least make use of 8 bits of dynamic range. Then, the highest entropy you could get would be 8 bit per pixel – in case all pixel values are present equally likely. That's how I know your entropy calculation is broken – a black and white image only has two values, and the discrete uniform distribution over 2 values has exactly 1 bit of entropy, and any other distribution has less. You'll need to fix your entropy calculation!
Contrast
However, contrast doesn't care about the full distribution – there's different contrast definitions/measures, but they typically only consider the first and second statistical moment. Many distributions can have the same variance – but only one is the uniform distribution that maximizes entropy. From that alone you can see that while there is a relation between contrast and entropy, it's not a direct relation:
You just can't have high entropy with very low contrast. However, high contrast doesn't say much about entropy.
Applying image transformations
Consider your contrast enhancements:
$H(X|Y)$ is the amount of the uncertainty that remains about the original image when you apply the operation $Y=F(X)$.
When your operation $F$ is reversible (e.g. if you just invert the pixel values, or if you scale them by 2 without clipping because an original value was larger than half the maximum pixel value), then there's nothing unclear about the original image (just apply the "reverse" operation), and thus $H(X|F(X)=Y)=0$. If you think about it, it's pretty intuitive: A revertible operation doesn't change the probabilities, it just changes the values that have these probabilities.
Information in the face of change
We call the information we get about the source when looking at our transformed image the mutual information between $X$ and $Y$, $I(X;Y)$. It's the difference between the source $X$'s entropy and the conditional entropy $H(X|Y)$, i.e. $I(X;Y) = H(Y) - H(Y|X)$!

In the case that the operation is revertible, as you can see, $H(X|Y)=0\implies I(X;Y)=H(X)$.
And, since we know the deterministic operation $F$, the uncertainty about what $Y$ looks like if we know $X$ is $H(Y|X)=0$, too, so that $H(Y)=I(X;Y)$, and conclusively, $H(Y)=H(X)$. In other words, when you do an invertible operation, you can't increase or reduce entropy.
Now, there's the thing: Your contrast operations aren't revertible, simply because

your image has finite bit depth, so some rounding happens and different values in $X$ get mapped to the same value in $Y$,
your image has finite bit depth, so there is some clipping when you reach the extreme values.

That means that you get a $H(X|Y) > 0$; from that we can see that $I(X;Y)<H(X)$. However, since the operation you apply is still deterministic (i.e. you can calculate $Y=F(X)$ based on $X$, it's always the same), the uncertainty about $Y$ knowing $X$ is zero, $H(Y|X) = 0$, and therefore
$$H(Y)=I(X;Y)+H(Y|X)=I(X;Y) + 0 = I(X;Y) < H(X),$$
which is exactly what you see: all your contrast-stretching and -reducing operations decrease the entropy compared to the original image.
What we've learned today

Having high contrast is a necessary, not a sufficient condition for high entropy
Entropy is limited by the number of values occupied
Entropy is invariant to revertible (that implies deterministic) (image) operations
Entropy is reduced by deterministic non-revertible operations

Corollary:

if you want to increase an image's entropy, you'll need to do something non-deterministic, like adding noise.
Information theory is fun and you should definitely study it.

PS.
Here's an entropy function that makes sense; I don't know what kind of data type img is, but I assume you can just go through all its values (it's iterable).
from collections import Counter

def shannon_entropy(img):
    ld = lambda x: math.log(x, 2)
    value_counter = Counter(img)
    N = len(img)
    H = 0
    for count in value_counter.values(): #this gives the counts of values
        p = float(count)/N
        H = H - p*ld(p)
    return H

 ¹ yes, that is a pun.


Octave: Finding the mean fundamental frequency given an audio signal I have an audio .wav file of someone's speech of which I want to find the mean fundamental frequency. Till now, all I have managed to do is load the data into a vector and (roughly) filter the silence out. I referred to a MATLAB example of finding the fundamental frequency, but it lies way out of the range of the mean fundamental frequencies which I have in a dataset I am referring to. The code on MATLAB is so-

function F0 = fundafreq(data)
  % Copied from matlab site
  F0 = 0;
  Fs = 48000;  % My sampling frequency
  [b0,a0] = butter(2,325/(Fs/2));
  xin = abs(data);
  xin = filter(b0,a0,xin);
  xin = xin-mean(xin);
  x2 = zeros(length(xin),1);
  x2(1:length(data)-1) = xin(2:length(data));
  zc = length(find((xin>0 & x2<0) | (xin<0 & x2>0)));
  F0 = 0.5*Fs*zc/length(data);
end

This gives me the fundamental frequency of around 40 Hz whereas my data set has mean fundamental frequencies not more than 1. I am only a beginner when it comes to digital signal processing, so I am not completely aware whether this code is giving me the answer or whether I need further calculations.

Also, it would be preferable if I could calculate this without needing the sampling frequency.

Your code appears to be looking for zero crossings of the rectified signal, which won’t work reliably if your signal has strong harmonics or overtones of the fundamental frequency within your filter’s passband.  Strong overtones of low frequencies is very common in audio data.

As a likely better alternative, look up pitch detection/estimation algorithms, such as autocorrelation, ASDF and AMDF


Spectral Leakage: Number of Samples Is there a relationship for calculating the number of samples ($N$) that must be used in order to minimize spectral leakage of a discrete time signal.

For example, given a signal

$$x[n] = \sin(2\pi f_0 n T)$$

and $f_0 = 5, T = 0.02, n = 0, \ldots, N-1$, and no information regarding the sampling frequency, how do we choose the value for $N$?

I saw somewhere that $N$ would be equivalent to the sampling frequency, yet again, assuming that isn't available, is there a way of obtaining $N$?

The DFT of a discrete-time sinusoid will only have components at the sinusoid's frequency if there are an integer number of periods inside the window. I.e., you must figure out the period $M$ of the sinusoid and choose the number of points $N$ as an integer multiple of $M$ ($N=kM$). I'm sure you can figure out the period of the given sinusoid.

Also have a look at this answer to a related question.


Question regarding 16 QAM digital modulation Below is an example of a digital modulation system using 16-QAM. I have copied this example from MATLAB communications toolbox,

   M = 16;                     % Size of signal constellation
   k = log2(M);                % Number of bits per symbol
   n = 30000;                  % Number of bits to process
   numSamplesPerSymbol = 1;    % Oversampling factor
   dataIn = randi([0 1],n,1);
   dataInMatrix = reshape(dataIn, length(dataIn)/4, 4); % Reshape data into binary 4-tuples
   dataSymbolsIn = bi2de(dataInMatrix);
   dataMod = qammod(dataSymbolsIn, M);

I am wondering if

1) there is a need to reshape the bits into tuples of 4 bits

2) there is a need to do the mapping from binary to decimal and then do the QAM modulation, couldn't I just have done

dataMod = qammod(dataIn, M);

Thanks.

First of all, I recommend getting familiar with Matlab's doc command. doc qammod has the answers to your second question: qammod requires a vector of integers between 0 and M-1.

Regarding your first question: the reason to reshape your bits into groups of four is that 16-QAM transmits $\log_2(16)=4$ bits per symbol. So, your bi2de command requires binary numbers of four bits to get decimal integers betwen 0 and 15.

If you want to save lines of code, what you could do is generate integers between 0 and 15 directly, with dataSymbolsIn = randi([0 15], n/4, 1). However, doing this has one disadvantage: calculating the bit error rate at the receiver is more complicated if all you have are the integers. Let's say that you transmit a 15 and receive a 13 (obtained with qamdemod), how many bit errors there were? This is why it's common practice to start with actual bits in the transmitter and compare them one-to-one with the received bits.


Time Series processing using fft I have a set of real data (timestamp and value) with an unstable step between samples (5sec, 30sec etc.).
The data is the % of fillage of a vehicle's tank through time.
Due to the harsh volatility of the data, I need to perform fft to transfer my time series to the frequency domain, select a cutoff point to remove all the noise and then transfer back to the time domain.
Is my thought process correct?
I want to use matlab but I cannot figure out what I have to do with my data. The first 10 rows (I have a document with 10ths of thousands of data) are as follows:

I have produced some dummy timestamp data with the same values and with time step at 0.01sec.
After writing the code in matlab I get the following:

My code is:
load input.txt;
plot(input);
figure;
Fs = 1/0.01;
Ts = 1/Fs;
dt = 0:Ts:5-Ts;
x = input(:,1);
y = input(:,2);
nfft = length(y);
nfft2 = 2*nextpow2(nfft);
ff = fft(y,nfft2);
fff = ff(1:nfft2/2);ff
plot(abs(fff));

If you want some data don't hesitate to contact me.
Thank you in advance!

Is my thought process correct?

No. Frequency domain filtering is difficult and even if you get it right, I sincerely doubt that low-pass filtering will solve your problem
Your data is very noisy. I looks like there are two noise sources: one with a "quantization" step of about 20 plus some smaller noise overlaid. Hard to tell without further analysis.
You haven't told us yet, what you actually want to get out of the analysis, so it's hard to give advise. Low-pass filtering will smear out your "refuel" transients a lot. I'm guessing, your best shot would be to build a parametric model and do a piece-wise least squares error fit, but again, you need to tell us what you actually want to achieve.


What is effect of fractional sample delay effect on FIR Filters? In the DSP lecture slides we have this slide:

And

So why this part of FIR Type 2 and Type 4 filter impulse response equation bolded in this slide?

Thanks for your attention.

To answer your primary question, what is the effect of fractional delay- in this case the fractional delay when it would exist (for Type 2 and Type 4 filters only) would be M.5 samples where M is an integer. (So would always be on a half sample increment).

In the picture below the blue line represents the sampled waveform that would pass through the filter, with the blue dots representing the input of the filter, while the output would have samples on the red dots (with the additional M sample delay). This is the effect of a fractional delay.

The presentations are simply showing you the delay of a linear phase filter; in all cases the delay is $N/2$ where N is the order of the filter, which is 1 less than the length L. If the filter has an odd length it will have an even order as in the case of Type 1 and Type 3 filters and the delay will be an integer number of samples. If the filter has an even length it will have an odd order as in the case of Type 2 and Type 3, and therefore the delay will be a half sample more than an integer number of samples.

Consider how this works which also illuminates why these four filter types are linear phase.

The frequency response for a FIR filter is determined by replacing evaluating $H(z)$ as z is swept over the unit circle. ($z=1$ corresponds to "DC" while $z=1$ corresponds to $\pm F_s/2$ where $F_s$ is the sampling rate.) To keep z on the unit circle, $z=e^{j\omega}$ as $\omega$ is swept from $-\pi$ to $\pi$ corresponding to $=F_s/2$ to $F_s/2$.

Therefore, we are evaluating $H(z) = H(e^{j\omega})$

Also note Euler's identities:

$e^{j\phi} + e^{-j\phi} = 2\cos(\phi)$

$e^{j\phi} - e^{-j\phi} = 2j\sin(\phi)$

Note that a real waveform such as $\cos(\omega n)$ has a phase that can only be $0$ or $\pi$ radians (mod $2\pi$).

Similarly the imaginary waveform such as $2j\sin(\omega n)$ has a phase that can only be $\pm pi/2$ (mod $2\pi$).

It is only complex waveforms such as $e^{j\omega n}$ with both real and imaginary components that can have arbitrary phases. And note that the magnitude of $e^{j\omega n}$ is always one, so this expression describes phase only.

From this, the approach is to factor the linear phase frequency response into the form $A(\omega)e^{j\phi(\omega)}$, where $A(\omega)$ is real and only describes the magnitude while $e^{j\phi(\omega)}$ describes the phase:

Given the general form of the z-transform for an FIR filter as

$$H(z) = \Sigma_{n=0}^{N-1}h[n]z^{-n}$$

And the frequency response as evaluated on the unit circle as:

$$H(e^{j\omega}) = \Sigma_{n=0}^{N-1}h[n]e^{-j\omega n}$$

If we factor out of this $e^{-j\omega N/2}$ the result is symmetric terms of $e^{-j\omega n}$, resulting in sine and cosine terms.

Here is a simple example of a 5 tap unity gain filter (coefficients = [1 1 1 1 1]). The is a length = 5 filter that is 4th order (even).

The frequency response of this is given as

$$H(e^{j\omega}) = 1 + e^{-j\omega}+e^{-j2\omega}+e^{-j3\omega}++e^{-j4\omega}$$

as $\omega$ is swept from $0$ to $2\pi$ or alternatively $-\pi$ to $\pi$ depending on how you want to present the frequency axis (corresponding to DC to $F_s$ or $-F_s/2$ to +$F_s/2$.)

If we factor out the linear phase component, we easily get the desired $A(\omega)e^{j\phi(\omega)}$ form since we can replace the remaining symmetric exponential terms with sines or cosines using Euler's identities. In this case, in order for the exponential exponents to be symmetric resulting in real sinusoidal terms, we need to factor out what will result in an integer delay of samples:

$$H(e^{j\omega}) = e^{-2\omega}(e^{-2\omega} + e^{-j\omega}+e^{j\omega}+e^{j2\omega})$$

Notice the symmetric exponential components such that we can now convert them easily using Euler's identities to cosine terms in this case:

$$H(e^{j\omega}) = e^{-2\omega}[(e^{-2\omega}+e^{j2\omega}) +  (e^{-j\omega}+e^{j\omega})]$$

$$= e^{-2\omega}(2\cos(2\omega)+ 2\cos(\omega))$$

So here we see that the term $(2\cos(2\omega)+ 2\cos(\omega))$ contributes nothing to the phase except for $0$ and $\pi$ radian phase changes, so only changes the magnitude in the frequency response. Similarly  the term $e^{-2\omega}$ only contributes to phase as the magnitude is always 1. Notice as the frequency term $\omega$ is changed, the phase is changing in direct proportion, in this case at a negative rate of twice the frequency-- this is consistent with the delay of two samples, and is "linear phase", as the phase is changing linearly with frequency.

We can only extract a single linear phase term and get real sinusoids when the filter coefficients are symmetric or anti-symmetric as in the four linear phase types. Further we see how when we have an odd number of taps with asymmetric coefficients how the center tap must be zero as well as other interesting properties that are revealed. Such as Type 2 filters must always have a zero at z=-1 (low pass) and similarly Type 4 filters must always have a zero at z = 1 (high pass).

It may be helpful to know that time delay is the negative slope (derivative) of phase with respect to frequency--- so if the phase versus frequency response is decreasing linearly versus frequency- this would correspond to all frequency components experiencing the same delay and thus there would be no distortion introduced.

So the reason for showing the fractional delay in those slides (and specifically it will be a 1/2 sample fraction) for the type 2 and type 4 filters is to simply show that this will always be the case for those filter types, and why. So the waveform at the output of a filter that has a half-sample fractional delay will have sample locations that are at the midpoints of the input samples for the underlying waveform that could also be represented in continuous time.


Inverse Short time fourier transform I am I am working with 'wav' files in audio signal processing. I am able to find out STFT by using the function " spectrogram(file name, window length). Now I want to apply Inverse Short time fourier transform on the above signal. I tried "invspectrogram". but its showing error as Undefined function for input argument type of double. So please tell me how to proceed.

A magnitude-only spectrogram is not invertible. The spectral phase required for reconstruction has been lost.

A complex FFT (a 1D column of spectrogram data before the absolute magnitude is taken and colored) might be (within numerical limitations) invertible by using an IFFT, but that won't recover the original signal if the STFT is using windowed FFTs.


wavelength filtering (photos) Preamble: I don't know anything about photo-editing or signal processing. But I am a composer who happens to be working on a project in which sound / light overlap.

Question: Is it possible to filter a given image through a spectrum of pre-determined wavelengths? (i.e. Passing an image through a filter in which only 480nm, 560nm, 678nm, 714nm are able to pass through? The numbers here are completely arbitrary.)

Many thanks.

If the image is given in the conventional format of sampled pixel values of an sensor matrix, then no it's not possible to apply any software means of EM wavelength/frequency processing on those pixel data. The reason is clear: the sensor pixels are not capturing time variation of electromagnetic wave, they only capture the static (or average) intensity at a point and collectively the intensity distribution along the sensor surface points. Wavelength (and frequency) sensing at a point requires a time domain capture, such as performed by the RF radio antenna or a microphone (sensing the acoustic pressure wave), which therefore enables subsequent frequency/wavelength processing possible.

Therefore practically there remains only a hardware lens filter possbility which is indeed a trusted parctice in the field of photography.

Look for image filters from optics part producers, they can supply any type of filter you want (they are a litle pricey however)


Algorithms for computing FFT in parallel I am trying to parallelize the computation of an FFT on terabyte-sized signal files. Right now such an FFT using an open-source library takes many hours, even running through CUDA on the fastest GPU I have. The framework I am trying to adapt to this process is Hadoop. In very basic terms, Hadoop distributes a problem over any number of server nodes in the following manner:

  • You split your input file into (key, value) pairs.
  •   These pairs
  are fed into a “Map” algorithm, which transforms your (key, value)
  pairs into some other (key, value) pairs based on what you put inside
  the Map.
  •   The framework then collects all the (key, value)
  outputs from the Maps and sorts them by key, as well as aggregating
  values with the same key to a single pair,  so you end up with (key,
  list(value1, value2, ..)) pairs
   •  These  pairs are then
  fed into a “Reduce” algorithm, which in turn outputs more (key, value)
  pairs as your final result (written to a file).

There are many applications for this model in practical stuff like processing server logs, but I am having a hard time applying the framework to chopping up an FFT into “map” and “reduce” tasks, especially since I am not really familiar with DSP.

I will not bother you with the programming mumbo jumbo, as this is a DSP Q&A. I am, however,  confused on what algorithms exist for computing FFTs in parallel; Map and Reduce tasks can’t (technically) talk to each other, so the FFT must be split into independent problems from which the results can somehow be recombined at the end.

I have programmed a simple implementation of Cooley-Tukey Radix 2 DIT that works on small examples, but using it for recursively calculating odd/even index DFTs for a billion bytes will not work.  I have spent a few weeks on reading many papers, including one on a MapReduce FFT algorithm (written by Tsz-Wo Sze as part of his paper on SSA multiplication, I cannot link more than 2 hyperlinks) and the “four-step FFT” (here and here), which seem similar to each other and to what I am trying to accomplish. However, I am hopelessly bad at mathematics, and applying any of those methods by hand to a simple set of something like {1,2, 3, 4, 5, 6, 7, 8} (with all imaginary components being 0) gives me wildly incorrect results. Can anyone explain an efficient parallel FFT algorithm to me in plain English (one that  I linked or any other) so that I may try and program it?

Edit: Jim Clay and anyone else who may be confused by my explanation, I am trying to do a single FFT of the terabyte file. But I want to able to do it concurrently on multiple servers in order to speed up the process.

I think your main problem isn't how to parallel the algorithm (which actually can be done) but it's numeric precision.
FFTs of a size that large are numerically quite tricky. The FFT coefficients are of the form
$e^{-j\cdot 2\cdot \pi \cdot \frac{k }{N}}$ and if N is very large the coefficient calculation gets noisy. Lets say you have $N = 2^{40}$ and you use 64 bit double precision arithmetic. The first 1000 coefficients have a real part that's exactly unity (although it shouldn't be that way) so you will need higher precision math, which is very inefficient and cumbersome to use.

You will also amass lots of rounding and truncation errors as the sheer number of operations that go into a single output number is very large too. Due to the "every output depends on every input" nature of the FFT, error propagation is rampant.

I'm not aware of an easy way to work around that. Your request is an unusual one. Most applications that do spectral analysis of large data sets do a running analysis where you don't have that problem. Perhaps if you can describe your application and it's constraints a but more, we can point you to a more suitable solution.


When doing a Hilbert-transformation, why not simply multiplying by an exponential? When doing the Hilbert-transformation, which is supposed to be a phase-shift by $\\-\frac{\pi}{2}$, why not simply use $\\e^{-j\frac{\pi}{2}}$ instead of $\\H(f) = -j \cdot \textrm{sgn}(f)$?

In contrast to Jason R's answer I claim that the Hilbert transform is a phase shift by $-\pi/2$ for real-valued signals. By definition, a phase shifter shifts the phase of a sinusoidal signal by some given phase $\phi$:

$$x(t)=\cos(\omega_0t)\quad\Longrightarrow\quad y(t)=\cos(\omega_0t+\phi)\tag{1}$$

Since

$$\cos(\omega_0t)=\frac{e^{j\omega_0t}+e^{-j\omega_0t}}{2}\tag{2}$$

and

$$\cos(\omega_0t+\phi)=\frac{e^{j\omega_0t+\phi}+e^{-j\omega_0t-\phi}}{2}\tag{3}$$

it is clear that in order for the signal $x(t)$ to experience a phase shift by $\phi$, the positive frequencies need to be shifted by $\phi$ and the negative frequencies by $-\phi$. For $\phi=-\pi/2$ that's exactly what a Hilbert transformer does.

A constant transfer function like $H(f)=e^{-j\pi/2}=-j$ doesn't really do much apart from scaling the input signal, in this case by the imaginary unit. So the signal $x(t)=\cos(\omega_0t)$ is "transformed" into $y(t)=-j\cos(\omega_0t)$. Of course, for the special case of the complex-valued input signal $x(t)=e^{j\omega_0t}$ multiplication by $-j$ achieves a phase shift by $-\pi/2$, but usually we want to shift the phase of real-valued signals, and in that case we need a frequency-dependent transfer function.


Flipping the impulse response in convolution During convolution on a signal, why do we need to flip the impulse response during the process?

Adapted from an answer to a different question (as mentioned in a comment) in the hope that this question will not get thrown up repeatedly by Community Wiki as one of the Top Questions....

  There is no "flipping" of the impulse response by a linear
  (time-invariant) system.
  The output of a linear time-invariant system is
  the sum of scaled and time-delayed versions of the
  impulse response, not the "flipped" impulse response.

We break down the input signal $x$ into a sum of scaled
unit pulse signals. The system response to the unit pulse signal
$\cdots, ~0, ~0, ~1, ~0, ~0, \cdots$ is
the impulse response or pulse response
$$h[0], ~h[1], \cdots, ~h[n], \cdots$$
and so by the scaling property the single input value $x[0]$,
or, if you prefer
$$x[0](\cdots, ~0, ~0, ~1, ~0,~ 0, \cdots)
=  \cdots ~0, ~0, ~x[0], ~0, ~0, \cdots$$
creates a response
$$x[0]h[0], ~~x[0]h[1], \cdots, ~~x[0]h[n], \cdots$$

Similarly, the single input value $x[1]$ or creates
$$x[1](\cdots, ~0, ~0, ~0, ~1,~ 0, \cdots)
=  \cdots ~0, ~0, ~0, ~x[1], ~0, \cdots$$
creates a response
$$0, x[1]h[0], ~~x[1]h[1], \cdots, ~~x[1]h[n-1], x[1]h[n] \cdots$$
Notice the delay in the response to $x[1]$.  We can continue further
in this vein, but it is best to switch to a more tabular form
and show the various outputs aligned properly in time. We have
$$\begin{array}{l|l|l|l|l|l|l|l}
\text{time} \to & 0 &1 &2 & \cdots & n & n+1 & \cdots \\
\hline
x[0] & x[0]h[0] &x[0]h[1] &x[0]h[2] & \cdots &x[0]h[n] & x[0]h[n+1] & \cdots\\
\hline
x[1] & 0 & x[1]h[0] &x[1]h[1] & \cdots &x[1]h[n-1] & x[1]h[n] & \cdots\\
\hline
x[2] & 0 & 0 &x[2]h[0] & \cdots &x[2]h[n-2] & x[2]h[n-1] & \cdots\\
\hline
\vdots & \vdots & \vdots & \vdots & \ddots & \\
\hline
x[m] & 0 &0 & 0 & \cdots & x[m]h[n-m] & x[m]h[n-m+1] & \cdots \\
\hline
\vdots & \vdots & \vdots & \vdots & \ddots
\end{array}$$
The rows in the above array are precisely the scaled and
delayed versions of the impulse response that add up to
the response $y$ to input signal $x$.
But if you ask a more specific question such as

  What is the output at time $n$?

then you can get the answer by summing the $n$-th column to get
$$\begin{align*}
y[n] &= x[0]h[n] + x[1]h[n-1] + x[2]h[n-2] + \cdots + x[m]h[n-m] + \cdots\\
&= \sum_{m=0}^{\infty} x[m]h[n-m],
\end{align*}$$
the beloved convolution formula that befuddles generations of students
because the impulse response seems to be "flipped over"
or running backwards in time. But, what people seem to forget
is that instead we could have written
$$\begin{align*}
y[n] &= x[n]h[0] + x[n-1]h[1] + x[n-2]h[2] + \cdots + x[0]h[n] + \cdots\\
&= \sum_{m=0}^{\infty} x[n-m]h[m],
\end{align*}$$
so that it is the input that seems "flipped over" or running
backwards in time!
In other words, it is human
beings who flip the impulse response (or the input) over when computing
the response at time $n$ using the convolution formula, but
the system itself does nothing of the sort.


Determine sequence's value from equally null circular convolution DFT Can anyone help with this exercise from my Signal Processing exam? The text is:

  Given two sequences $x$ and $y$ whose supports are respectively $[0,1]$
  and $[0,3]$ and their circular convolution's DFT is equally null, if
  the first value of $x$ sequence is $x(0)=2$ and the first value of $y$
  sequence is $y(0)=1$ then obtain $x(1), y(1), y(2), y(3)$

I have tried to solve that supposing respectively that $X(k)=0$ for a first case, $Y(k)=0$ for a second case and then both equal to zero, and I have found two solution for each case. The problem is that my professor tell me that it has just one solution and it can be solve using DFT properties in, more or less, 1 minute, without solving linear system or doing hypothesis on the DFTs but I really don't understand what is referring to because I don't find anything on my notes or on my book.

Is there anyone that have 1-minute solution? :S

Your professor is wrong, the solution is not unique. There are two (real-valued) solutions:

x1 = [2, 2, 0, 0]; y1 = [1, -1, 1, -1]

with DFTs

X1 = [4, 2-2i, 0, 2+2i]; Y1 = [0, 0, 4, 0]

(the product of which is obviously zero)

and

x2 = [2, -2, 0, 0]; y2 = [1, 1, 1, 1]

with DFTs

X2 = [0, 2+2i, 4, 2-2i]; Y2 = [4, 0, 0, 0]

(the product of which is also zero).

This is how you obtain these two solutions: first of all, note that the DFT of the circular convolution of the two sequences is equal to the product of the DFTs of the two sequences, where the individual DFTs are of length $4$. The product is zero for all indices $0\le k<4$ if for each index at least one of the two DFTs is zero. Since I assume that we're talking about real-valued sequences, we only need to consider the indices up to and including Nyquist (i.e. index $k=2$), because the last index $k=3$ is redundant due to the symmetry of the DFTs.

Since none of the two sequences can be identical to zero ($x[0]$ and $y[0]$ are both given and non-zero), only specific DFT bins can be zero, not all of them. For the length $4$ DFT of the zero-padded length $2$ sequence $x[n]$ it is easy to see that its DFT is given by

$$\begin{align}X[0]&=x[0]+x[1]\\
X[1]&=x[0]-jx[1]\\
X[2]&=x[0]-x[1]\\
X[3]&=x[0]+jx[1]\;\;(=X^*[1]\;)\end{align}\tag{1}$$

from which it is clear that with $x[0]\neq 0$, only $X[0]$ or $X[2]$ (but not both) can be zero. $X[1]$ can never be zero (for real-valued $x[1]$), from which it follows that $Y[1]$ must be zero. Since $Y[k]$ cannot be zero for all $k$, $X[0]$ or $X[2]$ must be zero. This leaves us with two possible choices for $x[1]$:

$$x[1]=x[0]\rightarrow X[2]=0\tag{a}$$

or

$$x[1]=-x[0]\rightarrow X[0]=0\tag{b}$$

In the first case, we need a DFT $Y[k]$ which is zero for all $k$ except for $k=2$, in the other case $Y[k]$ must be zero for all $k$ except for $k=0$. Note that $k=0$ corresponds to DC, and $k=2$ corresponds to Nyquist. The DFT of a sequence is zero everywhere except at DC if the sequence is constant, and it is zero everywhere except at Nyquist if the sequence is alternating. So for choice $(a)$ we need to choose $y[n]$ to alternate, and for choice $(b)$, $y[n]$ needs to be constant. These are the only two real-valued options.

EDIT: If you allow the sequences $x[n]$ and $y[n]$ to be complex-valued, there are two additional solutions. From $(1)$ it can be seen that for $x[0]\neq 0$, $X[k]$ can only be made zero for one index $k$. This means that $Y[k]$ must be zero for all indices except for that one index where $X[k]$ is zero. Therefore, $Y[k]$ must be a spectral line, i.e. its corresponding sequence $y[n]$ must be a complex exponential on the DFT ($N=4$) grid:

$$y[n]=e^{j\pi nl/2},\quad l=0,1,2,3\tag{2}$$

For $y[n]$ as given by $(2)$, $Y[k]$ is zero for $k=l$. The solutions $l=0$ (constant sequence), and $l=2$ alternating sequence are real-valued and were given above. The two other solutions are

y3 = [1, i, -1, -i]

($l=1$) and

y4 = [1, -i, -1, i]

($l=3$). The DFT of the corresponding sequences $x[n]$ must have a zero at $k=l$. From $(1)$ we get

x3 = [2, -2i, 0, 0]
and
x4 = [2, 2i, 0, 0]

The DFTs are

X3 = [2 - 2i, 0, 2 + 2i, 4]; Y3 = [0 4 0 0]

and

X4 = [2 + 2i, 4, 2 - 2i, 0]; Y4 = [0 0 0 4]

Both products are clearly zero.

In sum, there are four solutions to the problem, two of which are real-valued.


How to design digital filter from transfer function I have two transfer functions and would like to implement an IIR with coefficients.

What is the process to get from my complex transfer function to IIR coefficients?

\begin{align}
F_W &=\frac{1+j\displaystyle \frac{\omega r}{c}-\frac 13\left(\frac{\omega r}{c}\right)^2}{1+\displaystyle \frac 13 j\frac{\omega r}{c}}\\
F_{XYZ} &=\sqrt{6}\frac{1+\displaystyle\frac 13j \frac{\omega r}{c}-\frac 13\left(\frac{\omega r}{c}\right)^2}{1+\displaystyle \frac 13 j\frac{\omega r}{c}}\\
\end{align}

These are the Gerzon/Craven transfer functions for Ambisonic B-Format audio

Following OP comment requesting more details to Maximilian Matthé's answer:

Just like mentioned, the key to your question is the Bilinear Transform, where you substitute $s$ (or $j\omega$ in your case) with $\frac{2}{T}\frac{z-1}{z+1}$

This transform will have for effect to convert your transfer function from the continuous $s$ domain to the discrete $z$ domain.

I assume that you are fimiliar with the $z$ transform and that you know that a delay in the time domain is equivalent to multiplying by $z^{-1}$ in the z domain.  $Z\{f[x-n]\} = Z\{f[x]\}z^{-n}$

The implementation of a discrete filter (FIR or IIR) will happen in the time domain. Your software will implement something like

$$y[n] = b_{0}x[n]+b_{1}x[n-1]+b_{2}x[n-2]...+b_{N}x[n-N]-a_{1}y[n-1] - a_{2}y[n-2]...-a_{N}y[n-N]$$

In my previous equations :

$x$ is the input
$y$ is the output
$n$ is the index of the actual value in your buffer. $n-1$ being the previous value and so on.

If you bring that equation in the $z$ domain, you will have an equivalent like that.

$$Y = b_{0}X + b_{1}Xz^{-1}+b_{2}Xz^{-2}...+b_{N}Xz^{-N} -a_{1}Yz^{-1}-a_{2}Yz^{-2}...-a_{N}Yz^{-N}$$

After some algebra, you'll get the following form, which correspond to a transfer function (ratio of output on input)

$$\frac{Y}{X} = \frac{b_{0}+b_{1}z^{-1}+b_{2}z^{-2}...+b_{N}z^{-N}}{1+a_{1}z^{-1}+a_{2}z^{-2}...+a_{N}z^{-N}}$$

This transfer function is what Maximilian Matthé answer refers to. In it's condensed form :

$$\frac{Y}{X}=\frac{\sum_{i=0}^{N}b_{i}z^{-i}}{1+\sum_{i=1}^{N}a_{i}z^{-i}}$$

Designing a discrete filter correspond to determining the coefficients of that transfer function. e.g. Finding the values of $a$ and $b$.

The value of your coefficients will change according to your sampling frequency, and you can confirm that considering that the Bilinear Transform introduce the sampling time $T$ in your equation.

You might want to consider the fact that the Bilinear transform is an approximation. If you look at Wikipedia's page, you will find this :

To compensate this approximation, you can substitute $\omega$ with

$$\omega_d = \frac{2}{T}\tan^{-1}(\frac{2\omega}{T})$$

This substitution is called frequency warping and will have a slight effect on the frequency response of your filter, but will most likely be negligible.

Finally, if you have access to matlab and the signal processing toolbox, you might want to have a look to bilinear() which will do all of this in a single command.


How to measure the agreement between to curves? I have values (plotted below) of expected RSSI values over time that I would like to compare with my measured RSSI values.  What I was looking for was a way to quantify it so I can change parameters and be able to compare/contrast different approaches.

It is a hard problem in my mind because I don't know how to compare the signals and yet take into account the large-scale (overall shape) and small-scale (individual fluctuations) of signal.

For instance, here is a plot of one set of signals:

In the image I can see that the red measure signal roughly follows the model, but it also does an OK job of simulating some of the sinusoidal qualities of the model (in some places).  Any thoughts?

<>
In response to pichenettes' comments (which seem reasonable), I took a diff of the two values and plotted the abs(fft(diff)) and got this:

I am not sure what to make of that though.  Since we don't have any actual freqs, I am not sure how to scale the axis, and then if I did, what metric would you use?

If the signals are not aligned, yet you have a clue that they more or less "mean" the same thing, or refer to similar data, you could use dynamic time warping (DTW) algorithm to obtain a better correspondence (then simply taking the value at the same location). You could measure RMS, MSE or whatever you like, using those correspondences. For DTW, you might want to check:
http://en.wikipedia.org/wiki/Dynamic_time_warping

A good way to elaborate this idea is utilized in Earth Mover Distance (EMD), which computes the minimal effort to bring the signals to alignment, as a distance measure. EMD is presented here:
http://ai.stanford.edu/~rubner/emd/default.htm

EMD gives you a direct distance, which you might use for further analysis.


What Is the Difference between RLS, LMS and Wiener Filter? When Is One Preferred Over Another? I'm dealing with a channel equalization problem where the channel is modeled as a WSS process.
I understand LMS utilities a Wiener-like approach, ie it converges to the optimal (wiener) solution.

I understand RLS converges to the normal solution, assuming the statistics of the WSS process are unknown.

I would like to know:

What class of predictors can Wiener, LMS, and RLS be classified within? Is it correct to define Wiener as 'Optimal (in MSE sense)', LMS as 'Stochastic gradient predictors', and RLS as a 'Linear predictor'? And what class of filters?
When is one solution preferred to the other and why?
In particular, when defining the cost function with the MSE, what can we achieve? And when we define it simply with the squared error?

All three are Estimators / Predictors.
All of them try to estimate the coefficients of Linear Filter which minimizes an MMSE Cost Function.

The Wiener filter assumes all data is given and sets the way to calculate the optimal solution.

The LMS and RLS are sequential / on line methods to solve the same problem and given the data is stationary they all will converge to the same solution.

The LMS / RLS can adjust their solution according to data and hence are suitable in cases the model changes in time.

The easiest to implement with the least computing resources needed is the LMS and it is the one you should start with. The LMS works on the current state and the data which comes in.
The RLS, which is more computational intensive, works on all data gathered till now (Weighs it optimally) and basically a sequential way to solve the Wiener Filter.

So, I'd start with the LMS.
If the model is constant and the LMS Convergence rate is good enough for my problem I'd stick with it.
If it seems I need faster convergence rate then I'd go with the RLS.


Up sampling vs down sampling,? Refrence wikipedia article for upsampling and wikipedia article for downsampling

Upsampling , interpolation and expansion are described as synonyms

While

Downsampling , decimation and compression are described as synonyms

Are they mentioned properly as synonyms? Are those 3 terms referring to same operation?

Within the context of Digital Signal Processing...
Upsampling and interpolation are synonyms, expander is not.
Specifically, an upsampler and an interpolator include an expander followed by an anti-imaging lowpass filter. The expander, is the initial stage of an interpolator or an upsampler.
The following is an expander:
$$x[n] \longrightarrow \boxed{ \uparrow 2} \longrightarrow y[n] = \begin{cases} { x[n/2]  ~~, ~~~n = 0,2,4,... \\ ~~~~ 0~~~~~~~ , ~~~n= 1,3,5,...} \end{cases} $$
The following is an upsampler / interpolator:
$$x[n] \longrightarrow \boxed{ \uparrow 2} \longrightarrow w[n] \longrightarrow \boxed{ h[n]}  \longrightarrow y[n] $$
Similarly, downsampler and decimator are synonyms, compressor is not.
The following is a compressor:
$$x[n] \longrightarrow \boxed{ \downarrow 2} \longrightarrow y[n] = x[2n]$$
and the following is a downsampler / decimator:
$$x[n] \longrightarrow \boxed{ h[n]}  \longrightarrow w[n] \longrightarrow \boxed{ \downarrow 2}  \longrightarrow y[n] $$
Specifically, a downsampler and a decimator include an anti-aliasing lowpass filter followed by a compressor. The compressor, is the final stage of a decimator or a downsampler.
Note that you may omit the anti-aliasing lowpass filter when downsampling (or decimating) due to the fact that the input signal is already bandlimited to the band of interest. In such a case, apparently the expander, the downsampler, and the decimator become the same operations; but this is not the general case.
There's no fixed definition of such things but the above is very common.
Authors might prefer using them in different contexts as long they do explicitly describe what they are doing. Informally I may use upsampler and expander in similar places, but the distinction should be clear in my mind.


Ideal Geometric Arrangement of Microphone Array I'm going to have to give a bit of context for this question to make sense.

I am working on a project which includes audio source localisation in 3-D space through TDoA (Time Difference of Arrival - from here called $\tau_{ij}$, between sensor $i$ and $j$) measurements. In other words I have $m$ sensors that each have a slightly different recording of the audio event. Via cross-correlation methods I can estimate the TDoA (or delay) between each pair of sensors.

With each pair of sensors, $i$ and $j$, we can point to the source position through the following optimization problem:

$$\delta_{ij} = v \tau_{ij}$$
$$h(P_{E},P_{i},P_{j}) = \|P_i - P_E\|_2 - \|P_j - P_E\|_2$$
$$P_E^* = {arg\,min}_{P_E}\sum\limits_{i \neq j}\left(\delta_{ij}-h(P_{E},P_{i},P_{j})\right)^2$$
where $v$ is the speed of sound in the selected medium (consider it ideally constant), $P_E^* = (x_E,y_E,z_E)$ is the source position estimate (after optimization), while $P_i = (x_i,y_i,z_i)$ and $P_j = (x_j,y_j,z_j)$ are the sensor positions.

Note the first two equations. The first one is the distance between the two sensors, calculated from the TDoA estimate. The second one is a hyperbolic surface (see image below, in a 2-D example) that also represents the distance between the sensors. One way to see this is noting (suppose perfect measurements):

$$\|P_i - P_E\|_2 - \|P_j - P_E\| = v \tau_{iE} - v \tau_{jE} = v (\tau_{iE} - \tau_{jE}) = \delta_{ij}$$

In other words, each pair of sensors generate a different hyperbolic surface. We then try find the position by minimizing the function (third equation), by varying $P_E$. For sure there will always be a $P_E$ where $\delta_{ij} = h(P_E,P_i,P_j)$, but it won't necessarily nullify the same equation for different $i$ and $j$. Therefore optimization is necessary.

In a simpler 2-D example, the hyperbolic curves above represent points where the TDoA between two sensors is constant.

Finally, my question is:
What is the best geometric arrangement of the $m$ sensors in 3-D space, so that measurements are as precise as possible? For example, a line of sensors is  probably not a good idea as the the hyperbolic curves will be parallel to each other.

This will somewhat depend on the spectrum and properties of the signals and what your dominant source of noise, distortion, jitter, non-linearietes etc. are. Ideally you can resolve every point in space with 3 non-coplanar microphones. A very simple solution is one with 4 microphones: one in the origin and three one unit step displaced in x, y, and z dimension. The arrival difference between X and origin pair gives you the x-coordinate, etc.

You can improve this method by a simple least square approach that takes into account all microphone pairings (6 total). The fit error is metric of how much noise you have in the system.

Microphone spacing depends on spectral content of the signals you are tracking. If it's fairly narrow band you probably want them to be about half a wavelength apart so you maximize the phase difference between adjacent signals.
A good robustness metric is something like dx/dt, i.e. you wiggle the arrival time a little and see how much the calculated position moves. For almost coplanar microphones that is a large number which indicates that this setup is very sensitive to noise.


Modeling the human perception of noise Add together some sinusoids with irrational or "nearly" irrational frequency ratios, and you get something pseudoperiodic. If there aren't too many, you can hear that it has only pitched components, even if you can't isolate them by ear.

If there are enough of them playing together, it will sound like steady noise.

The same goes with random transients. If they're separated in time well enough, you can hear them as distinct sounds. If they become very dense it sounds like steady noise.

Is there a model that predicts when a long-playing signal sounds like steady noise, and when it sounds merely rough or dissonant?

As with any perceptual model, the answer is not exactly cut and dry.  In my research, I’ve found that many of the hearing models I’ve encountered are based on band pass filter banks.  These banks are numerous and tightly spaced.  The more adjacent banks that are stimulated, the ‘noiseyer’ the signal would sound.

There are examples of inharmonic signal analysis that take noise into account, such as: https://pdfs.semanticscholar.org/b5bf/98425051f47b0685e75da986ee166e41d6ec.pdf.  I’ve never tried to implement such a structure, but as I understand it, you analyze the peaks in the output of a frequency transform to create a ‘deterministic’ signal.  You then subtract the deterministic signal from the original to get the ‘residual’ noise.  You would then perform further analysis on the residual.  I assume that this isn’t exactly what you were looking for in terms of a model, but I figured it may be informative.


Order analysis or order tracking Last week I've been trying to implement the Order Analysis in MATLAB in vain. I've read a lot of docs about it like this one ,  but I still can't figure it out.
This will be my last shot.

Say I have recorded 10sec audio file of a rotating machine, the sampling frequency is 96Khz. In this 10sec interval, the speed of the machine will continuously increase, say from 800rpm to 1500rpm. Is there any way to calculate the Order Analysis for this machine?
I'll be grateful for any hint!

Well my friend, if you read the following B&K tutorial about the Order Tracking, then you should encounter so called Vold-Kalman Order Filtering. Personally I don't like it (rather prefer Chirp-Fourier based methods for extracting the response of the system).

You can find plenty implementations in the internet. When you look for MATLAB, then please refer to:

Vold-Kalman order tracking code

Second generation Vold-Kalman Order Filtering.

And most certainly I suggest you to search the bottom of the prof. Jiri Tuma website: click!

For more theory please take a look at - I'm sure you will like it:

Characteristics of the Vold-Kalman Order Tracking Filter

Vold-Kalman Order Tracking Filtration

Main Principles and Limitations of the Current Order Tracking Methods

Realization of the Vold-Kalman Tracking Filter - A Least Squares problem

Advanced Vold-Kalman Filtering Order Tracking - do not be scared, it's not in chineese ;)


Why does alpha of 0.5 in detrended fluctuation analysis indicate randomness? I'm trying to get an intuitive understanding of the different coefficients in detrended fluctuation analysis (DFA).
It is used to detect fractal patterns in time series and it yields a coefficient, Alpha, which, as I understand is a Hurst coefficient mapped for different scales — the slope of which is the alpha.

What I don't understand is if alpha shows a correlation between the frequency fluctuation and the scale (i.e. the bigger is the scale, the higher are the frequency fluctuations at that scale) — why is it that a time series is scale-free when alpha is around 1 or random when alpha is around 0.5?
We have a correlation in both cases, so where do these precise values come from? Why 0.5 and not 0.3? Why 1 and not 0.9?
Maybe somebody has an explanation for this as I haven't found anything precise in scientific literature.

The application of the DFA method starts with a conversion of a bounded time series under test into an unbounded process, the elements of the unbounded process being partial sums of the bounded time series sequence.
Let we have a (discrete-time) realization of WGN of length N. It is a bounded time series. The starting conversion procedure of the DFA method transforms this WGN time series sequence into a realization of the random-walk process by the definition of the random walk process. For the RMS deviation's log-log plot of this random-walk process realization (this RMS deviation named the fluctuation in the DFA parlance), we draw a straight trend line using the least squares method. The slope of this trend line approaches 1/2 when N tends to infinity.
Hence the exact equality of DFA's alpha to 1/2 for the random process. You may have failed to find an explanation in scientific literature, because it should be read from textbooks, for example, read the Wikipedia article on DFA.
UPDATE
For an uncorrelated random walk of length N, fluctuation F(N), introduced in the DFA Wiki article, grows like a square root of N, i.e., its Hurst exponent is 1/2. Recall that the "fluctuation" of the uncorrelated time sequence such as a white noise realization does not depend on the sequence length by definition, its Hurst exponent is 0, and a white noise realization has its self-affinity at negative maximum -- the stretches of steady increase/decrease must be compensated by the stretches of the opposite direction, because the mean is zero. The alphas between zero and 1/2 correspond to an "anti-correlated" behavior from the POV of self-similarity. Likewise, the alphas between 1/2 and 1 correspond to correlation, and the alphas greater than one bear witness of the presence of non-stationary, unbounded process. The alpha=1 process is the stationary process, correlated at the maximum amount of self-affinity, it is the "1/f noise", because its PSD(f) is proportional to the inverse of frequency. The variance of an uncorrelated random walk is proportional to the square root of the series length; for the time series characterized by the Hurst exponent H, the variance and higher order moments are proportional to the series length raised to power H. This statistics results in the scaling behavior
$$
x(at) \sim a^Hx(t)
$$
where x is a random variable of the time series generating distribution, a is the scaling coefficient, H is the Hurst coefficient.
With H=1, the scaling behavior is closest to the concept of pure self-similarity, x(at) ~ ax(t), reproducing the system behavior at every time range.
These explanations do not prove the self-similarity claims, as mathematics requires it. But it would greatly help you learn the signal processing theory and master developer skills, if you compute and simulate the problems and processes. You can create an array of data and initialize it with a WGN samples in Matlab or Python; transform this array into a random walk process realization; compute trend lines with least squares and plot log-log graphs. Then, you can generate 1/f noise samples (or whatever noise spectrum you want) and plot trend lines for these distributions. Changing time scales, you will see for youself the claimed self-similarity properties. And in the process, you will master the necessary math.


FFT scaling in GNU Radio In GNU Radio Companion, I have created a very simple flow graph where I generate a 1 kHz, 2.0 volts peak-to-peak sine wave, and feed it into an FFT sink and a time domain sink.

The reference dB in the FFT sink is set to 0, and the reference scale is set to 2.0.

As such, I would expect to see 0 dB on the FFT display, but I instead see −10 dB.

Why is this? Is it a problem in my understanding, or a problem in GNURadio?

The first thing to remember is that the FFT sink does measure actual PSDs, it performs an FFT, which may be scaled however, and displays its mag-square. So when you say "2 Volts p2p", this is incorrect, as it implies certain things, such as Ohm's law, which don't work that way in DSP.

Now, for your result. First, the window you've chosen (probably Blackman-Harris) will attenuate the main lobe height by a couple of dB. If you had chosen a boxcar window (rectangular), you would have noticed the peak would be at -6 dB. Half the power! And now, things fall into place: In a real-valued FFT (as you are running), the power is distributed between the two peaks that are generated at plus and minus the sinusoid's frequency. And, sure enough, if you run the whole thing in complex mode using a rectangular window, you end up at 0 dB.
If you had chosen the QT Frequency Sink, that defaults to showing the entire Nyquist zone and you would have seen the 2nd sinusoid.

Now, is GNU Radio broken for chosen such normalizations? Not really, since they are always arbitrary. But it is always important to realize that the frequency sinks are not true power spectrum analyzers.


Need help to understanding concept of particle filter I get confuse with the basic theory of particle filter for tracking single object..
What I want to ask is:

In some paper, they explained about velocity in state vector.. Are we
randomly put number to velocity or get it from some formula?
In pdf by eiji ota, he desrcibe system model as Xn=f(Xn-1)+Wn where Wn is noise.. Why need noise in particle filter?
And can you give me simple explanation in prediction step how the particles predict at t=1? I'm still blank in there..

The following papers are good resources:

Gordon, N. J.; Salmond, D. J. and Smith, A. F. M. (1993). "Novel approach to nonlinear/non-Gaussian Bayesian state estimation". IEE Proceedings F on Radar and Signal Processing 140 (2): 107–113.
Arulampalam, M.S.; Maskell, S.; Gordon, N.; Clapp, T.; (2002). "A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking". IEEE Transactions on Signal Processing 50 (2): 174–188.
Doucet, A.; Johansen, A.M.; (December 2008). "A tutorial on particle filtering and smoothing: fifteen years later". Technical report (Department of Statistics, University of British Columbia). (Available here)

I'd recommend you read at least the first one.

A particle filter works on a model specified by 2 equations (a hidden Markov model):

1) The state equation
$x_{k+1} = f_k (x_k,w_k)$ which tells you how the state (minimal set of variables to describe a system) evolves as being driven by noise.

2) The measurement equation $y_k = h_k (x_k, v_k)$ which maps the state to what you observe.

Given the observations $y_0, \ldots, y_n$, estimate the distribution of $x_n$ is the problem the particle filter tries to solve. Once you have the approximated distribution of $x_n$ (which is the distribution of the particles), you can estimate $x_n$ (say, by averaging the particles).

The particle filter essentially starts with a bunch of samples (called particles), evolves the state by running each particle through the state equation and re-samples the particles based on the observation you see in order to make the distribution of particles consistent with the observations.


How can I create a histogram showing the RMS amplitude for each frequency bin of an audio file Using Python or Matlab, how can I create a histogram in which each bin is equal to a proportional frequency range (e.g. one octave) and the data added to each bin is a count of the number of samples that have any amplitude within that frequency range?
Here's what I have so far:
length = data.shape[0] / Fs
print(f"length = {length}s")

RMS = lambda x: np.sqrt(np.mean(x**2))

sample = np.arange(int(length))
RMS_of_sample = np.zeros(sample.shape)
for ns in sample:
    RMS_of_sample[ns] = RMS(data[ns*Fs:(ns+1)*Fs])

plt.hist(RMS_of_sample, label="Left channel")
plt.show()

But this seems to create a bin equal to the full spectral range, thus only one bar appears in the histogram. I'm also not sure that the data in the bar actually represent counts of amplitude occurrences.

If what you want is to determine the amount of energy your signal has in particular frequency bands you can perform a FFT of your signal. For that you need to make sure that you have amplitude data and not energy data. So you would want to have sound pressure Pa and not Pa^2. After having calculated a FFT with that the values have to be multiplied by $1/\sqrt{2}$ to get to the RMS-format.
Depending on your time data length and sampling time you get a specific frequency resolution of your signal in the frequency domain. With $N$ time samples and a sampling time of $\Delta t$ you get a frequency resolution of $1/(\Delta t \cdot N)$ Herz. Now if you want to know the energy contribution of a specific frequency range you can calculate the RMS-value of that frequency range using the earlier calculated linear spectra in RMS-format. According to Parseval's theorem one can use the following equation to calculate the RMS-value of a specific frequency range:

Therein the $A_k$s are the amplitudes of different frequencies in the spectra:

I used that approach to calculate the Energy of 10 Hz bins of a signal containing two sines with amplitude $2$ and Frequencies $5$ and $45$ Hz. Sampling frequency is 0.0098.
Here is the time Signal:

The resulting linear Spectra with Peak to Peak values is:

From that I calculated the RMS-Values of $\Delta 10$ Hz bins. So the first bar indicates the RMS-Value of the signal from $0 - 10$ Hz. The second from $10 - 20$ Hz and so on.

If some noise is added the spectra changes to:

The RMS-values of the $\Delta 10$ Hz bins become then:

I further suggest looking at Root Mean Square and Overall Level
Here is the Matlab Code I used to do the above (probably still containing errors though):
clear all
close all

t=linspace(0,10,1024); % time vector
y=2*sin(5*2*pi*t)+2*sin(45*2*pi*t)+randn(1,length(t)); % time signal
Fs=1/(t(2)-t(1));
Fn=Fs/2;
N=length(y);
yf=fft(y,N);
df = Fs/N; % Frequency Resolution

plot(t,y)
grid on
title('Time Signal')
ylabel('Amplitude')
xlabel(['Time Resolution: ',num2str(1/Fs),' s'])

amplH = abs(yf);
amplitudengang = fftshift(amplH/N);

x_fn = 0 : df : Fn-df;
x_fa = 0 : df : Fs-df;

figure
stem(x_fa-Fn, amplitudengang, 'b.-')
axis([-Fn Fn 0 max(amplitudengang)])
title('Two sided Spectra')
ylabel('Amplitude')
xlabel(['Frequency Resolution: ',num2str(df),' Hz'])
grid

amplitudengang=[amplitudengang(513) amplitudengang(514:end).*2];

figure
stem([0:df:(Fn-df)], amplitudengang, 'b.-')
axis([0 Fn 0 max(amplitudengang)])
title('Single Sided Spectra')
ylabel('Amplitude')
xlabel(['Frequency Resolution: ',num2str(df),'Hz'])
grid

amp10=amplitudengang(1:round((10-df)/df));
amp20=amplitudengang(round((10-df)/df):round((20-df)/df));
amp30=amplitudengang(round((20-df)/df):round((30-df)/df));
amp40=amplitudengang(round((30-df)/df):round((40-df)/df));
amp50=amplitudengang(round((40-df)/df):round((50-df)/df));

rms10=amp10./sqrt(2);
rms20=amp20./sqrt(2);
rms30=amp30./sqrt(2);
rms40=amp40./sqrt(2);
rms50=amp50./sqrt(2);

rms10=sqrt((rms10(1)/2)^2+sum(rms10(2:(end-1)).^2)+(rms10(end)/2)^2);
rms20=sqrt((rms20(1)/2)^2+sum(rms20(2:(end-1)).^2)+(rms20(end)/2)^2);
rms30=sqrt((rms30(1)/2)^2+sum(rms30(2:(end-1)).^2)+(rms30(end)/2)^2);
rms40=sqrt((rms40(1)/2)^2+sum(rms40(2:(end-1)).^2)+(rms40(end)/2)^2);
rms50=sqrt((rms50(1)/2)^2+sum(rms50(2:(end-1)).^2)+(rms50(end)/2)^2);

figure
stem([5 15 25 35 45],[rms10,rms20,rms30,rms40,rms50])
axis([0 Fn 0 max([rms10,rms20,rms30,rms40,rms50])])
title('RMS values per 10Hz')
ylabel('Amplitude')
xlabel(['Frequency in Hz'])
grid on

Edit
To reduce leakage the time data needs to be windowed. After that the resulting spectra needs to be multiplied by the matching correction factor:


frequency domain to time domain with magnitude and phase I have a signal magnitude and phase in frequency domain. I need to have it in time domain but I really have no idea how to do it. I heard something about mirroring the signal but I'm kinda new at this. I have tried to do something like this in matlab but I'm clearly missing the point. And how do I define the time vector?
z = A .* exp(i*phase);
X = ifft(z);

I've attached the signal plot
https://i.stack.imgur.com/FfbTV.jpg

Okay, I don't use MATLAB, but I think this Python example represents what you may be asking.

import numpy as np

A = 3 + 4j
k = 5
N = 128

Z = np.zeros(N, dtype='complex')
Z[k]   = A
Z[N-k] = np.conj( A )

z = np.fft.ifft(Z)

print( z )

This will generate a real valued pure tone (all the imaginary values are zero to the limit of precision) signal with 5 cycles in the frame.  The "mirroring" means that for a real signal the DFT is a conjugate mirror image, i.e. $Z[k] = Z^*[N-k]$.  It is common convention to use X, but I like Z too.


Question about the derive of filter roll-off I don't understand the follow step in my notes:

It is about the filter roll-off.

$H(w)={1\over 1+j(w/w_c)}$

$\Rightarrow$
$|H(w)|= {1\over \sqrt(1+(w/w_c)^2)}$

where $H(w) = {v_o \over v_i}$

$w_c$ = cut-off frequency

For any arbitrary rational complex function $H(\omega) = \frac{N(\omega)}{D(\omega)}$,

$$
|H(\omega)| = \frac{|N(\omega)|}{|D(\omega)|}
$$

This follows straightforwardly from the properties of the polar (i.e. magnitude/phase) representation of complex numbers. So, for your example:

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

The magnitude of an arbitrary complex number is:

$$
|a + jb| = \sqrt{a^2+b^2}
$$

So for your example:

$$
\begin{align}
|H(\omega)| &= \frac{1}{\sqrt{1^2 + \left(\frac{\omega}{\omega_c}\right)^2}} \\
&= \frac{1}{\sqrt{1 + \left(\frac{\omega}{\omega_c}\right)^2}}
\end{align}
$$