Average FFT Magnitude in bins I am analyzing sounds of daily activities recorded by a smartphone. For example walking, getting up from bed, falling, running etc.. (one at the time).
Let’s take one them. My goal is to

convert from time domain to frequency domain and then
divide it in bins of 10 Hz wide
calculate the average FT magnitude in each of these bins.

Questions:

When I use fft Matlab function, Do I have to choose NFFT>N and power of 2 (Y = fft(signal,NFFT)/N;), or can I just use Y = fft(signal); ?
In my case, is signal = detrend(signal); necessary before doing FFT?
In order to have 10 Hz wide bins, using 4096 Hamming Window is the way to do so?
Do my plots need further processing in order to achieve my goal?

[signal,fs]=audioread(filename); % fs = 44100
N = length(signal); % N = 94144

%% Amplitude Spectrum

Y = fft(signal);
Y_mag = abs(Y);
Fbins = ((0: 1/N: 1-1/N)*fs).';
figure;
subplot(4,1,1);
plot(Fbins,Y_mag);
title('FFT MAGNITUTUDE REPONSE');
xlabel('FREQUENCY (HERTZ)');
ylabel('MAGNITUTDE');

%% Single-Sided Amplitude Spectrum
NFFT = 2^nextpow2(N);
Y_NFFT = fft(signal,NFFT)/N;
f_NFFT = fs/2*linspace(0,1,NFFT/2+1);
subplot(4,1,2);
plot(f_NFFT,2*abs(Y_NFFT(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

%% Pwelch defaut values
subplot(4,1,3);
pwelch(signal,[],[],[],fs);
title(sprintf('Pwelch defaut values'));

%% Pwelch Hamming Window 4096
len = 4096; % different sizes of Hamming Window
subplot(4,1,4);
pwelch(signal,len,[],len,fs);
    title(sprintf('Pwelch Hamming Window Size :: %d',len));

Thank you very much

When I use fft Matlab function, Do I have to choose NFFT>N and power
  of 2 (Y = fft(signal,NFFT)/N;), or can I just use Y = fft(signal); ?

Not necessarily, but it is better to have it in power of 2 for faster computations.
Also NFFT>N will make your response more clearer (just for visual purpose) as the frequency resolution doesnot changes with NFFT>N. You can also look into zero padding as it is what you will achieve with NFFT>N.

As also mentioned in the document of FFT:

The execution time for fft depends on the length of the transform. The
time is fastest for powers of two and almost as fast for lengths that
have only small prime factors. The time is typically several times
slower for lengths that are prime, or which have large prime factors.

  In my case, is signal = detrend(signal); necessary before doing FFT?

It depends on whether you want to measure the frequency component of 0-10 Hz range (specifically from your description)
If DC component is of no need then you can perform signal = detrend(signal);

  In order to have 10 Hz wide bins, using 4096 Hamming Window is the way
  to do so?

NO, first of all windowing is used to decrease the side loobes which occurs due to spectral leakage. There is a effect of response broadening due to side level suppression but I think this is not what you needed. Also do remember that windowing will be applied on raw time signal before zero padding.
You have to play with number of samples to have your desired resolution

  Do my plots need further processing in order to achieve my goal?

I cannot guess this, but might be from the above answers you can estimate what has to be done in order to achieve your specific goal.

Hope that this answer might help.


How to apply neural network parallelly on each part of image I have an image which is divided into four equal size square blocks. I want to apply neural network for denoising. Usually, I apply on the whole image. But i was thinking that is it possible to divide the image into equal parts and apply NN on each part parallelly.

Anyone, please help me how can I do that

The question has multi-fold answers:

First, you can simply crop the image 4 times and then pass them to the network. To parallelize you can use OpenMP. Of course the size of the input layer should match the size of the cropped patch. And just make sure your network doesn't use global variables.

Next, I don't think it's a good idea to divide the image into parts. Another approach would be to adjust the network connections such that different parts are not connected in the first layers, but accommodate the information in the deeper ones. In short, the desired way should be encoded into the network, rather than some ad-hoc external procedures. This way your network can be trained end-to-end.


Abs(poles) < 1 by what margin for a stable filter? I checked the literature for recent algorithms used to design a digital filter that is a minimax approximation of a desired frequency response. All the articles I found work out examples where all the poles have magnitude less than 0.92, or less that 0.89, etc. I havn't seen a published example with a pole having magnitude 0.95 and certainly not 0.9995. If a filter is implemented with double precision machine arithmetic, how close to the unit circle can a pole get to get a useful filter? Would it be a bad idea to have a pole with magnitude 0.95?

It would most likely depend upon the order of the filter, but having a pole at position $z_p$ where $|z_p| = 0.95$ shouldn't pose a stability problem if you're using double-precision floating-point arithmetic with a filter of reasonable length.

Since I don't know anything about your specific application, one other effect of the pole magnitude that might be relevant is the resulting filter's impulse response. Recall that poles in a system's transfer function correspond to decaying exponential terms in the system's impulse response. As Wikipedia notes:

$$
a^nu[n] \Leftrightarrow \frac{1}{1-az^{-1}}
$$

where $u[n]$ is the discrete unit step function, used here to express that the impulse response is causal ($0\ \forall\ n < 0$). Therefore, if your filter has a pole at $z = a$, then there will be a corresponding $a^n$ term in the filter's impulse response.

For small $|a|$, this term will decay out in a relatively small number of samples. As $|a| \rightarrow 1$, the amount of time (measured in samples) required for the exponential function to decay increases.
When you reach the "critically stable" point of $|a| = 1$, the exponential never decays and the system impulse response does not decay to zero.
If $|a| > 1$ (i.e. the pole lies outside the unit circle), then the exponential function diverges and the impulse response blows up to infinity; this is why a discrete-time system is not BIBO stable if it contains poles outside the unit circle.

Given the above, the other concern that you might have is the overall effective time duration of the filter's impulse response. Although, as its name would suggest, an IIR filter theoretically has an impulse response of infinite length, in practice the response will usually decay to a negligible level after a finite time period. If your application is sensitive to this characteristic, then it makes sense to choose pole locations that are further from the unit circle. There will be corresponding tradeoffs in frequency-domain characteristics, as placing poles near the unit circle can help to make transition regions sharper and more narrow, as is often desirable.


Autocorrelation for the sum of a deterministic and stochastic signal so I'm working on a problem where I need to find a PSD for the sum of a deterministic and stochastic function. So assume $x(t)$ being this sum:

$x(t) = d(t)+s(t)$

where $d(t)$ is the deterministic function and $s(t)$ is the stochastic function. How do I write down the autocorrelation $R(\tau )$ and the PSD for $x(t)$ that we calculate from Fourier transforming $R(\tau )$?

The random process $x(t)$ doesn't have a power spectral density (PSD) in the conventional sense because it is non-stationary, due to the presence of the deterministic time-dependent function $d(t)$.

You can compute the auto-correlation function of $x(t)$ as follows:

$$\begin{align}R_x(t_1,t_2)&=E\{x(t_1)x^*(t_2)\}\\&=E\{[d(t_1)+s(t_1)][d^*(t_2)+s^*(t_2)]\}\\
&=d(t_1)d^*(t_2)+d(t_1)E\{s^*(t_2)\}+d^*(t_2)E\{s(t_1)\}+E\{s(t_1)s^*(t_2)\}\\&=d(t_1)d^*(t_2)+d(t_1)\mu_s^*(t_2)+d^*(t_2)\mu_s(t_1)+R_s(t_1,t_2)
\end{align}\tag{1}$$

where $\mu_s(t)$ is the (generally time-dependent) mean of the random process $s(t)$, and $R_s(t_1,t_2)$ is its auto-correlation function.

If you assume that $s(t)$ is wide-sense stationary (WSS), and if you further assume (for ease of notation) that $d(t)$ and $s(t)$ are real-valued, then Eq. $(1)$ simplifies to

$$R_x(t,\tau)=d(t)d(t+\tau)+[d(t)+d(t+\tau)]\mu_s+R_s(\tau)\tag{2}$$

where I've used the usual substitution $t=t_1$ and $\tau=t_2-t_1$.
Note that $(2)$ generally still depends on the absolute time $t$, and, consequently, $x(t)$ is non-stationary and has no PSD.

As already pointed out in a comment, only if $d(t)=d$ is constant (and $s(t)$ is WSS) does the random process $x(t)$ become WSS with auto-correlation function

$$R_x(\tau)=d^2+2d\mu_s+R_s(\tau)\tag{3}$$

and its PSD becomes (by taking the Fourier transform of $(3)$)

$$P_x(\omega)=2\pi d(d+2\mu_s)\delta(\omega)+P_s(\omega)\tag{4}$$

where $\delta(\omega)$ is the Dirac delta impulse.

EDIT: As mentioned in a comment by Dilip Sarwate, $P_s(\omega)$ also generally includes a Dirac impulse if the mean $\mu_s$ is non-zero. If you define a zero-mean random process $\tilde{s}(t)$ by

$$s(t)=\tilde{s}(t)+\mu_s\tag{5}$$

then $P_s(\omega)$ can be written as

$$P_s(\omega)=P_{\tilde{s}}(\omega)+2\pi\mu_s^2\tag{6}$$

Using $(6)$, Eq. $(4)$ becomes

$$P_x(\omega)=2\pi (d+\mu_s)^2\delta(\omega)+P_{\tilde{s}}(\omega)\tag{7}$$

which explicitly shows all Dirac impulses in the PSD of $x(t)$.


Can the Cramér-Rao Lower Bound be constant? In many text books and research papers, I have seen that for the plot of Cramér-Rao Lower Bound vs Signal to Noise Ratio, the CRLB decreases with increasing SNR.

Is it possible that the CRLB remains constant for all SNR?
If so, what is the implication?

Prologue

The CRLB of an unbiased estimator $\hat \theta$ for the paramater $\theta$ is

$$
\text{var}(\hat \theta) \ge - \frac1 {\text E_x\left\{ \frac{\partial^2 (\ln f(x;\theta))}{\partial \theta^2} \right\}}
$$

with $f(x;\theta)$ being the likelihood function of $\theta$ under the observation $x$. It is the inverse of the so-called Fisher information.

Answer

  Is it possible that the CRLB remains constant for all SNR?

Yes. For that, the expectation (over all $x$, hence the index)

$$\text E_x\left\{ \frac{\partial^2 (\ln f(x;\theta))}{\partial \theta^2} \right\} $$

just needs to be constant.

The "easiest" explanation would be if noise just didn't have anything to do with the likelihood – ie. noise is just not actually observed by the estimator.

But that's kind of least exciting result – "well, if I add a lot of noise to something and then observe something completely unrelated, the noise has no effect" isn't that surprising.

A lot more complicated, is the option that the Fisher information is actually constant, although $f(x;\theta)$ actually depends on the noise.

The expectation operator $E_x$ is definable as

$$
\text E_x \{g(x)\} = \frac1{\mu (\xi)} \int\limits_{x\in\xi} g(x)\,d\mu\text,
$$

i.e. as the integral over all possible $x$ divided by the measure $\mu$ of the set $\xi$ from which $x$ can be chosen.

Assuming smoothness of $x$ and compactness of $\xi\subset \mathbb R^1$ (just to make all the derivates easier to handle and to avoid cascading integrals), this specializes to:

$$\begin{align}
\text E_x \{g(x)\} &= \frac1{x_{max}-x_{min}} \int\limits_{x = x_{min}}^{x_{max}} g(x)\,dx\\
\text E_x\left\{ \frac{\partial^2}{\partial \theta^2}{ (\ln f(x;\theta))} \right\}  \cdot\left({x_{max}-x_{min}}\right) &= \int\limits_{x = x_{min}}^{x_{max}} \frac{\partial^2}{\partial \theta^2}{ (\ln f(x;\theta))} \,dx\\
&= \left[ \frac{\partial}{\partial \theta}{ (\ln f(x;\theta))}  \right]_{x = x_{min}}^{x_{max}}\\
&= \left[ \frac1{f(x;\theta)}  \frac{\partial f(x;\theta)}{\partial \theta}  \right]_{x = x_{min}}^{x_{max}}\\
&= \frac1{f(x_{max};\theta)}  \frac{\partial f(x;\theta)}{\partial \theta}_{x=x_{max}}
- \frac1{f(x_{min};\theta)}  \frac{\partial f(x;\theta)}{\partial \theta}_{x=x_{min}}\\
&= c,&c\text{ const.}\\
\end{align}$$


Automatic Cropping of Arbitrary Shapes I have an arbitrary shape defined by a binary mask (gray = shape, black = background).

I would like to find a largest possible rectangle containing only gray pixels (such rectangle is pictured in yellow):

The shape is always "one piece" but it is not necessarily convex (not all point pairs on the shape's boundary can be connected by a straight line going through the shape).

Sometimes many of such "maximum rectangles" exist and then further constrains can be introduced, such as:

Taking the rectangle with its center nearest to shape's center of mass (or center of image)
Taking rectangle with aspect ratio nearest to a predefined ratio (i.e. 4:3)

My first thought about the algorithm is the following:

Compute distance transform of the shape and find its center of mass
Grow square area while it contains only shape's pixels
Grow the rectangle (originally a square) in width or height while it contains only shape's pixels.

However, I think such algorithm would be slow and would not lead to optimal solution.

Any suggestions?

There is a code on Matlab Fileexchange that is relevant to your problem:
http://www.mathworks.com/matlabcentral/fileexchange/28155-inscribedrectangle/content/html/Inscribed_Rectangle_demo.html

Update

I wrote this tutorial article on computing largest inscribed rectangles based on the above link from Atul Ingle.

The algorithm first searches for largest squares on a binary mask. This is done using simple dynamic programming algorithm. Each new pixel is updated using the three neighbors already known:

squares[x,y] = min(squares[x+1,y], squares[x,y+1], squares[x+1,y+1]) + 1

The sample binary mask and computed map look like this:

Taking maximum in the map reveals the largest inscribed square:

The rectangle searching algorithm than scans the mask two more times looking for two classes of rectangles:

width greater than square's size (and height possibly smaller)
height greater than square's size (and width possibly smaller)

Both classes are bounded by the largest squares since no rectangle at a given point can have both dimensions greater than the inscribed square (though one dimension can be larger).

One have to choose some metric for rectangle sizes, like area, circumference or weighted sum of dimensions.

Here is the resulting map for rectangles:

It is convenient to store position and size of the best rectangle found so far in a variable instead of building maps and then looking for maxima.

The practical application of this algorithm is cropping non-rectangular images. I have used this algorithm in my image stitching library SharpStitch:


FFT equivalent for generalized unitary transforms The DFT has the FFT, Hadamard transform has the Fast Hadamard Transform and so do a number of other unitary transforms (operators). Is there or has there been an attempt at creating FFT style algorithms for generalized unitary transforms?

It's all about structure. One early paper on this is  A Unified Treatment of Discrete Fast Unitary Transforms, 1977:

  A set of recursive rules which generate unitary transforms with a
  fast algorithm (FUT) are presented. For each rule, simple relations
  give the number of elementary operations required by the fast
  algorithm. The common Fourier, Walsh-Hadamard (W-H), Haar, and Slant
  transforms are expressed with these rules. The framework developed
  allows the introduction of generalized transforms which include all
  common.transforms in a large class of "identical computation
  transforms". A systematic and unified view is provided for unitary
  transforms which have appeared in the literature. This approach leads
  to a number of new transforms of potential interest. Generalization to
  complex and multidimensional unitary transforms is considered and some
  structural relations between transforms are established.

Among the most common, the discrete sine (DST), cosine (DCT), Hartley, wavelet transforms and many others (Walsh, Hadamard, Paley or Waleymard, triangle, jacket, slant, Hermite) have fast counterparts.

One global initiative towards a systematic construction was termed Algebraic Signal Processing Theory. Let me mention two papers laying some foundation using basic algebraic structures:

Algebraic Signal Processing Theory: Foundation and 1-D Time, 2008

  This paper introduces a general and axiomatic approach to linear
  signal processing (SP) that we refer to as the algebraic signal
  processing theory (ASP). Basic to ASP is the linear signal model
  defined as a triple ($\mathcal{A}$, $\mathcal{M}$, $\Phi$) where
  familiar concepts like the filter space and the signal space are cast
  as an algebra $\mathcal{A}$ and a module $\mathcal{M}$, respectively.
  The mapping $\Phi$ generalizes the concept of $z$-transform to
  bijective linear mappings from a vector space of signal samples into
  the module $\mathcal{M}$. Common concepts like filtering, spectrum, or
  Fourier transform have their equivalent counterparts in ASP. Once
  these concepts and their properties are defined and understood in the
  context of ASP, they remain true and apply to specific instantiations
  of the ASP signal model. For example, to develop signal processing
  theories for infinite and finite discrete time signals, for infinite
  or finite discrete space signals, or for multidimensional signals, we
  need only to instantiate the ASP signal model to a signal model that
  makes sense for that specific class of signals. Filtering, spectrum,
  Fourier transform, and other notions follow then from the
  corresponding ASP concepts. Similarly, common assumptions in SP
  translate into requirements on the ASP signal model. For example,
  shift-invariance is equivalent to $\mathcal{A}$ being commutative. For
  finite (duration) signals shift invariance then restricts
  $\mathcal{A}$ to polynomial algebras. We explain how to design signal
  models from the specification of a special filter, the shift. The
  paper illustrates the general ASP theory with the standard time shift,
  presenting a unique signal model for infinite time and several signal
  models for finite time. The latter models illustrate the role played
  by boundary conditions and recover the discrete Fourier transform
  (DFT) and its variants as associated Fourier transforms. Finally, ASP
  provides a systematic methodology to derive fast algorithms for linear
  transforms. This topic and the application of ASP to space dependent
  signals and to multidimensional signals are pursued in companion
  papers.

Algebraic Signal Processing Theory: Cooley-Tukey Type Algorithms for DCTs and DSTs, 2008

  This paper presents a systematic methodology to derive and classify
  fast algorithms for linear transforms. The approach is based on the
  algebraic signal processing theory. This means that the algorithms are
  not derived by manipulating the entries of transform matrices, but by
  a stepwise decomposition of the associated signal models, or
  polynomial algebras. This decomposition is based on two generic
  methods or algebraic principles that generalize the well-known
  Cooley-Tukey FFT and make the algorithms' derivations concise and
  transparent. Application to the 16 discrete cosine and sine transforms
  yields a large class of fast general radix algorithms, many of which
  have not been found before.

Aside, when those transforms are considered in the filter bank framework, other fast versions exist, based on polyphase decomposition or lifting transforms.

Additional references: I do not know of a full book on  Algebraic Signal Processing Theory. The above link has a handful of references. For book-style, you have PhD theses:

A. Sandryhaila, 2010, Algebraic Signal Processing: Modeling and Subband Analysis
M. Püschel, 2008, DFT and FFT: An Algebraic View,  Computer Algebra Handbook, Foundations, Applications, Systems, Eds. J. Grabmeier, E. Kaltofen, V. Weispfenning

On filter banks:

L. Liu,  On Filter Bank and Transform Design with the Lifting Scheme
K. Soman et al., 2020, Insight into Wavelets: from Theory to Practice


How to measure the time dependent correlation of two signals I have two signals (time series) shown in the plot below. Just by looking at the figure, we can see that the two main peaks of both signals are very closely aligned (correlated), however the red signal has additional features elsewhere which don't match the blue curve.

I am looking for some function which can tell me how correlated these two signals are at a given time. IE: at t = -18 and t = +10 where the peaks are located, the correlation should be very close to -1, while at t = -40, the correlation is probably close to 0.

Specifically, I want a robust way to tell that if I see a peak in the blue curve, then there is a corresponding (negative) peak in the red curve.

I've looked into Cross correlation, however this seems suitable for time series which simply have a constant time lag with respect to each other. In my problem the correlation between my two signals is changing with time. I've thought about dividing up the signal into small bins and computing the correlation of each bin, but I was hoping there is a better way.

Your signals look like they have about the same ordinal scale, and little lag. Hence, you can perform local correlations: take a segment of both signals centered at a time location, you can normalize them in amplitude  and dispersion (divise by mean and standard deviation, or median and median absolute deviation). Apply a smoothing window: if you do not do that, remember you are already applying a rectangular window. Then do a cross-correlation, and repeat over time.

Such a process creates a 2D image from two 1D signals. It  could be called windowed cross-correlation (see also Windowed cross correlation in Matlab), or local cross-correlation. From the latter:

  We address the problem of capturing and tracking local correlations
  among time evolving time series. Our approach is based on comparing
  the local auto-covariance matrices (via their spectral decompositions)
  of each series and generalizes the notion of linear cross-correlation

If you are unsure about scales, or want of more generic setting, some perform cross-correlation in a multiscale wavelet domain. You can check the concept of wavelet coherence (Matlab), illustrated by Torrence et al. Interdecadal Changes in the ENSO–Monsoon System
, with Matlab software


Influence of image size to edge detection in OpenCV I need to scan an image in OpenCV and obtain the lines. For this I use the canny edge detector and then the propabilistic Hough transform, with which I actually get the lines.
The problem now is that I get many short lines instead of a few continuous ones. The image I am analyzing has the size 1024x768 and I am wondering if this could be an error source. Theoretically making the image smaller should lead to fewer lines. Might that help me in my task? Or are there any other factors which I should take into account when extracting lines?

You should play around with the many parameters of the Canny() and HoughLinesP() functions before resorting to changing the image size.

The prototype for Canny() is:

void Canny(InputArray image, OutputArray edges,
           double threshold1, double threshold2,
           int apertureSize=3, bool L2gradient=false )

You should play with the threshold1, threshold2, aperturSize and L2gradient parameters until you get an output from Canny() that accurately marks the pixel involved in the lines in your image.  The threshold2 parameter controls the intensity above which a gradient value will be considered definitely part of an edge.  Setting it higher will give you fewer edges.  The threshold1 parameter controls the intensity above which a pixel will be considered to continue an edge from a neighboring pixel.  Setting it lower will give you longer edges.  The ratio threshold2/threshold1 should be at least 2.0 and can be as large as 3.0, and maybe larger.  Setting L2gradient to true will make the algorithm slightly slower but will use a more accurate gradient calculation from the first derivatives.  Setting apertureSize to 5 or 7 might smooth out some noise in the first derivative calculation.

Only once you have Canny() accurately finding the edges in your image should you move on to running HoughLinesP().  HoughLinesP() also has a bunch of knobs that you need to tune to get decent results.

void HoughLinesP(InputArray image, OutputArray lines,
                 double rho, double theta,
                 int threshold,
                 double minLineLength=0,
                 double maxLineGap=0 )

rho and theta can be set larger and smaller depending on how many lengths and angles you want to allow.  Setting rho and theta larger will give you fewer line segments, but they will be less accurate.  threshold can be increased to reduce the number of line segments that you get.  minLineLength can be increased to reject short line segments and maxLineGap can be increased to make it easier for short lines to get connected together.

Setting Canny()'s apertureSize parameter larger is somewhat like working with a smaller image, because a larger aperture is doing more blurring/smoothing of the details of the image.  Setting HoughLinesP()'s rho and theta parameters larger is somewhat like working with a smaller image because you are lumping together small differences in angle or length that would be indistinguishable in the smaller image.

The OpenCV documentation for HoughLinesP() has an example with:

Canny( src, dst, 50, 200, 3 );
...
HoughLinesP( dst, lines, 1, CV_PI/180, 80, 30, 10 );


What unit is the z axis of a power spectrogram of audio? the question might seem dumb but since I tried it out and got values well above 2000 (watts?) it got me thinking... Check this out:

winLen = 10*1e-3;
fftOverlap = 0.5;
[x, Fs] = audioread('speech_5.wav');
x = x(:,1);
hopSize = winLen * Fs;
nfft = hopSize / fftOverlap;
overlap = nfft - hopSize;
w = hann(nfft);
[S, f, t] = spectrogram(x, w, overlap, nfft, Fs);
S = (abs(S)).^2 + eps; %make it a power spectrogram

If you look into S (the spectrogram) now, you see a lot of values above 1000, in my case (this is clean speech) even 6000. It got me thinking, what kind of unit this is? Intuition says its watts, but 6000 watts of power in a certain frequency in a 10ms time frame seems very unlikely to me... I mean my hair dryer has 2000 watts... Can someone please help me out here?

Additionally, if you see anything in my code that could be improved, please share! I will be more than thankful for any input!

The short answer is that without a calibrated data acquisition chain from the original sound as it was picked up from at a microphone and converted to a digital file, it doesn’t have any units.   Its just numbers that fit into 16 bits (or 24 or 8) of dynamic range.

The number you have is probably proportional to something in nature if the original sound was natural but something coming out of Garage Band is synthetic.


Decimation FFT Result In my Signal Processing class we just learned decimation however after messing around in Matlab trying an example of decimation, I just can't understand what's happening.

I get the introducing twi zeros between samples and the after filtering steps however I don't get how the original signal comes to look like that in the decimated step(EDIT: M=3).

These are the frequency spectrums and the OP is seeing the effects of decimation (upper figure) and interpolation (lower figure). Note the frequency plots are the magnitude of the DFT result, so the left half of the spectrum represents the positive frequencies from $0$ (DC) to $F_s/2$ (where $F_s$ is the sampling rate) and the left half represents the negative frequencies from $-F_s/2$ to $0$ (nearly so, the bin at "0" only appears once as the first bin).  We start with the spectrum in the upper left, which gets decimated to be the spectrum in the upper right. This decimated waveform is then interpolated as shown in the lower two plots recovering the original signal in the spectrum in the lower right. Imperfections in the interpolation filter result in the degraded noise floor in the final result (spectral leakage), suggesting the filter used was far from ideal for interpolation (especially given this isn't a log magnitude plot!).
Decimation by D is filtering the signal first to remove aliases at higher frequencies and then keeping only every D'th sample to result in a lower sampled waveform. Given the first (upper left) spectrum has no higher frequencies or was already filtered, and the result on the right is what would happen if only every 3rd sample was taken (2 samples removed between every 3rd sample.). Decimation without the filter operation is simply "Downsampling".
Interpolation by I is first inserting $I-1$ zeros in between every sample which creates 2 image replicas of the spectrum in between the original images (which is what we see happen in the plot on the lower left), which must then be filtered to complete the full interpolation operation (the filter if done right will grow those zeros to the ideal interpolated sample). Interpolation without the filter operation is simply "Upsampling".
I explain this further at this link along with how to design a better interpolation filter: The condition in order to not loose information after upsampling and downsampling
To understand these multi-rate concepts and how the aliasing occurs, it is very helpful to thoroughly understand how the frequency spectrum for sampled waveforms is periodic as explained in this post showing the A/D process and the resulting periodicity in frequency, starting with knowing that convolution in one domain is multiplication in the other:  Higher order harmonics during sampling


Doppler effect on phase I am struggling to find or derive the relationship between Doppler frequency and the associated phase shift.  Let's say I have a signal transmitted at a given frequency $f_c$.  I then down convert from $f_c$ to baseband at the receiver, but due to speed between the transmitter and receiver, I receive a frequency shifted by a Doppler frequency $\Delta f$.

If I wish to generate a matching signal to that being received, I can use

$$
\cos(2 \pi \Delta f t + \phi) + i \sin(2 \pi \Delta f t + \phi)
$$

The part I'm struggling to get to grips with is how does $\phi$ relate to $\Delta f$?  Is it simply the integral of $\Delta f$ with respect to time?

$$ \phi = \int_0^t \Delta f dt $$

You are correct that the total accumulated phase offset due to Doppler could be calculated by integrating the Doppler shift $\Delta f(t)$ over time. In addition to Doppler shift, however, you would also need to take into account any frequency difference between the transmitter and receiver's reference oscillators.

However, I would question the value of trying to do so. In most communications systems, there will be a random, unknown phase offset between the transmitter and receiver oscillators. This makes absolute phase measurements meaningless unless you have some reference. In real systems, phase synchronization is either achieved by exploiting some known characteristic of the signal (e.g. the known constellation of a PSK signal), or in some rare cases, via the explicit synchronization of the transmitter and receiver local oscillators.


Using DTW on the MFCC's to match a spoken word to a template I am trying to code a simple algorithm in MATLAB that would be used to detect a single word.  The goal would be to have a user record the word one time to act as a template, then make the user repeat the same word to try and detect it. Searching on this website has answered many of my questions, but I still have some interrogations.

So far my algorithm calculates the MFCC's in the following way:

Extract the data from a .wav file (recorded at $8000\textrm{ Hz}$)
Build frames of $64\textrm{ ms}$ with a $50\%$ overlap, applying a Hann window
Calculate the FFT of all these frames
Calculate the Power Spectrum using the FFT results ($P(k) = \lvert X(k)\rvert^2$)
For each frame, extract 26 coefficients by multiplying the Power Spectrums with the 26 Mel filters and summing the results for a given Mel filter
Calculate the $\log_{10}$ of the coefficients
Apply a dct on a frame's 26 coefficients

My goal was to use the MFCC's of the template as well as the MFCC's of the repeated work and compare them using a DWT algorithm. My DWT algorithm is already programmed and functional.

However, my algorithm does not work very well and it seems like certain parameters affect the results quite a lot. Here are my questions:

Are these steps enough to detect a spoken word?
Would it be better to use a large number of fixed templates instead of having the user pre-record a template? Which method should result in better recognition performance?
If the person repeats the word using a different distance from the microphone, the increased/decreased values of the MFCC's are enough to mess up the DTW results. Is there a smart way to normalize the MFCCs to try and cancel the effect of the microphone distance?
Some websites recommend using only the 13 MFCC's with the smallest values.

Why is that?
Also, are they talking about the smallest magnitudes, or the smallest values?
Assuming 13 MFCC's are very big negative numbers, while 13 other MFCC's are small positive numbers, which set would I keep?

EDIT: Also, the first coefficient of each frame is always much bigger than the other coefficients. I would say it's magnitude is bigger my a factor of 100. Obviously, when I caculate the DTW using Euclidian distance, I'd say this coefficient is the only relevant one since it is so much bigger than all the other values. Should it be discarded?

Yes, this should be enough for a basic isolated word recognition system. Probably not something for a commercial product, but good enough for a university project or demo...
It would be better to ask the user to record a word and match against this, rather than attempt to match against a large database of utterances of the same word by different speakers. The larger the database the more likely you'd find an utterance matching the user's voice, but matching against a large database would be computationally expensive!
The first MFCC captures the signal energy - you might try to remove it when doing the matching (your feature vector will thus contain MFCC #2 to 13) - this will make recognition more robust to changes in amplitude.
You might have misunderstood something. Do you have a link to such a website?

A general proporty of MFCC is that as the coefficient index increases, the standard deviation of this coefficient decreases (the DCT can be thought of as an approximation of a PCA) - so the lower coefficients will have greater magnitude. You can compensate for that by estimating the standard deviation of each MFCC coefficient on a training set, once and for all, and dividing each coefficient by its standard deviation (or a power of it < 1.0) prior to the DTW step.


filling holes of a 3D image in Matlab I want to use morphology operator for filling holes of a 3d image. The 3d image consists of 116 2-D image slices of the same size (512*512 pixels). I tried using imfill in Matlab for the 3d image but the result is not I want. One of the slices of the 3d input image(a) and the corresponding slice from the 3d output image(b) result of using imfill  is shown here . The small holes are filled and other holes are not filled.

When I tried using imfill slice by slice the result is correct. Is there any way for filling holes in 3d images? How can I do this?

imfill would work with images and not volumes. However Matlab has many basic tools to fill the holes, that is to interpolate missing data points. What you are after is interpolation and Matlab has a function called interp3.


How to decompose raw acceleration signals into gravity and dynamic acceleration? I have a data from $3$-axis $50$  Hz accelerometer, acceleration over each of $3$-axis: $acc_X,\  acc_y,\  acc_z$.

I want to decompose acceleration into gravity and movement components.

First question:

Should I introduce some notation of magnitude, for example $mag = \sqrt{acc_x^2 + acc_y^2 + acc_z^2}$ and then apply all transformations to the magnitude rather than component-wisely?

The second question is the decomposition itself.

One solution I've found is to apply lowpass Butterworth filter. However, I don't know how to interpret ButterWorth filter return values, namely numpy arrays $a$ and $b$, in terms of gravity and movement component.

So, how can I implement decomposition using the Butterworth filter?

Assuming the platform does not rotate, then the gravity component will always point in the same direction.

Let's say $acc_x = A_x + g$ assuming that the $x$ axis is pointing upward, and $A_x$ is the acceleration due to motion relative to Earth's surface in the vertical direction.

So it's just a matter of estimating $g$ and subtracting: $A_x = acc_x - g$.

To estimate $g$, if you can place your platform in a motionless state, then measure and average several (hundreds) samples.

If you can not set your platform in an motionless state, but you know that the platform will only vibrate and never drift away, then you can still estimate and cancel out $g$ during sampling, since $g$ will be the 0-frequency component of the acceleration.

You can use averaging (moving window) or, as you suggest, any low-pass filter (such as Butterworth). E.g. if your platform vibrates at frequencies greater than 10Hz, then design a filter with cut-off frequency of 1Hz. And discard the first 10 seconds of sampling as the filter output converges to $g$.

Make sure that the cut-off frequency of your filter is way below the smallest frequency of your motion, so that you don't filter out your signal.

You can use this method independently on the 3 axes, so you also cancel out any DC offset of your sensors, or any misalignment of your platform.


Matrices of complex numbers multiplication I'm trying to implement the multiplication of two matrices something like this picture in c langage.
I want to read the numbers from a text file of x and store it later in an array the code that i tried is below
So the problem that I have is when i use the  static method Complex x[4] = { { 1,0 },{ 0,-1 },{ 1,1 },{ 4,-1 } }; it's work but when i want to read it from input.txt Complex x[4] it does not work what i missed.
typedef struct Complex
{
    double real;
    double image;
}Complex;

int main()
{
   Complex x[4] = { { 1,0 },{ 0,-1 },{ 1,1 },{ 4,-1 } };
   Complex y[4] = { { 0,1 },{ 1,-1 },{ 2,-3 },{ 4,0 } };
   FILE *input;
   input = fopen("input.txt", "r");
   if (sensorInput == NULL ) {
        printf(stderr, "Cannot open file\n");
        return 1;
    }
    else {
        for (i = 0; i<4; ++i)
        {
                fscanf(sensorInput, "%lf + j*%lf\t", &x[i]);
                printf("\n");
        }
    }

}

input.txt
1,0
0,-1
1,1
4,-1

The line
                fscanf(sensorInput, "%lf + j*%lf\t", &x[i]);

should be
                fscanf(sensorInput, "%lf,%lf\t", &x[i].real, &x[i].imag);

Notice that C++ provides complex type, they support the common complex operations, and provide an io, and they can read complex numbers writen as (1,3) directly file_in >> x[i], once you have a istream
Note that this type of question you should ask in stackoverflow instead.


How to remove a costant group delay? I have to subtract, from a signal x[n], its version filtered thorugh an LPF FIR filter, like this:

xf=filter(h_lpf,1,x);
y=x-xf;

How can I remove the group delay introduced by the filter h_lpf into the signal xf?

I know it is equal to $\tau_{grp}=\frac{N-1}{2}$, where N is the number of taps.

In my case N is odd, so $\tau_{grp}$ is an integer.

In this page:

http://it.mathworks.com/help/signal/examples/practical-introduction-to-digital-filtering.html#zmw57dd0e2544

in the section (Compensating for Delay Introduced by Filtering) there's a method, which I here adapt for my case: (note: "x" is the input signal, and D is the group delay):

xf = filter(h_lpf,1,[x zeros(1,D)]);      % Append D zeros to the input data
xf = xf(D+1:end);                         % Shift data to compensate for delay
y=x-xf;

Could anyone explain it to me? Could I use it for my problem?

If the low pass filter has a linear phase response, then the group delay it introduces is indeed constant, and it equals $\tau_g=(N-1)/2$ samples, with $N$ being the number of taps. If $N$ is odd, then the group delay is an integer number of samples, and in this case you just have to shift the output signal $\tau_g$ samples to the left as compared to the original signal:

$$y[n]=x[n]-x_f[n+\tau],\quad \tau\in\mathbb{N}$$

This is exactly what the Matlab code in your question does. Don't forget to make sure that the number of filter taps $N$ is odd!

Alternatively, you could delay the original signal $x[n]$ by the filter's group delay before subtracting from it the filtered signal.


Transmission frquency in MIMO-OFDM system What I understand in MIMO system, the user data is broken down into several chunks and transmitted over the antennas (say 3 Tx antennas) by using same frequency Fc. In order to use the interference as an advantage in this case.

And what I know of OFDM, the large BW is divided into several smaller ones, with each having a specific frequency. For example, 2.4GHz BW is divided into 14 channels, each channel has a different central frequency (say the 6th channel has a frequency of 2.437 MHz).

Now the question is while we combine this 2 method of MIMO-OFDM for transmission, does that transmit in same frequency from each antenna? Or at different frequencies?

For a given subcarrier, QAM symbols are transmitted on several antennas (MIMO).
And an antenna can transmit several subcarriers (OFDM). One can imagine that QAM symbols are modulated on a 3-dimensional grid frequency-space-time.


Amplitude in frequency domain of the same signal measured at different bandwidth (f_max - f_min) and different FFT lines Fourier series in complex form

Fourier series in trigonometry form

Here, I'm measuring the same signal (an impulse force - or at least as close as possible to an impulse, represented by a knocking force). The work I'm doing is parameter study (i.e.: studying the effect of different bandwidths and numbers of FFT lines).
This is on the same number of FFT lines but different bandwidth

And this is the same bandwidth and different FFT lines

The main conclusion I can draw is that changing one parameter (while keeping the other is the same) would change the amplitude of the frequency spectrum. I can explain this via math above (smaller frequency resolution leads to longer measuring time, and thus lower Fourier coefficients).
But what does it mean in a practical and realistic manner? As in for actual sampling process, how would changing the parameters (measuring time, sampling points...) affect the frequency spectrum of the same signal?
Edit: The smaller frequency resolution would eventually lead to lower amplitude in frequency spectrum is considered as analytical math by my instructor, while the actual system works with numerical math. So a more "physical" (or "what actually happens") answer is needed.

Are you using the correct FFT form for spectral analysis?  You need to divide by N when mapping time->frequency domain with an N-point FFT for a power spectrum. If you want power spectral density you have to factor in the bin bandwidth and window effective noise bandwidths.
The amplitude of a spectral line represents the power of a single signal frequency, power spectral density reflects power per unit bandwidth, which is different.


Spectral Analysis of a Time Series with Missing Data Points I use a PC to record time series of some physical property. The problem is that, for some reason, I did not record the time series as a whole, rather I record first segment, then second, third, etc. Each segment is 500,000 points, representing ~30 min. However, they are not strictly consecutive: between each segment and the other there will be a 'lag', or missing points, equal to a few seconds.

Since the number of missed points is too small compared to the total number of samples for each segment, can I simply concatenate the segments and treat them as a single time series? At the end I will be interested in doing FFTs, histograms, etc.

Given $ \left\{ x \left[ n \right] \right\}_{n \in M} $ where $ M $ is the set of indices given for the samples of $ x \left[ n \right] $.
The trivial solution (Which it would be great to have a faster more efficient solution is what I'm looking for) would be:
$$ \arg \min_{y} \frac{1}{2} \left\| \hat{F}^{T} y - x \right\|_{2}^{2} $$
Where $ \hat{F} $ is formed by subset of columns of the DFT Matrix $ F $ matching the given indices of the samples, $ x $ is the vector of the given samples and $ y $ is the vector of the estimated DFT of the full data of $ x \left[ n \right] $.
The solution is then given by the Pseudo Inverse (Least Squares Solution):
$$ y = { ( \hat{F} \hat{F}^{T} ) }^{-1} \hat{F} x $$
In practice, the matrix will be very poorly conditioned hence solution must be generated using the LS Solution using the SVD.
A sample code is shared on GitHub Repository.
Result of the code:

Remark: Borrowed from Estimate the Discrete Fourier Transform / Series of a Signal with Missing Samples.


Why is $y(n) = x(n) + x(n-1)$ a low pass filter? I have found that this is a simple, but bad lowpass filter:

$$y(n) = x(n) + x(n-1)$$

However, I can't understand why it is a lowpass filter. What is its cutoff frequency?

What you have here is the equivalent of a moving-average filter. Specifically, it is a filter of order 1, whose impuse response is

$$h(n)=\delta(n)+\delta(n-1)$$

Taking its $Z$-transform, we get

$$H(z)=1+z^{-1}=\frac{z+1}{z}$$

There is a pole at $z=0$ and a zero at $z=-1$. Plotting the magnitude of the frequency response $H(\omega)\triangleq H(e^{-\imath \omega})=2\vert \cos(\omega/2)\vert$, you get the following curve

As you can see, this clearly is a low-pass filter. You can easily calculate the cut-off frequency from here on.


Product of power signal and energy signal I tend to believe that if we multiply an energy signal with a power signal the product is an energy signal.

My reasoning is that the energy signal will have finite bandwidth and consequenty the product will have finite bandwith and will be an energy signal.

Can anybody confirm that I am right? Or am I missing something?

Thanks in advance.

If you are not looking for a rigorous proof, then the following would suffice for a particular class of power signals which has the property that their magnitude is bounded; i.e, $ |x_p(t)| \leq K$ for all $t$.

Numerous examples exist for such power signals like sinusoidal waves, complex exponentials, all kinds of periodic continuous-time signals which are defined via convergent Fourier series sums... Then for such a power signal $x_p(t)$ and a given energy signal $x_e(t)$ you can show that their product $y(t) = x_p(t) x_e(t)$ will be an energy signal as follows; given

$$ E_x = \int_{-\infty}^{\infty} |x_e(t)|^2 dt $$

and $ |x_p(t)| \leq K $ for all $t$ then the energy of y(t) is

$$ E_y = \int_{-\infty}^{\infty} |y(t)|^2 dt = \int_{-\infty}^{\infty} |x_p(t) x_e(t)|^2 dt \leq \int_{-\infty}^{\infty} K^2 |x_e(t)|^2 dt = K^2 E_x$$

Hence $y(t)$ will also be an energy signal.


Does gray coding bring improvement to SER of QPSK, compared with natural coding? I understand gray coding contributes to BER of QPSK, yet really don't understand why it brings improvement to SER. According to the relationship of SER and BER, it should be improved. Actually, when mapping the symbol to constellation points, the minimum distance between two different points remains the same, and that induced the same SER. It really bothered me much.
Any comments would help!Thanks!

The central point as crystallized in the comments is that you think SER is derived from BER.

It's not!

SER is short for "symbol error rate", and that's it. That's the underlying thing.

A symbol error is when one symbol was sent, and a different one was received. It doesn't matter whether the symbol you've sent was "labeled" 00, or 11, or 10; or whether the received symbol was defined to be 00, 11 or 10, or yoghurt for that matter; what matters is that not the symbol that the sender sent was decided for in the receiver.

Thus, the SER doesn't care about what binary values are assigned to the symbols. (and that assignment can be Gray, or something else.)

By the way, I'm not sure what you mean with "natural encoding": I know that term from one-dimensional encoding problems (e.g. ADC values), where natural is just binary counting up:

value    |  0  1  2  3
---------+----------------
encoding | 00 01 10 11

If you extend it two a second dimension by just forming coordinate vectors, you get gray encoding for QPSK:

    \ dim0  0   1
      -------------
dim1 |
     |
   0 |     00  01
     |
   1 |     10  11

So, I'd argue Gray and natural are the same here, and using anything different would just be ill-advised .


A very basic question about processing high frequency signals A standard frequency for a modern computer processor is 1-3 GHz.
The frequency of a GPS signal is ~1.5 GHz.
My question is how is it possible for a computer to do anything at all with a signal that it can only sample once or twice in a period?

In digital signal processing, one almost never deals directly with high-frequency signals. The reason is that the frequency band that a signal occupies is completely irrelevant to the information that the signal conveys; the information exists in the envelope of the signal, which usually can be handled either by a generic CPU, a DSP, FPGA, or ASIC.

So, let's say that your GPS signal is $s(t)=x(t)\cos(2\pi f_c t)$, where $f_c\approx1.5\,\text{GHz}$. The actual information is in the envelope $x(t)$; the purpose of the cosine is simply to shift ("upconvert") $x(t)$ to some appropriate band. In the receiver, the signal $s(t)$ is "downconverted" to some lower (or even zero) frequency $f_I$, but this is done in the analog domain (i.e. by RF electronics), to obtain $s_1(t)=x(t)\cos(2\pi f_I t)$. The signal $s_1(t)$ has a small bandwidth, so it requires a small sampling rate, that can be handled in the digital domain.

Note that for quadrature signals the math is very slightly more complicated; see this question and the complex envlope article in wikipedia.

Note also that a similar process exists in the transmitter: the upconversion to RF is done in the analog domain, after the discrete envelope has been converted to an anlog signal.


What Are Linear and Circular Convolution? I have some basic understanding of signals and convolution. As far as I know it shows the similarities of two signals. Could I get some explanation in plain English of:

what are the linear and circular convolution
why they are important
practical situation where they are used

Linear convolution is the basic operation to calculate the output for any linear time invariant system given its input and its impulse response.

Circular convolution is the same thing but considering that the support of the signal is periodic (as in a circle, hence the name).

Most often it is considered because it is a mathematical consequence of the discrete Fourier transform (or discrete Fourier series to be precise):

One of the most efficient ways to implement convolution is by doing multiplication in the frequency.
Sampling in the frequency requires periodicity in the time domain.
However, due to the mathematical properties of the FFT this results in circular convolution.

The method needs to be properly modified so that linear convolution can be done (e.g. overlap-add method).


SPL values from FFT of microphone signal I am trying to get a frequency response vs dB SPL plot for recorded sounds. I have a microphone with a flat response, a sensitivity of 4mV/Pa, and 1 Pa = 94dB.

To get the dB SPL, I have the following formula:

$dB SPL = 94 + 20 \cdot \log_{10}(AI/0.004)$

where AI in the input voltage from the microphone, and assuming no gain.

This seems to work to get the overall dB SPL level for the signal, using peak voltage as AI.

My question is, how can I apply this to an FFT of the input signal, to get the intensity level per frequency?

I saw this post, and am not sure how to adapt this to input voltage.

Here is some python code to demonstrate what I have thus far, my goal is to apply this to ultrasonic signals that contain many different frequency components.

import numpy as np

fs = 5e5
duration = 1
npts = fs*duration
t = np.arange(npts, dtype=float)/fs
f = 1000
amp = 1.414*0.004
signal = amp * np.sin((f*duration) * np.linspace(0, 2*np.pi, npts))

rms = np.sqrt(np.mean(signal**2))
dbspl = 94 + 20*np.log10(rms/0.004)

sp = np.fft.fft(signal)
freq = np.arange(npts)/(npts/fs)

sp_db = 94 + 20*np.log10(???)

Let's assume that signal you are analysing is sinusoid with amplitude $A$:

$x=a\sin{2\pi f_{0} t}$

Its RMS value of amplitude is then: $\dfrac{a}{\sqrt{2}}$ as you noticed in your code. Before performing DFT it is good practice to window signal as it decreases leakage, etc.

After transforming this signal into the frequency domain, you will obtain a two-sided spectrum. Then we have to take the magnitude of these complex values. Probably you would notice that all of them are very high. That's because energy is not normalised. Normalisation can be done via dividing values of magnitude by energy of your window. As people tend to take signal "as it is" (which is equal to applying rectangular window) the spectrum is then divided by a number of samples $N$. For different types of windows this factor that is equal to the sum of all window samples.

Additionally, we are only interested in one half of a spectrum, therefore the amplitude of all samples must be multiplied by $2$ to compensate the loss of energy, except DC component (which appears only once) and Nyquist bin (if it's present).

Lastly, as you want to analyse RMS of the signal, you must understand the following. All values of single-sided spectrum are of height $\dfrac{A_{i}^{2}}{2}$, which is equal to $ \left( \dfrac{A_{i}}{\sqrt{2}}\right)^{2}$, where $\dfrac{A_{i}}{\sqrt{2}}$ is the RMS magnitude amplitude for the i-th frequency component. Translating this into Python code you obtain:

import numpy as np
import matplotlib.pyplot as plt

fs = 5e5       # Sampling frequency [Hz]
duration = 10  # Duration of the test signal [s]

N = int(fs * duration)
t = np.arange(N, dtype=float) / fs

f = 1000  # Frequency of the test signal [Hz]
sensitivity = 0.004  # Sensitivity of the microphone [mV/Pa]
p_ref = 2e-5 # Reference acoustic pressure [Pa]
amp = sensitivity * np.sqrt(2)  # Amplitude of sinusoidal signal with RMS of 4 mV (94 dB SPL)
signal = amp * np.sin(2 * np.pi * f * duration * t)  # Signal [V]

# Calculate the level from time domain signal
rms_time = np.sqrt(np.mean(signal**2))
db_time = 20 * np.log10(rms_time / sensitivity / p_ref)

# Apply window to the signal
win = np.hamming(N)
signal = signal * win

# Get the spectrum and shift it so that DC is in the middle
spectrum = np.fft.fftshift( np.fft.fft(signal) )
freq = np.fft.fftshift( np.fft.fftfreq(N, 1 / fs) )

# Take only the positive frequencies
spectrum = spectrum[N//2:]
freq = freq[N//2:]

# Since we just removed the energy in negative frequencies, account for that
spectrum *= 2
# If there is even number of samples, do not normalize the Nyquist bin
if N % 2 == 0:
    spectrum[-1] /= 2

# Scale the magnitude of FFT by window energy
spectrum_mag = np.abs(spectrum) / np.sum(win)

# To obtain RMS values, divide by sqrt(2)
spectrum_rms = spectrum_mag / np.sqrt(2)
# Do not scale the DC component
spectrum_rms[0] *= np.sqrt(2) / 2

# Convert to decibel scale
spectrum_db = 20 * np.log10(spectrum_rms / sensitivity / p_ref)

# Compare the outputs
print("Difference in levels: {} dB".format(db_time - spectrum_db.max()))

plt.semilogx(freq, spectrum_db)
plt.xlim( (1, fs/2)  )
plt.ylim( (0, 120) )
plt.grid('on')
plt.xlabel("Frequency [Hz]")
plt.ylabel("SPL [dB]")

plt.show()

  Difference in levels: -1.4210854715202004e-14 dB

As you can see, this is producing consistent results for both time and frequency domain. Good luck!


Two-View 3D reconstruction using the sparse Levenberg–Marquardt algorithm I have trouble implementing the Levenberg-Marguardt algorithm as described in the book Multi View Geometry in Computer Vision. To be more specific I have trouble calculating the partial derivatives that are needed in order to pursue with the main algorithm implementation.

Consider having the projection vector in two views $\mathbf{X}=[\mathbf{x},\mathbf{y},\mathbf{x}',\mathbf{y}']$ where the $x,y$ are the projected coordinates in the first view and $\mathbf{x}′,\mathbf{y}′$ are the coordinates in the second view respectively. Also lets say we are given an initial estimated $\mathbf{A}$ Matrix with 12 parameters describing the rotation and the translation of the second view with the respect to the first view and an initial estimation for the point in the 3D coordinates $\mathbf{b}=[\mathbf{X},\mathbf{Y},\mathbf{Z}]$.

How can the matrices $\partial\mathbf{X}/\partial{A}$ and $\partial\mathbf{X}/\partial{b}$ can be calculated ?

PS: I tried searching for numerical differentiation methods but I had hard time understanding them i.e finite difference.

First, please consult this tutorial about the LM minimization and this document on the bundle adjustment.

Along with many other derivative based descent algorithms, the Levenberg-Marquardt algorithm relies on the partial derivative matrix, a.k.a. Jacobian Matrix, which is the matrix of all first-order partial derivatives of a vector-valued function:

$$
\begin{align}
J_i &= \frac{\partial f(x_i, \beta)}{\beta} \\
\mathbf{J} &= \begin{bmatrix}
    \dfrac{\partial f_1}{\partial x_1} & \cdots & \dfrac{\partial f_1}{\partial x_n}\\
    \vdots & \ddots & \vdots\\
    \dfrac{\partial f_m}{\partial x_1} & \cdots & \dfrac{\partial f_m}{\partial x_n} \end{bmatrix}
\end{align}
$$

Now how you compute Jacobian depends on what you actually want from the minimization. You could for example do a full bundle adjustment where Jacobian might have an involved expression. But, let me give an example of the simplest scheme I could think of.
Triangulation from multiple views is generally written in the form of minimization of the reprojection error. Here is the cost function for the two view problem:

$$
C(X) = d(x,\hat{x})^2 + d(x',\hat{x}')^2
$$
where $x$ and $x'$ are points in 2D views and $\hat{x}$ and $\hat{x}'$ are the projections of 3D points. In the matrix notation, the error function reads:

$$
\begin{pmatrix}
\epsilon_1 \\
\vdots \\
\epsilon_m \\
\end{pmatrix} =
\begin{pmatrix}
\|x_1^a-P_a X^1\|+ \|x_1^b-P_b X^1\|\\
\vdots \\
\|x_m^a-P_a X^m\|+\|x_m^b-P_b X^m\| \\
\end{pmatrix}
$$
where $P_a, P_b$ are the projection matrices of cameras $a$ and $b$. Since we are searching for the 3D points only $\{\mathbf{X}_i\}$, the Jacobian is $Mx3$:

$$
\begin{align}
\mathbf{J} &= \begin{bmatrix}
    \dfrac{\partial \epsilon_1}{\partial x} & \dfrac{\partial \epsilon_1}{\partial y} & \dfrac{\partial \epsilon_1}{\partial z}\\
    \vdots & \ddots & \vdots\\
    \dfrac{\partial \epsilon_m}{\partial x} & \dfrac{\partial \epsilon_m}{\partial y} & \dfrac{\partial \epsilon_m}{\partial z} \end{bmatrix}
\end{align}
$$

Of course this problem requires a good initial solution, which could be obtained from a linear closed form solution. The quality of the non-linear minimizer is highly dependent on the initialization and the quality of 2D points.

For the rest, and different parameterizations, check out the links.

Here are some remarks:

If you are not looking for re-inventing the wheel, then I would recommend Ceres Solver, which includes Automatic Differentiation. This saves the headache of Jacobians (and is much more accurate than numerical derivatives - numerical approximations are highly discouraged). Ceres also includes samples for bundle adjustment.
If you are also optimizing for the pose (camera orientation), then I do not recommend operating on the 12 parameter vector of rotation and translations. Rather choose a representation such as quaternions or the exponential map (twist coordinates). As this preserves the manifold-ness, you will get much higher quality solutions, without ambiguities. Some examples are given in the links I provided.


Filters using equiripple approximations (different from type I lowpass) In Oppenheim & Schafer's Discrete-Time Signal Processing, only type I lowpass filters are explained in detail in the chapter "Optimum Approximations of FIR filters".
I want to know what happens not only with the other three types of generalized-phase filters, but also understand how the equiripple approximation works with bandpass and highpass filters.
In the previously mentioned book, it says:

Specifically, we will show that for type I lowpass filters:

The maximum possible number of alternations of the error is $L+3$.
Alternations will always occur at $\omega_p$ and $\omega_s$.
[...] the filter will be equiripple, except possibly at $\omega=0$ and $\omega =\pi$.

I have a bunch of questions, for example:

Is there a maximum possible number of alternations for every type (II, III, IV) and kind (lowpass, highpass, bandpass) of filter?
Do the other two properties (alternations at edges and equiripple-ness except possibly at $0$ and $\pi$) apply to other types and kinds of filters?
Do other types and kinds of filters have other special properties (i.e., properties that differ from the ones in the quote I wrote)?

It would be great if you could give me some insight in the subject, or if you could recommend a book to which I could fall back on.

The only place where I found the answer to this question was in Theory and Application of Digital Signal Processing by Rabiner & Gold.

In its chapter about theory and approximation of FIR filters, there is a section where all of my questions are explicitely answered:


Creation of Brownian noise as random walk I am trying to create brown noise by integration of white noise (random walk).

I created Gaussian white noise. The histogram shows Gaussian distribution and Welch (and also FFT) shows flat spectrum.

Now I try to integrate it

$x(k+1) = x(k) + q(k)$

where $q$ stands for white noise and $x$ for resulting Brownian noise.

But the resulting spectrum is

According Wikipedia I would expect that there should be straight line in the plot.

What I am doing wrong?

So this question does not remain as unanswered:

Making your horizontal scale into log-scale should give you the decade drop and yield the linear amplitude graph.


Why downsample exhibit wider bandwidth when doing discrete time fourier transform? Based on Nyquist sampling theorem, sampling would have a convolution in frequency domain, and reasonably the bandwidth of each convolution would be as that of the original signal.
Nevertheless, the discrete time Fourier transform of a signal ,which has been  downsampled , has the $M$ multiples bandwidth of each convolution.($M$ is the downsample factor). I could not have got the idea, what is the true relation between time/frequency domain with underlying discrete time Fourier transform?

Sampling gets convolution in frequency domain:

Downsampling gets more wider bandwidth:

You don't get "more" bandwidth. That's simply a misunderstanding of your x-axis. Your signal after downsampling occupies a $M$ times larger part of the Nyquist rate, but that's because the Nyquist rate got reduced by a factor of $M$, not because the signal got wider!


How does a digital filter work? To introduce my situation: I'm developing a digital synthesizer in a form of a C++ library, working with low level APIs like WASAPI, ASIO, ALSA etc. It's probably not very practical and I'm mostly "reinventing the wheel" but my intention is to learn about digital synthesis in depth. So far I have successfully implemented basic concepts like oscillators and modulation of their properties. The next logical step is a filter.

So my question is: How does a digital filter work on this low level? How exactly does it modify the individual samples?

I understand, that this involves a lot of math. That's not a problem for me. I only need a good starting point (some sources to learn) and an intuitive explanation as all I was able to find were either analogue explanations or just formulas explained with a lot of advanced terminology that I'm not familiar with.

Consider a moving average over N samples- this is a simple FIR filter where each new output is the average of the past N samples. It is easy to see how high frequency noise can be filtered out (so is a low pass filter), and the longer time duration we include in the averaging window the lower will be the frequency cut off (just compare a stock market 30 day moving average to a 1 day moving average).

A moving average is a poor low pass filter, having a frequency response that approaches a Sinc function, which rolls off relatively slowly in frequency. By doing a weighted moving average where different samples are given different weights in the averaging process, we can significantly improve the frequency response - and coming up with the correct weights is the science of digital filter design.

IIR filters are similar except we are performing the average with previous outputs instead of past inputs.


Numerically Stable IIR filter I am making an IIR Filter and I am running into problems with stability with higher orders. The values go off to infinity as soon as it starts.

So far, I've been using MATLAB to generate my filters and then running them in an embedded system. I'm using a band-pass filter, state-space representation and single precision values. If I understand this correctly, the issues happens because the poles are getting too close to the unit circle, so the rounding errors on my coefficients or state space are getting out of hand.

I could reduce the order, but that feels wasteful when I have the computational power to run a higher order filter. Is there a filter design that will improve at a higher order while giving the poles a good margin?

Alternatively, Is there another trick I can use to mitigate this issue?

First of all, what is the order of your IIR filter? The highest order I have ever used was an order-10 IIR filter for a control loop application. I feel like it is unlikely that you need more that this.

Second, it is a good idea to split your filter in second-order-sections (SOS) and cascade them , this usually fix most issues.

https://www.dsprelated.com/showarticle/1137.php

You can also split your filter in second-order-sections and put them in parallel and add their ouputs too.


Why is In-phase Quadrature sampling not used to record and store digital audio? I/Q is a standard method to sample radio frequency and its use is even seen in some non radio signals like sonar. But it seems that it is not used to sample acoustic frequencies.

What property of sound makes plain in-phase PCM or PDM adequate for recording and storing audio–even for high-fidelity recording?

I have an intuition as to the answer, but I think satisfying explanation would be useful to those researching this area of signal processing.

It is because the audio signals are real and already at baseband. In contrast radio frequency signals are often represented as complex numbers once they are brought back to baseband. Real signals can be represented as a single stream of real numbers, while for complex numbers two streams of real numbers are required to represent them (as in $I+jQ$).

When the upper sideband of a radio waveform at a given carrier frequency does not match the lower sideband of the waveform in magnitude and opposite phase, then the waveform when brought to baseband must be represented with a complex signal (because at baseband the positive and negative frequencies will not be complex conjugate symmetric, which is the requirement for any real signal). A simple example of this is Quadrature Phase Shift Keying and Quadrature Amplitude Modulation where samples are mapped to locations on the complex (IQ) plane as part of the modulation.

It is helpful to have a full grasp of a complex frequency represented as $e^{j\omega t}$ vs $cos(\omega t)$ and what happens when you translate either in frequency to a real carrier: in the case of the exponential only a positive frequency exists so when translated using a real carrier frequency will be a single real tone at a higher frequency than the carrier (upper side band) and similarly $e^{-j\omega t}$ will be a single real frequency at a lower frequency than the carrier (lower side band). The sinusoid in contrast contains both exponentials (see Euler’s identify) and therefore will have both sidebands when translated. This shows how the complex form can be used to occupy half the bandwidth which is a motivation for doing so (in fact the higher modulation orders mentioned offer better spectral efficiency and are complex waveforms in order to do this.)

Similarly the audio signal could be made complex if it was to be converted to a carrier for the same reason (for example single-sided AM is one variant of an analog modulation that can be used with audio in complex form.)

But ultimately an audio signal that is going to one speaker is represented by a single stream of real signals so as explained earlier is not complex — but between the microphone and the speaker we can choose any way to encode the data including complex signals.


Interpreting FFT Coefficients from System Matrix That Is Originally Toeplitz (Not Circulant) If I have a measured signal $y$, true signal $x$, and a convolution matrix $A$ that is a Toeplitz but not circulant matrix, I can write the convolution as

\begin{equation}
y = A x \ .
\end{equation}

However, I would like to analyze the system using diagonalized coefficients of the Fourier transform. A re-formulation such as the one given below is often used to provide a means for fast computation:
\begin{equation}
\left[\begin{array}{c}
y\\
y'
\end{array}\right] =
C
\left[\begin{array}{c}
x\\
0
\end{array}\right] =
\left[\begin{array}{cc}
A & B \\
B & A
\end{array}\right]
\left[\begin{array}{c}
x\\
0
\end{array}\right] =
\left[\begin{array}{c}
A x \\
B x
\end{array}\right] \ ,
\end{equation}
where $C$ is a circulant matrix into which the Toeplitz matrix $A$ is embedded, using an additional matrix $B$ derived from the values of $A$.

Diagnoalized values of $C$ can be obtained by multiplication with Discrete Fourier Transform matrices. But are these values relevant for interpreting the behavior of the original system, $y = Ax$ - for instance, for use in construction of a Wiener filter for the original system?

If I understand your question, then the answer lies in the equivalence of linear and circular convolution - the linear convolution being implemented by multiplication with a toeplitz matrix, while circular convolution is implemented by multiplication with a circulant matrix.

If you have a length N signal and length M filter impulse response, the circular and linear convolutions are equivalent if you pad the signal and filter with zeros to length N+M-1 - that is, you can implement the linear convolution using the DFT. You are padding x with zeros above, and presumably the number of zeros is such that the length N+M-1 requirement is met.  Then in matrix form, the DFT matrix would diagonalize the corresponding convolution operator.

Hope this helps - I'm short on time right now and can't provide a more detailed explanation.


Meaning of Multiplying a complex number with its conjugate While reading Estimation of Symbol Timing Offset, i came across an equation like multiplying a signal with its conjugate.

$$
\hat{\delta} = \arg \max_{\delta} \sum_{i=\delta}^{N_g-1+\delta} \left|y^*(n+i) y(n+N+i)\right|^2
$$
I know multiplying

(5+4i) X (5-4i) = 25 - 20i + 20i + 16 = 25 + 16 = 41.

But am not getting the meaning of this 41. What it represents ???

Multiplying a complex number with it's conjugate yields the squared magnitude of that complex number.

However, in your formula, you are multplying the signal with the conjugate of the N samples delayed signal, take the squared magnitude of that (thus dropping phase differences between the signal and delayed signal, which would result from a frequency offset) and average it over Ng sample times. This means you are calculating the autocorrelation of the signal at a time delay of N sample times.

I assume the formula is from an algorithm for coarse symbol timing estimation for OFDM signals. The pseudo-random pilot sequence of length N is repeated twice. It's autocorrelation thus has a peak at N. Now, if you want to estimate the timing offset, you just take the time index of this peak.


Principle of Mapper in frequency diversity in OFDM I am trying to understand frequency diversity in OFDM. An extra block 'MAPPER' is introduced in conventional OFDM transmitter to apply frequency diversity. Can anyone explain the principle of this MAPPER with its working and mathematical formulation. The transmitter and Reciever block diagrams are as shown below:

In transmitter the N symbols from the series to parallel convertor are repeated M times at the mapper giving you M sets of N symbols. This is a guess, but I  assume the mapper is spreading the contents of each of these M sets out across the whole transmit frequency band so the repeated symbols are not concentrated in any one region. This way if any there is a null in any particular region of the band only a few of the repeated symbols would be affected. This would give you frequency diversity.

AT the receiver the contents of each of the M sets are weighted and combined in a kind of majority rule detection scheme to determine the estimated N symbols.


Does horizontal motion blur imply the shape of PSF is rectangular pulse? For an image has a simulated horizontal motion blur on it, is it reasonable for someone to say that the suggested shape for the point-spread function resulting in this motion blur is a rectangular pulse aligned in the direction of motion?

If so, what is an explanation for this?

Think a minute of what a motion blur means for the physical pixels: they are moved during the acquisition process.
Assuming left-to-right horizontal motion for simplicity, then what happens is that one pixel will measure the amount of light in its initial position, then the light in the position that was occupied before by its right neighbour, and so on until the frame acquisition ends.

Mathematically, the light measured by one pixel is then a mixture between the value in its original position and the values of its neighbours in their original positions.
Hence, up to a rescaling of the output value (that is unnecessary of you normalize the blur kernel), an horizontal motion blur can be modelled by an horizontal blur kernel.

Of course, this is grossly simplified (I did intentionally forgot about minimum exposure time for example and supposed instead that the measurement was instantaneous) that assumes that the scene is at a fixed distance from the camera and that all the image points are at the same distance, without independantly moving objects.


Separating/recovering base signal from two mixed signals, given phase information I have collected two signals, $B_1(x)$ and $B_2(x)$. The signals start and end at the same $x$-values. Each signal $B_i(x)$ contains:

a base signal $b(x)$, which is the same for both, and
a signal, either $w(x)$ or $w(x+\phi)$, which are identical except for the shift $\phi$ along $x$.

Signals $b(x)$, $w(x)$, and $w(x+\phi)$ are non-periodic.
Essentially, I have:
$B_1(x) = b(x) + w(x) \\
B_2(x) = b(x) + w(x + \phi)$,
with only access to $B_1(x)$, $B_2(x)$, and $\phi$.
My question: Is is possible to use the phase-shift knowledge from my two signals to recover the base signal $b(x)$? And if so, can you point me in the direction of a method that would work?
I've read about or tried to implement methods such as blind source separation, independent component analysis, and stationary subspace analysis, but none seem to take advantage of the phase information, or work at all for my data.
Note: If it helps, I am able to collect more data with different $\phi$-values, to obtain more $B_i(x)$ signals.
I've attached a figure with toy data (my real data is messier and longer). The top plot shows $b(x)$ (what I'd like to recover), the middle plot shows $w(x)$ and $w(x+\phi)$, and the bottom plot shows the measured signals $B_1(x)$ and $B_2(x)$, representing the combination.

Here's the best method I've found so far:
First, I shifted $B_2(x)$ by $\phi$; this shift aligns the underlying (but unknown) signals $w(x)$ and $w(x + \phi)$. Note that once $B_2(x)$ is shifted, the signals will not begin and end at the same $x$-values. I chopped off the values of each signal that didn't overlap. Parts (a) and (b) of the figure below show the original measured signals $B_1(x)$ and $B_2(x)$, and then the aligned signals, which I'll call $s_1(x)$ and $s_2(x)$.
The aligned signals $s_1(x)$ and $s_2(x)$ are not identical, because the base signal $b(x)$ has been added to them at a known phase shift. Therefore, finding the difference between $s_1(x)$ and $s_2(x)$ at each point, and then dividing by the phase shift $\phi$, will give us the numerical derivative of the base signal, $\frac{\triangle b}{\triangle x}$. In summary, the numerical derivative of the base signal is given by
\begin{equation}
 \frac{\triangle b}{\triangle x} = \frac{s_2(x)-s_1(x)}{\phi}.
\end{equation}
Then, to find the estimated base signal $b_{est}(x)$, we numerically integrate $\frac{s_2(x)-s_1(x)}{\phi}$. I corrected for integration drift by numerically integrating forwards and backwards, then averaging the two results.
Lastly, I shifted $b_{est}(x)$ forward by $\frac{\phi}{2}$. I think the reason for a shift of this value is because I computed the numerical derivative between points separated by a full phase, so the derivative approximate should be centered at half of that.
The results (base signal $b(x)$, estimated base signal $b_{est}(x)$, and error between the two) are shown in part (c) of the figure below.


Questions about log polar transform I read various papers about the log polar transform and its application on template matching with images and have some questions:

In "Image Registration Using Log Polar Transform and Phase Correlation to Recover Higher Scale" the authors say:
"In this algorithm first the sense image is downscaled by the factor of 2. Log-Polar algorithm gives best result for any value of rotation and translation but it will not work if the scale is more than 1.8 or 2, as the frequency components alters at high scale [7]." In the referenced paper "IMAGE REGISTRATION USING LOG-POLAR TRANSFORM AND  PHASE CORRELATION" the authors only mention at the end that their results are verified with a scaling factor up to 2.0, there is no mention at all about altering frequency components at high scale.
Then I read through http://cvpr.uni-muenster.de/teaching/ss08/seminarSS08/downloads/Wentker-Vortrag.pdf . He says that the polar transform is unable to detect scales > 2, the log transform doesn't have it. This clearly contradicts the first paper, which says the log polar transform can only go to scales < 2.
So does the log polar transform normally only work up to scale 2.0?
Regarding question 1, what exactly is the reason that it only works up to a scale of 2? The altering frequency components are mentioned, but I don't know what they really mean with that and it is never explained in any of the papers I found, only that this value exists.
In the presentation http://cvpr.uni-muenster.de/teaching/ss08/seminarSS08/downloads/Wentker-Vortrag.pdf page 22, we have the normal, polar and log-polar transformed image. The polar image has the black areas, which are parts that can't be mapped to the original image because it would be outside of its boundaries. However, these areas have nearly disappeared in the log polar image. Why is that? I tested various Logpolar implementations using OpenCV, Python, another example is Scale and Rotation invariant Template Matching and all the time, the log polar images looks more like the middle image and not like the right one.
What exactly is happening here? Why does the third image look so different, compared to other log polar images which look more like the second?

I hope you can answer my questions.

3.

  The polar image has the black areas, which are parts that can't be
  mapped to the original image because it would be outside of its
  boundaries. However, these areas have nearly disappeared in the log
  polar image. Why is that?

Look at the second (polar) and third (log-polar) images.  What's different about them?

In the 3rd image, stuff that was on the left is now in the center.  Stuff that was in the center is now on the right.  The "log" in log-polar means that you're stretching the image in a logarithmic way.  So the black out-of-image scallops are still there, they've just been squished very thin against the right side of the image, while details from the left side of the image have been stretched out to take their place.


inverse discrete FFT in python, multiple times? I was wondering what really happens when taking the inverse discrete FFT on some set of numbers, for 3 times? Because looking at it, it looks like we're getting an output that is identically with the input, and I don't get why.

Also I was trying to implement my own inverse discrete Fourier transform in python, following wikipedia, looking at the IFFT provided by octave:

and at the plot provided by my code:

It just doesn't make any sense.
Here is my code:

#!/usr/bin/python
import matplotlib.pyplot as plt
import numpy as np
import math

points = [24+0j, -18.39230+0.73205j, 8-3.46410j,
        -2+4j, 0+0j, 2.39230-2.73205j, -4+0j, 2.39230+2.73205j, 0-0j,
        -2-4j, 8+3.46410j, -18.39230-0.73205j]

def inverse_fourier_transform(points):
    inversed_points = []
    N = len(points)

    for k in range(0, N):
        s = 0.0
        for n in range(0, N - 1):
            s += points[n]*np.exp(1.j*2.0*np.pi*k*n/N)/N
        inversed_points.append(s)

    return inversed_points

inv_points = inverse_fourier_transform(points)
inv_points2 = inverse_fourier_transform(inv_points)
inv_points3 = inverse_fourier_transform(inv_points2)

x = []
y = []

for pt in inv_points3:
    x.append(pt.real)
    y.append(pt.imag)

plt.scatter(x, y, s=0.5)
plt.plot(x, y)
plt.show()

The FFT and IFFT are almost identical operation: other than scaling (depending on convention) and  complex conjugation they are the same.

Ignoring the scaling (or assuming $1/\sqrt N$). We have

Apply one time: spectrum
Apply second time: time inverse original
Apply third time: conjugate spectrum
Apply 4th time: original signal

For inverse inverse FFT: reverse the order of the steps.


What does the 'B' parameter represent in Wikipedia's page on window functions? I was trying to learn some more window functions and was referencing the Wikipedia article on the subject. I've noticed on a lot of the plots, there appears to be a capitalized and italicized letter B. What is this parameter? I can't find any reference to it on the page. link to wikipedia article And here's a photo that features two plots, each having an associated 'B' parameter with them. imgur What is this?

See paragraph from wiki. There is an explanation of B metric:

  Each figure label includes the corresponding noise equivalent bandwidth metric (B),[note 1] in units of DFT bins.

Note 1:

  Mathematically, the noise equivalent bandwidth of transfer function H is the bandwidth of an ideal rectangular filter with the same peak gain as H that would pass the same power with white noise input. In the units of frequency f (e.g. hertz), it is given by:
  $$
B_{noise}={\frac {1}{|H(f)|_{max}^{2}}}\int _{0}^{\infty }|H(f)|^{2}df.
$$


Fourier Like Spectral Analysis with Uneven Intervals and Redesigned DFT Matrix I intended to use a discrete Fourier transform (DFT) on a time series sampled at uneven intervals. What I did was to calculate a DFT matrix where the elements are the values at the uneven locations like this
$$
{\displaystyle {\begin{aligned}X_{k}&= \frac{1}{\sqrt{N}}\sum _{n=0}^{N-1}x_{n}\cdot e^{-{\frac {2\pi i}{time_{N-1} - time_0}}k \cdot time_n}\\&\end{aligned}}}
$$
where $time_n$ is the time location in the sequence and $time_{N-1} - time_0$ is the entire time span of all samples.
I realize that orthogonality of the basis vectors in the DFT matrix might get somewhat lost. Is anything more lost? Can I do something more to make the transformation more comparable to an ordinary DFT?
Here is an implementation of it you can try with your own data: David Jonsson
(bodavid) - Fourier Transform for Unevenly Sampled Data.
Orthogonality between basis vectors are visualized in a colorized matrix.

I have derived the DFT for data whihc is sampled in a non uniform manner:
The DFT Matrix for Non Uniform Time Samples Series
Problem Statement
We have a signal $ x \left( t \right) $ defined on the interval $ \left[ {T}_{1}, {T}_{2} \right] $.
Assume we have $ N $ samples of it given by $ \left\{ x \left( {t}_{i} \right) \right\}_{i = 0}^{N - 1} $. The samples time $ {t}_{i} $ is arbitrary and not necessarily uniform.
We're after the DFT of the samples $ \left\{ X \left[ k \right] \right\}_{k = 0}^{K - 1} $ as it was samples in a uniform manner (Implicitly means the samples in Frequency domain will be uniform as well).
Deriving the Connection
In the DFT Transform the connection between time and frequency is given by:
$$ x \left[ n \right] = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N} n } \tag{1}  $$
In $ \eqref{EqnIdft} $ we use $ n $ for modeling the sample index in time. We usually build samples in time as $ x \left[ n \right] = x \left( n {T}_{s} \right) $ where $ {T}_{s} $ is a uniform sampling interval.
Hence we could write:
$$ x \left( n {T}_{s} \right) = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N {T}_{s}} n {T}_{s}} \tag{2} $$
In $ \eqref{EqnIdft2} $ we added explicit scaling of time. This is a known property of Fourier transform family which scales the domain in order to normalize the transform.
Now, there is nothing which blocks us from using arbitrary time:
$$\begin{align*} \tag{3}
x \left( t \right) & = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k}{N {T}_{s}} t} && \text{} \\
& = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{j 2 \pi \frac{k {F}_{s}}{N} t} && \text{Since $ {F}_{s} = \frac{1}{{T}_{s}} $}
\end{align*}$$
As can be seen $ \eqref{EqnIdft3} $ makes sense as it goes through each element according to its frequency and sums to give the output at time $ t $. We can go step farther and generalize it for cases we don't have uniform sampling frequency.
The average sampling frequency is given by $ \bar{F}_{s} = \frac{N}{ {T}_{2} - {T}_{1} } $. Let's define $ T = {T}_{2} - {T}_{1} $ and we'll get:
$$ x \left( t \right) = \frac{1}{N} \sum_{k = 0}^{N - 1} X \left[ k \right] {e}^{ j 2 \pi k \frac{t}{T} } $$
Which is many ways resembles the DTFT Transform equation which does the same in the other direction, transforming uniform discrete samples in time domain to arbitrary frequency (Within a frequency interval) in Frequency Domain:
$$\begin{align*} \tag{4}
X \left( f \right) & = \sum_{n = 0}^{N - 1} x \left[ n \right] {e}^{-j 2 \pi f {T}_{s} n } && \text{} \\
& = \sum_{n = 0}^{N - 1} x \left[ n \right] {e}^{-j 2 \pi \frac{f}{ {F}_{s} } n } && \text{Since $ {F}_{s} = \frac{1}{{T}_{s}} $}
\end{align*}$$
We see the same scaling, $ \frac{f}{ {F}_{s} } $ which scales the continuous $ f $ relative to the interval of frequencies $ {F}_{s} $ which is equivalent to $ \frac{t}{ T } $ which scales $ t $ relative to the time interval of the continuous signal.
The Transform Matrix
So, given the set of time indices $ {\left\{ {t}_{i} \right\}}_{i = 0}^{N - 1} $ the transformation matrix, from frequency domain to time domain, is given by:
$$ D \in \mathbb{R}^{N \times K}, \;  {D}_{i, k} = {e}^{ j 2 \pi k \frac{ {t}_{i} }{T} } $$
The Model
In vector form the model is:
$$ x = D y $$
Where $ y \in \mathbb{C}^{K} $ is the vector of the frequency coefficients in uniform grid, $ x $ is the samples in time (Non Uniform, Or at least no assumption of uniformity) and $ D $ as defined above.
Since in our model we're after $ y $ the answer is given by:
$$ y = {D}^{\dagger} x $$
Where $ {D}^{\dagger} $ is the Pseudo Inverse Matrix of $ D $.
Implementation & Results
The code is as following:
subStreamNumberDefault = 79;

run('InitScript.m');

figureIdx           = 0;
figureCounterSpec   = '%04d';

generateFigures = ON;

%% Simulation Parameters

samplingFrequency = 101; %<! [Hz]
samplingInterval = 1 / samplingFrequency; %<! [Sec]
startTime = 1; %<! [Sec]
endTime = 4; %<! [Sec]
timeInterval = endTime - startTime; %<! [Sec]

numSamples = round(samplingFrequency * timeInterval);
numSamplesTT = round(1.2 * numSamples);

signalFreq = 2; %!< [Hz]

% The uniform time grid
vT      = linspace(startTime, endTime, numSamples + 1);
vT(end) = [];
vT      = vT(:);

% The non uniform time grid - Reconstruction
vTT = endTime * rand(numSamplesTT, 1);
vTT = sort(vTT, 'ascend');

% The non uniform time grid - DFT
vTD = linspace(startTime, endTime, (10 * numSamples) + 1);
vTD(end) = [];
vTD = vTD(sort(randperm(length(vTD), numSamples)));
vTD = vTD(:);

% The uniform frequency grid
vF      = (samplingFrequency / 2) * linspace(-1, 1, numSamples + 1);
vF(end) = [];
vF      = vF(:);

vK = [-floor(numSamples / 2):floor((numSamples - 1) / 2)];
vK = vK(:);

%% Generate Data

vX  = cos(2 * pi * signalFreq * vT);
vFx = fftshift(fft(vX));

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = subplot(1, 2, 1);
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['Reference Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);

hAxes           = subplot(1, 2, 2);
hStemObj = stem(vF, abs(vFx));
set(hStemObj, 'LineWidth', lineWidthNormal);
set(get(hAxes, 'Title'), 'String', {['DFT of the Reference Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Frequency [Hz]']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Magnitude']}, ...
    'FontSize', fontSizeTitle);

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Analysis - Reconstruction

mD = exp(1j * 2 * pi * (vTT / timeInterval) * vK.') / numSamples;

% Reconstruction according to the model
vY = real(mD * vFx);

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries     = plot(vTT, vY);
set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

%% Analysis - DFT of the Non Uniformly Sampled Data

vY  = cos(2 * pi * signalFreq * vTD);

mD = exp(1j * 2 * pi * (vTD / timeInterval) * vK.') / numSamples;
vFy = pinv(mD) * vY;

figureIdx = figureIdx + 1;

hFigure         = figure('Position', figPosLarge);
hAxes           = axes();
set(hAxes, 'NextPlot', 'add');
hLineSeries     = plot(vT, vX);
set(hLineSeries, 'LineWidth', lineWidthNormal);
hLineSeries     = plot(vTD, vY);
set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Time Index']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Sample Value']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

figureIdx = figureIdx + 1;

hFigure     = figure('Position', figPosLarge);
hAxes       = axes();
set(hAxes, 'NextPlot', 'add');
hStemObj    = stem(vF, abs([vFx, vFy]));
set(hStemObj, 'LineWidth', lineWidthNormal);
% hLineSeries     = plot(vTT, vY);
% set(hLineSeries, 'LineWidth', lineWidthNormal, 'LineStyle', ':', 'Marker', '*');
set(get(hAxes, 'Title'), 'String', {['DFT of the Uniform Signal & Non Uniform Signal']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Frequency [Hz]']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'YLabel'), 'String', {['Magnitude']}, ...
    'FontSize', fontSizeTitle);
hLegend = ClickableLegend({['Uniform Signal'], ['Non Uniform Signal']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

Results are:

Summary
In this post we derived how to estimate the Uniform DFT of a Non Uniform Time Series by solving linear system of equations.
The full code is available on my StackExchange Signal Processing Q32137 GitHub Repository.
Remark: Why Do We Need to Apply fftshift() on the DFT of the Signal?
Indeed in the Reconstruction part we use fftshift(). The shallow answer is easy, we also build the vector vK as symmetric around zero.
But there is a deeper reason for that. In the DFT when we use uniform sampling in Frequency Domain and Time Domain Magic happens without us seeing it explicitly.
When we defined the term $ \frac{k}{ N {T}_{s} } n {T}_{s} $ we replaces $ n {T}_{s} $ with $ t $ hence we prevent the term $ {T}_{s} $ to cancel itself. Now setting $ {F}_{s} = \frac{1}{{T}_{s}} $ means that we multiply by $ k $ and we get frequencies which are out of the Nyquist Frequency.
In most cases when we that happens the Modulo property of the exponent comes in and we get the correct negative value of the frequency in the range $ \left[ -\pi, \pi \right] $. Yet when $ t $ is arbitrary we can think that $ {F}_{s} $ is changing per sample which means when we go farther than $ \pi $ the modulo doesn't bring us to the correct answer.
First, as intuition, always think the DFT is defined on the $ \left[ -\pi, \pi \right] $ interval and it is continuous. So as long as you work on this range things works as intended. This intuition can come from the Fourier Series and Discrete Fourier Series (DFS).
Let's try explaining it using a concrete example. Let's examine the exponent term from the derivation:
$$ 2 \pi \frac{k}{N {T}_{s}} n {T}_{S} = 2 \pi \frac{k}{{F}_{S}} \frac{{F}_{s}}{N} n = 2 \pi \frac{k b}{{F}_{s}} n $$
Where $ b $ is the Bin Resolution in the Frequency domain. Now given the signal is:
$$ x \left( t \right) = \cos \left( 2 \pi f t \right) \Rightarrow x \left( n {T}_{s} \right) = \cos \left( 2 \pi f {T}_{s} n \right) \Rightarrow x \left[ n \right] = \cos \left( 2 \pi \frac{f}{ {F}_{s} } n \right) $$
For $ {F}_{s} = 100 $ [Hz] and $ N = 100 $ (Which means $ b = 1 $) we will have delta at $ k = 2 $ and $ k = 98 $. For $ k = 98 $:
$$ 2 \pi \frac{98}{{F}_{s}} n $$
This is clearly above the Nyquist frequency ($ \frac{{F}_{s}}{2} $) and only for $ {F}_{s} = 100 $ its modulo is $ -2 $ which is correct. But in the model above, since we have arbitrary $ t $ one could think we have changing $ {F}_{s} $ which means we don't get the correct value.
This means the actual equation should be:
$$ x \left( t \right) = \frac{1}{N} \sum_{k = \left \lfloor - \frac{N}{2} \right \rfloor }^{ \left \lfloor \frac{N - 1}{2} \right \rfloor } X \left[ k \right] {e}^{ j 2 \pi k \frac{t}{T} } $$


Locating reference audio in longer audio I'm trying to solve the following problem:

Let's say I have a reference sound and a recording that I know contains at least one occurrence of the reference sound. How can I determine the positions in the audio where the sound was played? It's possible that the longer recording could have noise, and I cannot assume perfect time alignment with the sound that's playing in the recording (meaning, if I knew the content of the sound played in the recording other than the reference sound, I cannot just subtract that from the recording itself).

The method I'm experimenting with is as follows:

Grab windows the size of the reference sound from the recorded sound
For each window, compute the N dominant frequencies (by finding peak magnitudes in the FFT results)
Assign points to the dominant frequencies of the reference sound (where the most dominant gets a score of N, and the least gets a score of 1)
Compare the window's dominant frequencies to these and score it

The next few parts are just experimental, and used to narrow down results.

Do some filtering on the matches (for instance, a match with a score lower than half the possible score gets filtered out)
Look at overlapping matches and compute the cross correlation with the reference sound (using every possible time shift, prefer the maximum result). Prefer the match which requires the smaller shift
Perform that shift on the match's endpoints

The problem is that I want to also detect if part of the reference sound was played, and step 6 would not allow that if the reference sound was played two times in quick succession where one was cut off.

It seems like there has to be a solution that's not as ad-hoc as this. I'm very new to the field of DSP, so there may be something simple that I don't know about.

Any ideas?

Thanks!

It seems that a matched filter is good enough for my use case.


Best way to remove unwanted region from image I am finding the techniques to remove unwanted regions (small dots) from image. I have an image that includes object and some unwanted region (small dots- see first image). I want to remove it. Hence, I use some morphological operator example 'close' to remove. But it is not perfect. Do you have other way to remove more clear (similar second image)? You can download example image at here

Looks like you are using Matlab. Try bwareaopen(I, N), where I is the original binary image and N is the estimated size of each unwanted connected region.

You can try edit bwareaopen for more details. Basically the algorithm tries to find the size of connected regions. Connected-component labeling with union-find algorithm is expected to get you there.