How to get Instantaneous Magnitude for a Instantaneous Frequency From FFT? My understanding if a frequency in my signal doesn't line up exactly with a bin, it is smeared over a few bins to the right and left.

How can the instantaneous magnitude be determined? (do I have to worry about the windowing function?)

Frist of all you need compute the Magnitude:

windowed = framed .* hann(length(framed))
Fourier = fft(windowed)
Mag = abs(Fourier)

And you need build a Window Kernel based in your Window Function!

For a Hann Window the kernel can be:

   Wk(k) = 0.5 * (sinc(k*(M/N)) / (1 - (k*(M/N))^2))

M=Window Size
N=FFT Size

Of course you can use the same size for both, N and M !

k = Bin Number from FFT / N
sinc = Normalized Sinc, visit http://en.wikipedia.org/wiki/Sinc_function

If you have already determined the Instantaneous Frequency, you need find the corresponding offsets for each Frequency and the Instantaneous Magnitude is found when you multiply the magnitude from FFT by your Hann Window Kernel index:

 Offset = abs(InstFreqs(i) / (Fs/N) - i)
 index = floor(abs(Offset) * N) + 1
 InstMag = Mag(i) * Wk(index)


Kalman Filter's Correlation Formula I'm reading a book where the autocorrelation of white noise is expressed as:

What is the term $Q(k)$ and why is is expressed as an average value of a dot product ?

I’m guessing your book is on Kalman filters and $u(k)$ is a term in a state variable equation of a linear system.  If so, it is a vector and the components of that vector at some time index $k$ can be correlated with each other but are independent with respect to $k$.

it isn’t a dot product, it is an outer product.


AWGN : Recombining AWGN to obtain new AWGN I am a newbie to AWGN.

My professor has given me a task of taking 4 different AWGN (0 mean and SD of 1) channels (Lets say each channel has 100 samples). Then I recombine these noise samples by adding all the noise samples together. Now I divide the AWGN channel by the square root of the channel count. In this case, I have 4 samples so I divide the newly generated channel by 2. Now as I understand this new channel has same mean and SD as each of the channels that were recombined to generate this new channel. I now use this new AWGN and combine it with my signal and study a Decoders performance using signal+newly created AWGN. I will be adding the noise to a BIAWGN channel. So my signal will be +1 or -1 and noise will be added on top of that

My questions are as follows:
1. What could be the motivation of doing this. Will this help me in any way?
2. Does anyone know of any research papers that employ this method when studying a decoders performance?

My guess is your professor is trying to teach you about how noise in a dual-pol QAM system affects decoder performance. In such a system, there are 4 uncorrelated AWGN noise channels (XI,XQ,YI,YQ). It could also be that your professor is trying to teach you about linear impairments and their properties.


Decimation and filtering in the frequency domain I'm working on a Software Defined Radio project where I'd like to low-pass filter and decimate an analytical signal (IQ) sampled at 96ksps. Let's say the low-pass filter has a cutoff at 5kHz and I'd like to decimate by a factor 4 so that I have 24ksps out.

The idea is to perform the filtering using fast convolution using the overlap and save method as described in this article [pdf]:

http://www.3db-labs.com/01598092_MultibandFilterbank.pdf

I'm wondering if there are any pitfalls to my approach:

Performing an N length FFT.
Then doing an N length circular convolution (my multiplying with the FFT of my filter of length P).
Then performing an N/4 IFFT back to decimate by 4 using the N/4 center taps of the forward FFT. Since my filter is a low-pass with a cutoff at 5kHz there should be very little energy outside the N/4 center taps of the FFT, and the P - 1 samples I need to discard should also lie outside the IFFT (if my filter is not too long).

EDIT:

This specific application is on a Raspberry Pi 3. After having given this some more thought I've realized it's not as clever as I first thought. I've stared too much on 0 centered FFTs and briefly forgot this is not the case. I would have to "remove" zeros in the middle of my FFT to make it shorter and then perform the IFFT.

What I will do is to do the fast convolution with an N length FFT and N length IFFT, and decimate when I copy samples to the output buffer.

I do a lot of decimation in the frequency domain. Little details are important.

I assume you already know the basic rules for fast convolution: the FFT length N is equal to the data blocksize L plus the length of the filter impulse response M minus 1. Each operation uses L samples of new data plus M-1 samples of data from the old block.

Ensure that the impulse response of your lowpass filter is shifted to the front of your time domain buffer AND properly windowed to M samples before you take the forward FFT to get the frequency domain representation of your filter. This keeps the result from wrapping around in the time domain when you take the inverse FFT. (Remember you're actually doing circular convolution when you want linear convolution.)

Kaiser is by far my favorite window because of its tuning knob. Use a large enough value to push the sidelobes down, as they will alias into your output. I typically construct a zero-phase brick wall in the frequency domain, take the inverse transform, apply the window and shift to the front of the buffer, then take the forward FFT.

Also make sure that the main lobe of your frequency response is still essentially zero past the Nyquist limit of your decimated output. I.e., don't make the Kaiser parameter too big.

If you're also doing frequency shifting by rotating the frequency bins, remember that you have to shift by a number of bins that corresponds to one data block size. For example, if L=M-1, then you can translate by any even number of bins.

Don't worry about FFT blocksizes that aren't powers of two. Just pick a convenient size that doesn't have any large prime factors and FFTW3 will perform well. I use L=M-1=3840 to decimate 192 kHz by 4:1 to 48 kHz with a block time of 20 ms. That's a FFT blocksize of L+M-1 = 7680 = 2^9 * 3 * 5.


Confusion regarding model order and lags I have similar questions as the one asked in these posts: https://stackoverflow.com/questions/47083890/fir-filter-length-is-the-intercept-included-as-a-coefficient-matlab/47085339?noredirect=1#comment81124362_47085339
and Terminologies - lags, order in time series model

But the answer is quite different from the one in books and other online resources such as this tutorial: https://onlinecourses.science.psu.edu/stat510/node/47

If the time series model is of order 2, then according to the answers to the earlier question, there should be 3 coefficients. But, the tutorial link says differently. In the link an AR(1) model has 1 coefficient.

In a research article titled, "X. Xu and J. Guo, "A Fast System Identification Method Based on Minimum Phase Space Volume," 2012 International Conference on Cyber-Enabled Distributed Computing and Knowledge Discovery, Sanya, 2012, pp. 523-526." the Authors have doen system identification of AR(2) model. There they have used 2 coefficients to express the model.

Based on my understanding, the number of coefficients that a system has is known as the length, $L$ and the order is $q=L-1$.

I want help to confirm what is the correct representation and terminology. In general for AR and MA models,

Confusion 1) For an AR(1) system of order $q=1$, should there be 2 coefficients or 1 coefficient? Same thing for MA(1)

According to the online course link for AR(1), there is only 1 coefficient.

Confusion 2) What is correct? Considering an AR(2) model having coefficients [a1,a2,a3] then can I express the model as: x[t] = a1*x[t] + a2*x[t-1] + a3*x[t-2]+ e[t] and for MA(2) of order 2 as x[t] = a1*e[t] + a2*e[t-1] + a3*e[t-2]

where e[t] is the excitation input driving signal.

What is the correct method? Please help, I am extremely confused.

Hi: As far as what you wrote in confusion 2):

The AR expression should not have an $x[t]$ on both sides. For example, below is an AR(2):

$x[t] = a1*x[t-1] + a2*x[t-2]+ e[t]$

Similarly, the MA expression should always have a coefficient of 1 on the  $e[t]$ term. So, for example, an MA(2) would be:

$x[t] = e[t] + a1*e[t-1] + a2*e[t-2]$

As Stanley said, it is only a matter of convention but the order of an AR(p) is generally thought of as $p$ so there are $p$ AR coefficients if the order is p. Similarly, for an MA(q) there are q MA coefficients so the order is considered to be q ( the $e[t]$ is not included in the order for an MA because there is no coefficient associated with $e[t]$).

Note that, for both MA and AR models,  the number of parameters to be estimated is not always equal to the order of the model because there may be non zero mean terms or trend terms. For example, below is AR(2) with unknown mean $\mu$:

$x[t] = \mu + a1*x[t-1] + a2*x[t-2]+ e[t]$.

Note that there are also ARMA models where the number of coefficients is then p+q. For example, an ARMA(2,1) is:

$x[t] = a1*x[t-1] + a2*x[t-2]+ b1*e_[t-1] + e[t]$.

But usually the term "order" is not used for ARMA models.

I hope this makes it clear.


Speech recognition using MFCC and DTW(Dynamic Time Warping)? I know basics about DSP, and now trying to complete a project on Speech recognition. So I read as many resources as I found, and got some ideas. They are like:

take a sample sound and break it into small frames

calculate MFCC for each frames( 13 of them)

use DTW to match the testing and template pattern.
Now, here are the things I am having trouble with-

13 MFCC's for n frames leaves me with 13*n MFCC's. How exactly do I use DTW here?

DTW gives me a distance between two patterns. How do I determine if they match or not?

At the very beginning let me warn you that DTW approach is suitable only for spoken word recognition. Nonetheless it is interesting as a basic exercise.

I assume that you have a database of training files (templates) and you already extracted MFCC's for those. Each training file contains one utterance for a given word class, i.e. "Hello", "Bye". I also imagine that each class has more than one example (plausible).

Here is what you should do during the recognition stage:

Once you have the recording for the whole word, you calculate the MFCC's for each frame. You want the recording of the whole word, not just part of it.
Now you calculate the distance between your recording and each of the templates in your database. In case of DTW you will be calculating the cost between each of 13 dimensional frames (simple distance metric/norm, i.e. Euclidean, Manhattan, etc.). Once the DTW algorithm is finished, you will end up with the distance value (upper right matrix entry as can be seen below) between your test sample and each of the templates.
The last step is to make a decision: to which class your test sample actually belongs to? You can either do it by picking the class of template with the minimum DTW distance. But what is even better, you can use the k-Nearest Neighbours for that.


$y[n] = x[n-1]$ is an all pass filter I'd like to know if I understand correctly why $y[n] = x[n-1]$ is an all pass filter.

By definition of the filter, we get the output the same of the input only delayed by one unit of time. Therefore, all frequencies that entered must pass (eventually).

Therefore, it's an all-pass filter.

Is that a correct reasoning?

Thanks

Your notation is a bit odd, as you are mixing the time-domain with the $z$-domain in the equation you wrote. If you refer to the system

$$y[n]=x[n-1] \quad \mathrm{or} \quad Y(z)=z^{-1}X(z)$$

then you are correct. To put it mathematically, the transfer function of the filter would be:

$$H(z)=\frac{Y(z)}{X(z)}=\frac{1}{z}$$

The frequency response is then

$$H(z) \Bigg|_{z=e^{j\omega}} =\frac{1}{e^{j\omega}}$$

You can see that the magnitude of $H(e^{j\omega})$ is $1$ for all frequencies. Therefore, it is an all-pass system.


2D Fourier downsampling How can one correctly downsample a Fourier 2D matrix of complex numbers ? More precisely, if I have a 2D DFT of an image:

x=phantom
X=fft2(x)

How should I do to in the frequency domain to obtain:

XHalf=fft2(imresize(x,0.5))

(And imshow(ifft2(XHalf)) shows the downsampled image of course).

Of course, doing imresize() in the frequency domain does not work.

If you do an imresize(X,0.5), after doing an ifft2(imresize(X,0.5)) the result is nonesense.

I tried simply cropping the Fourier "image": it does give something recognizable but it is clearly not simply that.

For a 2D FFT (fft2), first do an "fftshift" to center the low frequency portion of the 2D FFT in the center of the matrix. Select this inner region and then do an inverse 2D FFT (ifft2) of the fftshifted result.

The format of the result also needs to be 8-bit unsigned.

The image below depicts the low frequency portion of the frequency domain of the image in green centered on the origin after using the fftshift command (my apologies to our B/G color blind users!). After executing fft2, the pixel corresponding with the lowest frequency will be in the upper left hand corner. fftshift will swap quadrant 2 and quadrant 4, and swap quadrant 1 with quadrant 3 resulting in the lowest frequency pixel to be at the origin. This simplifies the selection of the lowest frequency quartile of the image and after processing the ifft2 of the ifftshifted result will provide the half dimensioned image as desired.

The devil is in the details to properly shift and select the image properly, so see the MATLAB/ Octave code below that properly does the proper selection as well as the scaling.

d = 2;   # decimation rate
dim = size(image);
m = floor(dim/2)+1;
n = floor(m/d);
fimage = fftshift(fft2(image));
# select middle
fimageh = fimage(m(1)-n(1):m(1)+n(1),m(2)-n(2):m(2)+n(2),:);
# ifft
image2 = uint8(real(ifft2(ifftshift(fimageh/d^2))));
imshow(image2);

Note that the code includes taking the real part, but if done correctly all the imaginary terms should be insignificant. However uint8 will not operate on a complex numbers.

Test image before (800x1200):

Test image after for d=2 (401x601):

image source: https://twitter.com/quirktree/status/1079911361122131968


Scaling Function in Wavelet Analysis I was watching a YouTube tutorial on wavelets https://www.youtube.com/watch?v=wgRgodvU_Ms, where the presenter showed that if we multiply a Gaussian with a sine wave (and follow some other conditions), we get an example of a Morlet wavelet.

In the case of Haar, the scaling function (left) and the wavelet (right) have a different shape. How does one design a corresponding scaling function once a wavelet is designed? For example, a Mexican hat wavelet is a second derivative of a Gaussian with its signs reversed. How does the corresponding scaling function look like.?

A pointer to an introductory reference would be appreciated as well.

There is a huge gap between the formalism and construction for the continuous wavelet analysis (Mexican hat or Morlet wavelet, the latter being only approximate or complex Gaussian wavelet) and the discrete  wavelet transforms (DWT). Let us make a long story short.

Continuous wavelet transforms (CWT) hinge on a continuous wavelet $\psi$, for which for any (nice) function $f(x)$ can be represented as (or analyzed) a weighted ($c_{a,b}$) sum of translates and dilations  $\psi\left(\frac{x-b}{a}\right)$. Whether the synthesis formula is equal to $f(x)$ requires $\psi$ to satisfy a moderately mild admissibility condition.

Then a question arose: how can we discretize all of the above to get tractable algorithms? In other words, are there $a_0$ and $b_0$ such that projections onto enumerable series

$$\psi\left(\frac{x-k b_0 a_0^j}{a_0^j}\right)$$

with $j$ and $k$ integers can still be perfect? The answer is a big yes, provided that $a_0b_0 < C_\psi$, a constant depending on $\psi$ (this is a king of Pauli-Heisenberg-Gabor inequality, result by Ingrid Daubechies, 1984). This becomes a wavelet frame: a redundant set of vectors that exactly recovers signals.

Now, with DWT, we want a base, i.e. a non-redundant frame. Something that keeps the size of the initial data. From the above harmonic analysis techniques, if you choose an admissible wavelet, there is very little chance that a basis could be derived from it. Haar or Frankling bases were known examples, and Yves Meyer in 1985 derived a non-trivial example, the Meyer basis.

To build discrete wavelet bases, it was discovered that a more efficient method was needed.
On a discrete space, as a given resolution, you except that there exists an orthogonal basis from an orthogonal template of a translated scaling function or father wavelet. Then, you embedded different resolutions,  and you get a 2-scale linear  equation linking  $\psi(t)$ and $\psi(t/2)$. Using orthogonal projections, you can build another  2-scale linear equation linking a mother wavelet at different dilations. Those equaltions tell you than you can implement the above with convolution and subsampling. Consequently, there are polynomial relations between scale. Solving such nonlinear equations yield a series of weights/linear coefficients providing a hierarchy of discrete wavelets that, plugged into the formalism, are genuine CWT. However, in most cases, they don't have nice close formulae like Gaussians or Morlet. The latter in turn have no discrete counterpart.

In Stéphane Mallat "A wavelet tour", chapter one on sparse representations says:

  Continuous versus Discrete and Finite Originally, many signal processing engineers were wondering what is the point of considering
  wavelets and signals as functions,since all computations are performed
  over discrete signals with conjugate mirror filters.Why bother with
  the convergence of infinite convolution cascades if in practice we
  only compute a finite number of convolutions? Answering these
  important questions is necessary in order to understand why this book
  alternates between theorems on continuous time functions and discrete
  algorithms applied to finite sequences. A short answer would be
  "simplicity." In L2(R), a wavelet basis is constructed by dilating and
  translating a single function $\psi$. Several important theorems
  relate the amplitude of wavelet coefficients to the local regularity
  of the signal $f$. Dilations are not defined over discrete sequences,
  and discrete wavelet bases are therefore more complex to describe.The
  regularity of a discrete sequence is not well defined either, which
  makes it more difficult to interpret the amplitude of wavelet
  coefficients. A theory of continuous-time functions gives asymptotic
  results for discrete sequences with sampling intervals decreasing to
  zero.This theory is useful because these asymptotic results are
  precise enough to understand the behavior of discrete algorithms. But
  continuous time or space models are not sufficient for elaborating
  discrete signal-processing algorithms.The transition between
  continuous and discrete signals must be done with great care to
  maintain important properties such as orthogonality. Restricting the
  constructions to finite discrete signals adds another layer of
  complexity because of border problems. How these border issues affect
  numerical implementations is carefully addressed once the properties
  of the bases are thoroughly understood.


Solve for minimum and maximum value of dct coefficient Solve for minimum and maximum value of dct coefficient

The max value for a DCT coefficient with this formula would be when $k_1 = 0$ and $k_2 = 0$.  Using this give us the $\cos(0) = 1$ so you simply sum up the pixel values.  Since the pixel values are $0-255$, then the max value is $64*255 = 16320$.
The min value can be computed by putting a pixel value of $255$ whenever the product of the cosines is negative, and a zero when it is positive.  You would have to iterate through the $k_1$ and $k_2$ values to find the smallest value.


Why correlation is not associative in the context of image processing? I can see most of the places it says "convolution is associative, while correlation, in general, is not".

Denote $*$ convolution operator, let's say you have an image $f$, which you need to convolve with $g$ and then with $h$ :
$$f∗g∗h=f∗(g∗h)$$

But when I see the difference between convolution and correlation, correlation is equal to convolution but after flipping the kernel or window.

Can somebody explain with an example why correlation is not associative?

Subtraction is also not associative: (a-b)-c does not equal a-(b-c) in general. The order of the flipping of the sign matters, just like it does for the flipping of the kernel.


Accessing Maximum Value from a Singular Value Decomposed Matrix I have read few image enhancement papers, where an image is enhanced in transformed domain using either DCT (discrete cosine transform) or DWT (discrete wavelet transform) with SVD (singular value decomposition). An example is given in  Enhancement of low contrast satellite images using discrete cosine transform and singular value decomposition, 2011.

In these works, once the image matrix is factored into three matrices (say $U$, $S$, and $V$) using the SVD, then only the maximum value of matrix S is taken for processing. Can anyone please tell that why only the  maximum value is taken for processing?

The SVD Decompose the image into the (One way to look at it) many matrices.
For instance, given an Image $ I $ its SVD is given by:

$$ I = U S {V}^{T} = \sum_{i=1}^{\textrm{rank}(I)} {s}_{i} {u}_{i} {v}_{i}^{T} $$

As you can see, the image is decomposed into linear combination of matrices (Given by $ {u}_{i} {v}_{i}^{T} $) where the weight of the $ i $ -th matrix is given by $ {s}_{i} $ which is the $ i $-th Singular Value (Usually they are organized in descending order).

The dominant features are usually goes with the large Singular Values (Moreover, the smallest one usually describe noise).
I guess methods want to focus on the large dominant features and avoid dealing with noise hence it makes sense to deal only with the dominant matrix in this decomposition.

By the way, I'd experiment with trying to work on the next ones as well.
In addition, on the recompositon process I'd neglect the last few to have "Simple Denoising".

One should be aware that the SVD of the whole image might be time and memory consuming when dealing with large images.

Moreover, it makes sense to stack few image patches into large matrix and then apply the SVD.
This might extract what is common to all patches (Linear Dimensionality Reduction).


Is there a relatively easy way to detect likely real-time peaks in discrete-time data? Let's say I have a set of data over time, t:

[0, 4, 6, 7, 7, 6, 4, 0]

It seems likely that this data would peak at t=3.5.

Is there a well-known algorithm for calculating this sort of peak?

The most straight-forward solution would be to interpolate the discrete values around the maximum you're interested in, and compute the location of the maximum of the interpolating function. I'd start by using a quadratic or cubic polynomial. This is a simple example of using an interpolating quadratic polynomial.


Computing real signal with minimum absolute values from even magnitude spectrum I want to derive a real audio signal from an arbitrary even magnitude spectrum.

The phase spectrum affects the values of the signal in the time domain; for example, a phase of 0 for all frequencies amounts to summing cosines in the time domain, which will result in a large peak at t = 0. I want to minimize the absolute values of the signal in the time domain without scaling the magnitude spectrum

I want this so that the audio signal can be multiplied by the largest scalar value possible without clipping. Is there a way to compute the phase spectrum that will do that? If there is a stupidly simply solution, I apologize in advance.

Turns out this problem is known as crest factor minimization. I'm using an algorithm by Yang et al. (Phyiol. Meas. 36, 2015). Crest factor minimization algorithms typically are not proven to work on all possible magnitude spectra and only rarely proven to work on magnitude spectra with certain characteristics. Also, some iterative algorithms like Yang et al.'s have been shown in later publications to get stuck at a local minimum, not the global minimum. However, Yang et al.'s algorithm is satisfactory for the magnitude spectra I've tried so far.


Obtaining normalized matrix for the Haar Wavelet Transform I've been reading this article: http://aix1.uottawa.ca/~jkhoury/haar.htm which explains the Haar Wavelet Transform.

At a certain point, the author says:

  ... Since the transformation matrix $W$ is the product of three other matrices, one can normalize $W$ by normalizing each of the three matrices. The normalized version of $W$ is ...

and provides final matrix.

But he doesn't provide the exact procedure for matrix normalization he used. Probably someone with a deep understanding of the Haar Wavelet Transform can explain me, how matrices were normalized so that final result was obtained.

I've been searching the Internet for matrix normalization and couldn't find anything suitable. Also contacted author of article but didn't get any response so far.

After a few quick calculations, it seems to me that the trouble comes from poor notations for the root in your reference. If you read, in the final normalized matrix, $\sqrt{8/64}$ and $\sqrt{2/4}$ instead of  $\sqrt{8}/64$ and $\sqrt{2}/4$ (along with the $\pm$ signs), then  the final result is correct.

The matrices $V_i$ are orthogonal. To normalize them, you only need to multiply them by a diagonal matrix $D_i$ made from the inverse of the norm of each column.

For $V_1$, the norms are all equal to $\sqrt{1/2^2+1/2^2}=\sqrt{1/2}$. Hence, you can multiply $V_1$ by the diagonal matrix
$$D_1= \operatorname{Diag} \{\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2}\}\,.$$
and get an orthonormal matrix. Similarly, you get that:
$$D_2= \operatorname{Diag} \{\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,1,\,1,\,1,\,1\}\,,$$
and
$$D_3= \operatorname{Diag} \{\sqrt{2},\,\sqrt{2},\,1,\,1,\,1,\,1,\,1,\,1\}\,.$$

Finally, because of commutativity, you get a matrix $D=D_1 D_2 D_3$:

$$D= \operatorname{Diag} \{2\sqrt{2},\,2\sqrt{2},\,2,\,2,\,\sqrt{2},\,\sqrt{2},\,\sqrt{2},\,\sqrt{2}\}\,.$$

If you now multiply the matrix $W$ (the one with the $1/8$)

with that matrix $D$, you get a final matrix $W$:

provided the notation interpretation for the square root. For instance, you can recover $1/8\times 2\sqrt{2} = \frac{\sqrt{8}}{\sqrt{8^2}} =\sqrt{8/64} $.


Understanding Soft Thresholding Operator I want to understand what is soft thresholding operator? Why we use this operator? I came across this term while I was reading the paper - A New Algorithm Based on Linearized Bregman Iteration with Generalized Inverse for Compressed Sensing (at page no 4.).

Any reference that can explain this will also help me.

Thank you for your help.

The soft-thresholding function finds the minimizer of an objective function  that involves data fitting in an $\ell_2$ sense as well as minimization of the $\ell_1$ norm (i.e. absolute value). The Lecture Notes - Penalty and Shrinkage Functions for Sparse Signal Processing gives a good discussion of how the soft threshold function is derived.


Training custom classifiers using trainCascadeObjectDetector Is there any good documentation about training your own classifiers(with examples included) different from this one provided by MathWorks.
Thanks in Advance

Check out this tutorial on how to use trainCascadeObjectDetector. It describes common pitfalls, and includes examples.

There is also a graphical tool for labeling training images on MATLAB Central, which makes preparing your training data much less painful.

By the way, the link you posted is actually the documentation for vision.CascadeObjectDetector, which uses the classifier produced by trainCascadeObjectDetector. The doc page for trainCascadeObjectDetector is here.


How do I get a "blurry line" in gimp? This is a close up of a Pokémon (Teddiursa):

How do I achieve this same effect in GIMP- i.e. a line that is kind of blurry on the edges? I'm sure there's a term for this, and I'm sure there is a really easy solution- I'm just not sure what to google.

Note: I may have asked this on the wrong site (I just searched image in the bar and this site came up). If I did, feel free to redirect me to where I should ask!

Thanks in advance!

It's called anti-aliasing and happens when you low-pass filter your signal (image) and then decimate it (scale it down). Hence, the way to get this is take a large picture with sharp edges, and image->scale image, using anything but "None" as method. Alternatively, blur filter (==low pass) or selective blur will also do.


wavelength / RGB valve per pixel Preamble: I want to stress that I don't know anything about photo-editing or signal processing. But I am a composer who happens to be working on a project in which sound / light overlap.

Question: Am I correct in understanding that a pixel of a given image contains a single RGB value?

Many thanks.

pixel data of color-images contain three (bytes of) intensity values denoted as R,G and B. Most common they are 8-bits but other formats are possible too.

Each intensity value at a pixel point represents the amount of light captured by a corresponding sensor at that point which is most sensitive to the given range of wavelengths (color); for example a RED value means that the associated RED sensor at that pixel captured that much of light in the wavelengths about $550$ nm to $750$ nm.

In a practical setting, the distinct wavelength captures can be realized by three different optical filters that separate the light spectrum into three regions associated with Red, Green and Blue. Then the incoming light (after being filtered) will be captured by three same type of sensors which are sensitive to the whole optical range (+ more possibly)


Average identical sweeps or lowpass filter sweeps first? There are 5-10 replicate sweeps of measurements that are meant to be averaged. I also typically lowpass filter the signal at 2000Hz. The target signals are in the 0-100Hz range, and have, for lack of a better term, varying variability (the 5-10 sweeps are usually very similar to each other, but sometimes they are pretty different, but averaging out the variability is fine since it's random).
My question is whether I should (1) lowpass filter the 5-10 sweeps individually before averaging them, or (2) average the 5-10 sweeps and then lowpass filter the resulting mean.
Is there a general rule that applies to this? I am new to signal processing, but I'm aware that order of operations matters.

If your are using a linear time invariant lowpass filter (such as a Finite Impulse Response (FIR) or an Infinite Impulse Response (IIR) filter), and have negligible rounding errors in each operations then doing the lowpass filtering before or after the averaging will yield the same results.
Note that from an implementation perspective, doing the averaging first would likely require fewer computational resources.


Extracting common signal from 4 sets of observations I am working on a signal processing assignment where I need to find out one common time domain signal from 4 observation. The math is like this:

\begin{array}{lcl} y_1(t) & = & a_1(t)*x(t) + b_1(t) \\
                   y_2(t)& = &  a_1(t)*x(t) + b_1(t) \\
                   . & & \\
                   . & & \\
                   y_n(t) & = & a_n(t)*x(t) + b_n(t)\end{array}

where:

$n$ is up to 4 in my case, and $x(t)$ is my signal of interest which I want to extract from observed signals $y_n(t)$.

What is the general direction I should take in doing background study for this kind of application?

Thanks,
K

Because you have a set of linear equations, why don't you set up a linear system (for each data point) and solve the system (using SVD for example) to obtain the solution?

Write in the the matrix form such as $Y=\beta X$ where $\beta$ are known and $X$ is unknown.


Are there standards/references for determining spatial scale in photographic images? EDIT, 12/12/20:
Images below are of the radial sinusoid pattern. Left side is the "unrotated" version. Right side is what the a photograph of the pattern would look like if the disc was rotated into the depth dimension, 45 degrees about the x/horizontal axis, and 20 degrees about the y/vertical axis. Note that the diameter of the disc on the left is equal to the length of the long axis of the disc (ellipse) on the right. And thus the spatial frequency of the sinusoid along the long axis is equal to that of the unrotated disc.

END EDIT
I'm tasked with estimating the spatial (pixels/inch) resolution at the location (depth) of an object in a digital photograph. I can create a small "standard" reference signal (akin to a QR code or bar code) that will be clipped to the object, such that I can assess the resolution via image processing, Fourier analysis, etc.

Does such a standard reference exists? Along with optimized algorithms for estimating the spatial resolution via image processing?
If not, I'm looking for ideas for patterns that might serve as such a reference, and the associated estimation methods.

Some caveats:

The object and the visual standard are not fixed in space, and may rotate (into the depth dimension), to some degree.
The visual standard needs to be limited to ~3 inches in diameter (or a ~3x3 inch square).

Given the above caveats, my early thinking has been to use a circular reference (a "disc"), displaying concentric rings of a fixed spatial frequency--e.g., a sinusoid that is radial symmetric. If the disc rotates in any direction, the concentric circles will look like ellipses in the photographic image, but the sinusoidal frequency along the long axis of the ellipses will remain constant. I thought a 2D FFT would allow me to obtain a good frequency estimate, but based on some testing, the circular nature of the pattern results in "artifacts" in the frequency domain, such that the estimate is not sufficiently accurate. Perhaps the concentric rings are an adequate reference/standard pattern, but a better frequency estimation method is needed (I've also tried direct fitting via gradient descent). Or maybe I need some other pattern altogether?
Suggestions welcome!

One possible answer to my own question is provided by the photographic scale reference detailed on this website.
https://smallpond.ca/jim/scale/

It doesn't quite provide the solution I was looking for, buy may help spur ideas. I really need something that can accommodate a good scale estimate even if the reference is rotated in space (about any of the three axes), and is of low pixel resolution.


What is the definition of a "complementary filter"? I'm writing my MSE thesis report on an inertial navigation problem I've been working on. In my work, I use a complementary filter to track the orientation of a device. Now, since my fellow students are likely not familiar with complementary filters or Kalman filters, I wish to give a brief description of what a "complementary filter" is, in general.

A definition, if you will.

So, what defines a complementary filter? Is it any filter on the form z = a * x + (1 - a) * y, where x and y are separate measurements of a single quantity?

It's fine to define it in terms of Kalman filters if that's appropriate.

Usually, a complementary filter (like a complementary function) complements another filter.  The two filters that are complementary to each other add to one.  Or, at least, add to an all-pass filter (which is what Linkwitz-Riley crossovers do.

so either

$$ H(f) + G(f) = 1 $$

or

$$ H(f) + G(f) = A(f) $$

for $H(f)$ and $G(f)$ being complements of each other and $|A(f)|=1$ is an APF.


What is difference between "Direct form" and "Normal Direct form"? What is difference between "Direct form" and "Normal Direct form"?
I was reading about filters implementations(structures) and found out both words "direct form"and "normal direct form" especially in proakis,2nd edition as shown highlighted in attached snapshot

The modifier normal is meant to distinguish the form from the transposed form. So there is a normal direct form I, a transposed direct form I, and similarly for the direct form II. As far as I know, this terminology is not universal, and I've only come across it in books (co-)written by Manolakis or Ingle.


How can i compute perimeter and area of a skeletonized circle? In opencv if i use findcontours function, i get contours and i can compute area, perimeter, convex hull.
if i use a Skeletonized image of a circle how can i compute perimeter,area, convex hull of the skeletonized circle?

UPDATE

image:

My skeletonize function:

import numpy as np
import cv2
import scipy.ndimage.morphology as m

def skeletonize(img):
    h1 = np.array([[0, 0, 0],[0, 1, 0],[1, 1, 1]])
    m1 = np.array([[1, 1, 1],[0, 0, 0],[0, 0, 0]])
    h2 = np.array([[0, 0, 0],[1, 1, 0],[0, 1, 0]])
    m2 = np.array([[0, 1, 1],[0, 0, 1],[0, 0, 0]])
    hit_list = []
    miss_list = []
    for k in range(4):
        hit_list.append(np.rot90(h1, k))
        hit_list.append(np.rot90(h2, k))
        miss_list.append(np.rot90(m1, k))
        miss_list.append(np.rot90(m2, k))
    img = img.copy()
    while True:
        last = img
        for hit, miss in zip(hit_list, miss_list):
            hm = m.binary_hit_or_miss(img, hit, miss)
            img = np.logical_and(img, np.logical_not(hm))
        if np.all(img == last):
            break
    return img

nomeimg = "JPEG/e(5).jpg"
img = cv2.imread(nomeimg,0)
ret,img = cv2.threshold(img,127,255,0)
element = cv2.getStructuringElement(cv2.MORPH_CROSS,(3,3))
img = 255 - img
img = cv2.dilate(img, element, iterations=7)
img = cv2.erode(img, element)

#element = cv2.getStructuringElement(cv2.MORPH_CROSS,(6,6))

skel = skeletonize(img)
skel = skel.astype(np.uint8)*255

It is not clear if the skeletonized image can be reversed to a contour image. Skeletonization is often done by calculating the distance transform and then thresholding the result or performing a non-maximum suppression. At this point you have already lost too much information and the features, you are asking for cannot be calculated.

If you have the distance transform, you can approximate those values:

You get the radius ($r$) of the circle from the height of peak which represents the circle (the skeleton of an ideal circle would be a dot)
Perimeter = $2\pi r$
Area = $\pi r^2$
Convex hull $\rightarrow$ draw a circle with radius $r$ at the position where the peak is


Using ROC to find stability of system in specific example I've started learning about finding the ROC from the transfer function, but I'm confused about an example.

$$H(Z) = \frac{2Z + 1}{Z^2 + Z - \frac{5}{16}}$$

I understand the poles lie at $z = \frac{-5}{4}$and $z = \frac{1}{4}$, and I assumed the system was non-causal because even though it has more poles than zeros, the unit circle does not lie to the exterior of the outermost pole(if I'm right), but because the ROC includes the unit circle, the system was stable. However, the solution to the problem states the system is both unstable and non-casual so I was hoping someone could explain this to me? I've only seen examples giving the ROC when both poles have been of the form $(z-a)$ however here one is $(z + b)$, so I thought the ROC was $\frac{5}{4} > |Z| > \frac{1}{4}$? I imagine this is where I have went wrong.

Thanks in advance.

because even though it has more poles than zeros

WRONG. The Z-transform transfer function will always have equal number of poles and zeros. Your poles are at $z = -5/4$ and $z = 1/4$. Zeros are at $z = -1/2$ and $z = \infty$. Since ROC will always be concentric circle region without including poles, there are 3 possible ROC for the given transfer function.
$$
ROC_1 = 5/4 \lt |z| \lt \infty
$$
$$
ROC_2 = 1/4 \lt |z| \lt 5/4
$$
$$
ROC_3 = 0 \le |z| \lt 1/4
$$
$ROC_1$ represents an unstable but causal system (since it can include $z \rightarrow \infty$ so terms like $z^{-k} \,\,with\,\, k \ge 0$ will only be present (only negative or zero powers of $z$). And because only negative powers of $z$ are present, system is non-zero for $n \ge 0$. So it is causal. Example - $H(z) = a_0 + a_1z^{-1}+a_2z^{-2}+a_3z^{-3}+a_4z^{-4}+...$

$ROC_2$ represents stable but non-causal system because it includes unit circle $|z| = 1$ but does not include $z \rightarrow \infty$. Terms with $z$ raised to positive power can be there along with negative power of $z$. So you can rule out $z=0$ and $z \rightarrow \infty$.Example - $H(z) = a_0 + a_1z^{-1}+a_2z^{-2}+a_3z^{-3}+a_4z^{-4}+a_{-1}z^{1}+a_{-2}z^{2}+a_{-3}z^{3}$

$ROC_3$ represents an unstable and non-causal system because it does not include unit circle and terms with only positive power of $z$ are present (because ROC includes $z=0$ and does not include $z \rightarrow \infty$. Example - $H(z) = a_0 + a_{-1}z^{1}+a_{-2}z^{2}+a_{-3}z^{3}+...$

In your case, the ROC might be $ROC_3$.


Disadvantages of wavelet transform I have a question related to wavelet transform: we know that while the Fourier transform is  good for a spectral analysis or which frequency components occurred in signal, it will not give information about at which time it happens. That's why the wavelet transform is suitable for the time-frequency analysis. It is also good for signal denoising, but of course it has some disadvantages.

So I would like to know what are main advantages of the wavelet transform? Is it good for spectral estimation; like finding amplitudes, frequencies and phases, or it just helps us to find discontinuous and irregularities of a signal?

Thanks in advance

If you consider the whole set of potential wavelet transforms, then you have a lot of flexibility.

For instance, should you use 1D continuous complex wavelet transforms, by analyzing the modulus and the phase of the scalogram, and provided you use well-chosen wavelets (potentially different for the analysis and the synthesis), and a proper discretization, you can:

find discontinuities and irregularities of a signal and its derivatives
find break point location by wavelet ridge extrapolation
denoise
perform matched filtering based on templates (with complex continuous or discrete dual-tree wavelet frames)
analyse (multi-)fractalty
analyse frequencies (with Gabor wavelets for instance)

Due to the redundancy, and the quantity of available wavelets (not the same is best for different purposes), they could  appear a little less efficient for the analysis of pure stationary and harmonics signals, for which Fourier is better suited.

The main drawbacks are:

for fine analysis, it becomes computationaly intensive
its discretization, the discrete wavelet transform (comp. efficient), is less efficient and natural
it  take some energy to invest in wavelets to become able to choose the proper ones for a specific purpose, and to implement it correcly.


FFT High Frequency Humps I am using Matlab to take the FFT of a series of pulses. Can anyone explain why there are multiple humps at high frequencies rather than a smooth roll-off?

If I read the plot correctly, I see tones repeating every 1 MHz, 3MHz, 5MHz, etc under the envelope of a $Sinc^2$ function with nulls every 200 MHz and an overall sampling rate of 2 GHz (You further clarified 10 GHz, so the full digital spectrum must not be shown?). That envelope can be the result of a 10 ns triangular pulse, and given the repetition rate, consistent with the plus minus pattern you described repeating every 1MHz.

Please include a zoom-in of the time domain of your waveform, specifically over a 50 ns time span showing the discrete samples. This will further clarify if this theory is correct of if you are also being affected by a zero-order hold.

Here is a plot of showing the frequency response of [1 1 1 1 1 1 1 1 1 1] convolved with itself, thus creating a triangular pulse 10 ns wide at a 2 GHz sampling rate), which is a $Sinc^2$.  The same result can be achieved at the 10 GHz sampling rate; being simply discrete samples of a 5 ns rectangular pulse - which result in a Sinc function in frequency with nulls at 200MHz spacing; since we see from the frequency roll-off in your plot what would be a $Sinc^2$ function, then two rectangular pulses each 5 ns wide convovlved in time (resulting in 10 ns triangular pulses) would then give the $Sinc^2$ function achieved:

Generally, you can determine the Fourier Transform of any repeating pattern by using the Fourier Transform of the base pattern as the envelope (in your case the triangle waveform which results in the $Sinc^2$ response as the envelope) and then within that envelope, due to the signal repetition, the only non-zero frequency components that can only exist will be at integer multiples of the repetition rate of the base pattern (I explain this in detail in another response and will update this with the link to that if I can find it). In your case it is further complicated by a positive triangle wave followed by a a negative triangle wave as your "base pattern" that repeats at the 1 MHz rate. To get the resulting envelope in this case, this can be easily conceptualized by first passing your triangle wave through a comb filter response as shown in the plot below (The comb filter is a delay of 500 ns and a subtraction resulting in nulls at DC, 2 MHz, 4 MHz etc).  Thus the original $Sinc^2$ envelope now has nulls at DC, 2 MHz, 4 MHz. When you repeat this base pattern at the 1 MHz rate, this would result in tones every 1 MHz, except the even ones are nulled due to the effect of the comb filter.

In the plot below, this process is shown both in the time domain and the associated response in the frequency domain for each portion of this process as previously described.

So with all that said, as long as your triangle pulses are 10 ns and not the 50 ns you first stated, the FFT you show is how it should look and not a problem with your settings. However you may want to consider windowing your signal prior to taking the FFT to reduce the sidelobes shown on each of the tones (that are an artifact of the FFT when windowing is not used). This will reduce the width of each "tone" but substantially reduce the sidelobes.


Maximum Ratio Combining vs. Selection Combining In literature, I am seeing that MRC almost always works better than SC. Also, in my simulations for Rayleigh and Rician fading channel, I am getting the same results. Could anyone please explain what are the reasons behind this? I tried to read some papers but did not get a clear idea.

Very quick explanation: SC just throws away $N-1$ of $N$ observations, whereas MRC combines all $N$.

MRC does to a signal in spatial domain what a matched filter does in frequency domain: it maximizes the inner product of signal vector and weights. The exact same derivation as for the matched filter's optimal performance under AWGN apply to the MRC under the same assumption.

You can also do the same consideration on a basic stochastic level.

Assume the combiner has correct knowledge of each observation's SNR.

Then, SC picks the best one. Derivation is simple: if channel realizations are Rayleigh distributed, then the best of these channels has $\sum_{n=1}^N \frac1N$ increased SNR in expectation, $N$ being the number of concurrent observations.

For MRC, you literally get the weighted average of the received signals, so your overall SNR literally is $\sum_{n=1}^N \frac {S_k}{N_k}$, which yields a diversity gain of $N$ (the maximum possible).

Regarding your "I've read papers": the 1959 Brennan paper is the only one you need to read for this :) But, honestly, this is knowledge that's 60 years old, so you don't have to read the original paper, and should probably rather read a good modern textbook or lecture slides explaining combining methods, as these tend to didactically process the paper's content to make it more accessible. A really quick google for "diversity gain MRC" instantly yielded a link to these lecture notes, and you could have found them, too! They're easy to read, I find after quickly scanning them.


Assessing relative wireless channel quality using Capacity and Condition Number CDFs I have Channel capacity analysis figures for two wireless channels show below:

My question:

While ergodic capacity in the first subplots are very clear in conveying the channel capacity, How to make an assessment of channel capacity based on P(Capacity < x) in second subplots and P(Cond. Num < x) in the third subplots for each of the above channels, individually and relatively?

What the CDF of the capacity tells you is how it distributes, which allows you to say something about the quantiles. Ergodic capacity tells you: what capacity will we see on average? Quantiles are more relevant in practice: what capacity can we guarantee with 95% certainty? 99%? 99.9%? For this, you need a CDF and more often than not, you'll see a heavy-tailed distribution there.

Now, what you would like your system to do is a CDF that's very steep (yours are looking pretty good), which implies a well localized PDF, i.e., a low variance. It means that you'll have the 20-something bsp/Hz not just on average but with a pretty high probability. If the CDF were less steep you might still have 20 on average but with variations between 10 and 30, which is not that useful in practice if there is a certain data rate you want to guarantee.

For the condition number: this is an important aspect in MIMO channels if you want to maximize data rate by applying some form of spatial multiplexing, i.e., sending different data streams from your transmit antennas either towards different users (like in a multi-user MISO oder MIMO system) or towards your different RX antennas in a P2P system. For this to work, the MIMO system needs to offer enough spatial degrees of freedom. It needs to be full rank for it to work properly, and it will work even better of all the spatial eigenmodes are equally strong, which means the channel's singular values being almost equal (think of the singular values as the "strength" of individual spatial eigenmodes). This ideal case implies a condition number close to 1. The further it is from one, the larger the eigenvalue spread, which means some spatial eigenmodes will be weaker and thus be able to sustain only a lower data rate.


Vein extraction from this image I have this image:

What methods should I use to extract the veins from this image.
Thank you.

I would like to direct you to 3 references:

  C. Steger: “Extracting Curvilinear Structures: A Differential
  Geometric Approach”. In B. Buxton, R. Cipolla, eds., “Fourth European
  Conference on Computer Vision”, Lecture Notes in Computer Science,
  Volume 1064, Springer Verlag, pp. 630-641, 1996.

  C. Steger: “Extraction of Curved Lines from Images”. In “13th
  International Conference on Pattern Recognition”, Volume II, pp.
  251-255, 1996.

  C. Steger: “An Unbiased Detector of Curvilinear Structures”. IEEE
  Transactions on Pattern Analysis and Machine Intelligence, vol. 20,
  no. 2, pp. 113-125, 1998.

In these works, Steger develops a subpixel curvilinear structure extraction algorithm, which is found to work very well for such images. I have also run the algorithm on your images and here you go:

Lines visualized on original image

Lines only

I set the line thickness to 2 for better visualization, normal contours are at sub-pixel level.


Meaning of frequency and bandwidth of a signal, despite the fact that we do not know the signal First of all, I am completely new to the domain of signal processing.
As far as I know, a signal can be represented with an infinite integral of infinitesimal complex exponentials, which is known as a Fourier Transform. From calculating the Fourier Transform of a signal, we can find the bandwidth of it as well.
We tend to say that voice signals have a bandwidth starting from 50Hz and ending approximately at 10kHz. At the same manner, we say that radiowaves and microwaves lie on a specific bandwidth. How is it possible to know the bandwidth of a signal, if we don't know the signal itself? I mean, when speaking, the output signal produced by my mouth is always different as time proceeds, but we know that the bandwidth of that signal is always at the aforementioned frequency interval. How is that possible to know?

Some signals are generated by processes that are physically capable of generating frequencies only in a specific frequency range. For example, the human vocal tract can only generate signals between roughly 50 and 10,000 Hz. This does not mean that every human can achieve that range, or that there are people who can generate frequencies slightly above or below. When it is said that human voice has a certain bandwidth, it does not mean that every instance of human voice will actually cover that entire range; it only means that you will not often find frequencies outside that range.
Other signals are limited in bandwidth by design. One example is voice in analog telephony; the voice is filtered to a range of approximately 300 to 3400 Hz, and it is impossible to find frequencies outside that range in the system. Another example is analog television: the signal is designed to have a bandwidth of 6 MHz, and except for equipment malfunction, you'll never find a TV signal with frequencies outside that range.


Most efficient way to find single dominant frequency (without amplitude) in analog signal I want to construct a reference locked sine oscillator. The multibit sine wave is calculated using DDS. The multibit reference signal is typically not sinusoidal, can have offset and substantial noise, but is expected to contain one dominant frequency.
If this dominant frequency is nearly known, a standard PLL approach will work to keep track of the dominant frequency in the reference signal.
However, I now suppose this cannot work to acquire initial reference lock if the dominant frequency is far away from the oscillator frequency. In this case, the PLL gets stuck on a locally dominant frequency.
What is the common approach to resolve this ? Calculating an DFT seems like overkill, since it produces a lot of unnecessary information. A recursive 4 sample Zoom FFT looks a bit more efficient. But maybe there is a way that directly modifies the PLL filter to capture further off frequencies ?
Please let me know if I missed crucial info.

One approach is to use a phase/frequency detector as the circuit determining error and to ensure that the initial loop bandwidth is wide enough to “see” both the weaker signal that would otherwise be a false acquisition as well as the strongest signal desired. The loop will lock onto the strongest signal within the loop bandwidth as long as it is 6 dB above the total power of the other tones (a phase vector diagram illustrates well the need for 6 dB margin as anything less causes abrupt 0/180 degree changes In phase such that a phase lock loop can’t converge to zero error). Once locked, the bandwidth can then be tightened for better SNR performance on the locked signal. This approach also leads to much faster acquisition time in general and is often referred to as “gear-shifting”.
If the signals are much wider than the acquisition bandwidth then more elaborate search and assessment methods are needed prior to setting the center frequency for acquisition. A separate frequency lock loop with wider bandwidth can be employed or other techniques that sweep over the entire bandwidth and access signal strength (such as the FFT or if processing time is available and resources scarce a linear frequency ramp can be employed). Note that the FFT approach can be significantly simplified since the number of bins (samples) needed in the FFT computation only need to be enough to ensure one bin for each step size equal to the loop acquisition bandwidth or less (so if the entire bandwidth of interest is ten times the loop bandwidth in its widest setting for acquisition or a real signal, then a 20 point FFT could be used to determine the location of the strongest bin).


Need to filter very small signal I am measuring ECG signals on subjects and I am having issues filtering a very small signal from the inherent noise of the collecting device. I am using an 18 bit ADC to digitize the signal
Here is an example of the noisy signal I get:

I would like it to look like this:

I have tried different filters and the best signal yet has been from the 8-point mean averaging of the signal:

However, I am trying to find the p-wave, which is a small upward deflection just before the large QRS complex of the pulse. This p-wave is buried in the noise and I can't seem to get it out. I have tried creating a median beat of 300 heartbeats, and by doing this I can show that the p-wave does exist, but I need to show that it exists for every heartbeat.
Median beat:

Is there a more elaborate way to tease out the p-wave using the median beat, or could you point me to other strategies for improving the signal to noise ratio?

The first thing you should do is examine your signal collection hardware and do everything humanly possible to improve the signal-to-noise ratio of the analog signal applied to your A/D converter. I realize this is easier said than done, but do your very best to reduce the noise in your analog signal.
Next, you might consider the following process:

Based on the following figure, examine your noisy signal samples and identify the $T_n$ time instants of your QRS complex spikes.
Create multiple signal segments, $S_n$, where each segment comprises $N$ samples just prior to and including the QRS complex spike sample.
Average the $S_n$ signal segments.
Examine the "averaged" $N$-sample segment to see if your can detect an averaged p-wave.


Understanding the Probability Distribution of Intensity Values From "Digital Image Processing 4th ed. - Gonzalez, Woods"

  Let $p_r(r)$ and $p_s(s)$ denote the PDFs of intensity values $r$ and $s$ in two different images. A fundamental result from probability theory is that if $p_r(r)$ and $T(r)$ are known, and $T(r)$ is continuous and differentiable over the range of values of interest, then the PDF of the transformed (mapped) variable $s$ can be obtained as:
  $$(1)\ \ p_s(s)\triangleq p_r(r)\left|\frac{dr}{ds}\right|$$

I should specify that $T(r)$ indicates a single-pixel intensity transformation. Now, I understand what the two PDFs indicate, but I fail to understand:

The derivation of (1)
The utility of (1)

EDIT: I found that $(1)$ is basically the change of variable formula for probability density functions, which holds when $T(r)$ is monotonic.

Let $R$ be the random variable "intensity" with probability density function $f_R(r)$ ($=p_r$). That implies that $R$ also has a cumulative density function $F_R(r) = P(R \le r) = \int_{-\infty}^r f_R(u)\,\mathrm du$.

Let $S = T(R)$ be another random variable, the transformed version of $R$, with the transformation mapping $s=T(R)$ smooth and monotonous¹. We're looking for the pdf of $S$, $f_S(s)$ ($=p_s$).

We take the convenient detour through $f_S(s) = \frac{\mathrm d}{\mathrm ds}F_S(s)$:

\begin{align}
\frac{\mathrm d}{\mathrm ds} F_S(s) &= \frac{\mathrm d}{\mathrm ds}P(S\le s)\\
&= \frac{\mathrm d}{\mathrm ds}P(S\le T(r))\\
&= \frac{\mathrm d}{\mathrm ds}P(T^{-1}(S)\le r)\\
&= \frac{\mathrm d}{\mathrm ds}F_R(T^{-1}(S)) & \text{chain rule of deriv.}\\
&= \frac{\mathrm d}{\mathrm d \left(T^{-1}(s)\right) }F_R(T^{-1}(S))\cdot \underbrace{\frac{\mathrm d}{\mathrm dS}T^{-1}(S)}_\text{deriv. of inv. function}\\
&=f_R(T^{-1}(S))\cdot \overbrace{\frac1{\frac{\mathrm d}{\mathrm dr}T(r)}}\\
&=f_R(R)\cdot \frac1{\frac{\mathrm d}{\mathrm dr}s}\\
&=f_R(R)\frac{\mathrm dr}{\mathrm ds} & \text{in original notation:}\\
&=p_r(r)\frac{\mathrm dr}{\mathrm ds}
\end{align}

I feel like I've forgotten the formal step that introduces the absolute value, but in case of $T(r)$, I required local invertability, so that kind of goes with that. If someone has a good notation for that, please do not hesitate to edit my answer. It's kinda late!

¹ that's an addition to the original problem, but for this way of derivation I need $T$ to be invertible. We can go a lot more measure-theoretic at this, but I don't really think that helps, so I'm going with this restriction, which I don't think is too harsh, as other mappings would lead to questionable results, just from a picture perspective.


Zeroing bins before IFFT Endolith's recent question on zeroing FFT bins got me thinking about a possible error in my workflow and I'd like to double check.

Consider a signal $x_i[m],\ m\in\{1,\ldots,M\}$, i.e. $M$ samples sampled at $f_s$ and the subscript denotes the $i^{th}$ such signal segment. I then filter this with a bandpass filter to between $0.1 f_s$ and $0.3f_s$ (at this stage, I was careful to not do a zero-the-bins kind of filtering) and Fourier transform it with an $N>2M$ long FFT. Let $X_i[k],\ k\in\{0,\ldots,N-1\}$ denote the Fourier coefficients and the filtered portion (including the roll-off frequencies on both ends) be the indices $0< k_l,...,k_u< N-1$.

My application now involves applying some function, $g$ to the Fourier coefficients $X_i[k_l],\ldots,X_i[k_u]$, averaging across the different segments ($i$) and inverse Fourier transforming the resulting vector. In math, I compute

$$X_{avg}[k]=\left\langle g(X_i[k])\right\rangle_i,\ k\in\{k_l,\ldots,k_u\}$$

where $\langle\cdot\rangle_i$ denotes averaging over the different segments. So far, so good.

Now, I didn't know what to do with the remaining Fourier coefficients $X_i[k'],\ k'\in\{0,\ldots,N-1\}\backslash \{k_l,\ldots,k_u\}$ before taking the IFFT. These coefficients are theoretically insignificant, since I've already filtered out content in them. So I went ahead and set $X_{avg}[k']=0$ and then computed the IFFT of $X_{avg}$. The results were along the lines of what I expected, but I didn't think of the aforementioned filtering issue, although I was aware of it (this is the same situation, with a bit of twists at places).

My question then is: Is it still a problem in this case i.e., when the bins that I "zeroed out" were technically low/insignificant? If so, how should I have handled it?

Filtering a signal by modifying the Fourier coefficients and then inverse Fourier-transforming it is a filter design method known as the "frequency-sampling" method.

Many research works on source separation for instance consider such filters because they are convenient (working on spectral characteristics to separate the signals - e.g. the fact that the frequency support of 2 sources is different) easy to implement and still lead to very convincing results (even the infamous "binary masks" :D).

You could compare this method to others, for instance FIR filter design by the window method : compute the ideal IIR impulse response and crop it or window it with a Hann/Hamming/your-preferred-window function. The window method corresponds to multiplying your impulse response by a window, in the time domain. It therefore results in a convolution between the involved Fourier transforms, leading to the well-known Gibbs effect.

For the frequency sampling method, you don't need to find the impulse response, but this comes at the cost that you do not control much of what is happening /between/ the frequency bins: in your example, that would correspond to the frequencies $\frac{\nu}{N} f_s$ where $\nu \notin \{0,1, 2,..., N-1\}$. The actual frequency response of such a filtering process is then an interpolation between the sampled points of the desired frequency response. This issue is not often discussed in the literature, but could lead to big cripples in the resulting frequency response.

There are other methods for designing filters, for instance for designing low/high pass or band-pass filters, but I have less practical experience there. I guess those are anyway not what you are looking for. I mentioned the window method only as a way of comparing what's happening: with the window method, you have trouble with the Gibbs effect, and with the frequency sampling method, you don't really know what's happening between the sampled frequency bins.

This said, for your practical problem, from my own experience, zeroing frequency bins is usually fine. Assuming the way you decide how to enhance this or that frequency bin makes sense (avoiding too many "jumps" in the frequency response, for instance), the result should not exhibit too many artifacts from that process. Note that setting the low values to zeros, instead of keeping the original ones, may however deteriorate what is happening in the frequency region of interest. In source separation as I have experience it, we usually try to minimize some criterion (e.g. mean squared estimation errors), and that can provide us with an optimal way of setting all the values of the Fourier transform. I would argue that this is a better way of dealing with the matter, instead of trying all sort of values for the "missing" ones until the result is acceptable (a method which is most likely to depend on the test signal, hence poorly generalizing to other signals).


Finding Fourier series coefficients for discrete time signal Let $x[n]$ be a periodic sequence with period $N$ and Fourier series representation $$x[n] = \sum _{k=<N>}a_ke^{jk\frac{2\pi}{N}n}$$ Determine the Fourier series coefficients for
$$y[n]  =
\begin{cases}
x[n],  & \text{if $n$ is even} \\[2ex]
0, & \text{if $n$ is odd}
\end{cases}$$

My try
$N$ even
We have $x[0], 0 , x[2] , 0 , \dots x[N-2] , 0 , x[N]$.
Since $x[0] = x[N]$ the sequence $y[n]$ is periodic with $N$. Also we have
$$e^{jM\frac{2\pi}{N}n}x[n] \leftrightarrow a_{k-M}$$Let $M = \frac{N}{2}$, then $(-1)^nx[n] \leftrightarrow a_{k-\frac{N}{2}}$.
In this case $$b_k = \frac{1}{2}(a_k + a_{k-\frac{N}{2}})$$
$N$ odd
By writing the first terms of $y[n]$, we see that $y[n]$ is periodic with period $2N$
$$x[0], 0 , x[2] , 0 , \dots x[N-1] , 0 , x[N+1] , 0 , x[N+3] , \dots , 0 , x[2N]$$Which is $$x[0], 0 , x[2] , 0 , \dots x[N-1] , 0 , x[1] , 0 , x[3] , \dots , 0 , x[0]$$
So we write
\begin{align}
b_k &= \frac{1}{2N}\sum _{n=<2N>}y[n]e^{\frac{-jk\pi n}{N}}\\
 &= \frac{1}{2N}     \sum _{n=0}^{\frac{N-1}{2}}x[2n]
                e^{\frac{-jk2\pi n}{N}} +
    \frac{(-1)^k}{2N}\sum _{n=1}^{\frac{N-1}{2}}x[2n-1]
                e^{\frac{-jk\pi (2n-1)}{N}}.
\end{align}
I don't know how to proceed further and how to relates $b_k$ and $a_k$.

You are asking for DFS (discrete Fourier series) coefficients $b_k$ of the periodic sequence $y[n]$, in terms $a_k$ of the periodic sequence $x[n]$. Since DFS and DFT (discrete Fourier transform) are the same things, I will instead write down an answer for the DFTs $Y[k]$ and $X[k]$ of the associated sequences, in the context of which $X[k] = a_k$ , and $Y[k] = b_k$ will be understood.
I would like to relate DFTs to DTFTs (discrete-time Fourier transform) through the sampling relation, which indicates that :

N-point DFT X[k] of the N-point sequence x[n], is the uniform samples of the DTFT X(w) of
x[n].

The DTFT of $x[n]$ is :
$$ X(\omega) = \sum_{n=0}^{N-1} x[n] e^{-j \omega n } \tag{a}$$
and the DFT of $x[n]$ is :
$$ X[k] = X(\frac{2\pi}{N}k) = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N}nk } \tag{b}$$
Let $x[n]$ and $y[n]$ be two sequences related by:
$$y[n]=\begin{cases}{ x[n] ~~~,~~~n:even \\  ~~~0 ~~~~~,~~~ n:odd }\end{cases}$$
You can also express $y[n]$ as :
$$ x[n] \longrightarrow \boxed{ \downarrow 2 } \overset{w[n]}{\longrightarrow} \boxed{ \uparrow 2 }    \longrightarrow y[n] $$
Applying DTFT relations to the sequences :
$$ W(\omega) = 0.5 \left( X\left( \frac{\omega}{2} \right) +  X\left( \frac{\omega-2\pi}{2} \right) \right) \tag{1} $$
and
$$ Y(\omega) = W(2\omega) \tag{2} $$
hence
$$ Y(\omega) = 0.5 \left( X\left( \omega \right) +  X\left( \omega-\pi\right) \right) \tag{3} $$
Eq.3 describes the relation (for any $N$) between the DTFTs $Y(\omega)$ and $X(\omega)$ of the sequences $x[n]$, and $y[n]$. To turn this into the equivalent relation between the N-point DFTs $X[k]$ and $Y[k]$, (and for your case a relation between DFS  $a_k$ and $b_k$) we shall sample the associated DTFTs. Then we have :
$$  \begin{align}
Y[k] &= Y( \frac{2\pi}{N} k ) ~~~,~~~k = 0,1,...,N-1 \tag{4} \\ \\
     &= 0.5 \left( X\left( \frac{2\pi}{N} k \right) +  X\left( \frac{2\pi}{N} k-\pi\right) \right) \tag{5} \\ \\
     &\boxed{ Y[k] = 0.5 \left( X\left[ k \right] +  X\left[k- \frac{N}{2} \right] \right) }\tag{6} \\
\end{align}$$
Eq.6 describes the requested relation between N-point DFTs $X[k]$ and $Y[k]$.
When $N$ is EVEN (say $N = 2 M$ ) the relation becomes :
$$ \boxed{ Y[k] = 0.5 ( X[k] + X[k-M] ) }\tag{7} $$
which is what you have arrived at $b_k = 0.5(a_k + a_{k -\frac{N}{2}})$.
However, when $N$ is ODD, there is a problem due to $N/2$ not being an integer. The sequence $X[k-N/2]$ is interpreted as an interpolation of $X[k]$ at the fractional indices $k-N/2$.
For $N=2M+1$ , we have $X[k-N/2] = X[k-M-0.5]$ , $k=0,1,...,N-1$ , which evaluates $X[k]$ at the points $m = 0.5,1.5,...,N-0.5$. To get those intermediate samples, we need an interpolation of $X[k]$ by 2.
Interpolation of $X[k]$ by 2, is achieved by zero padding $x[n]$ by $N$ samples, and computing its $2N$-point DFT :
$$ x_e[n] = \begin{cases}{ x[n] ~~~,~~~ 0 \leq n < N \\ ~~~0 ~~~~ , ~~~~ N \leq n < 2N }\end{cases} $$
Then we can see, by the DFT - DTFT sampling relation that :
$$ \mathcal{DFT}\{x_e[n]\} =X_2[k] = X(\frac{2\pi}{2N}k )= X(\frac{2\pi}{N} k/2 ) = X[k/2] \tag{8} $$
for $ k = 0,1,.,2N-1$. Note that DTFTs of $x[n]$ and $x_e[n]$ are the same.
Based on Eqs. 8, you can restate Eq.6 when $N$ is odd :
$$ \boxed{ Y[k]= 0.5 \big( X\left[ k \right] +  X_2\left[2k- N\right] \big) } \tag{9} $$
with your notation it becomes:
$$ b_k = 0.5 ( a_k + a^e_{2k-N} ) $$ where $a^e$ is the DFS of the sequence $x_e[n]$.
As you can see, you cannot express $b_k$ (with a simple formula) in terms of sole $a_k$ when $N$ is odd.


Phase correction using second order all pass filters/LMS with Wiener-Hopf equations I am working on phase correction algorithms. Can anyone help me understand how can the following diagram be used for phase correction:

Additionally, for a channel function suggested How determine the delay in my signal practically, I have tried the following:
Fs = 306.700*1e+3; % sampling rate
Fc = 100*1e+3; % central frequency
B = 4*1e+3;  % bandwidth
N = 3072;
T = (N-1)/Fs;
t = 0:1/Fs:T; %time

baseband_signal = chirp(t,-B/2,t(end),B/2)-i*chirp(t,-B/2,t(end),B/2,'linear',90);

% Modulation
modulated_signal = baseband_signal.*exp(i*2*pi*Fc.*t);
modulated_signal = real(modulated_signal);

% Transmitted signal
tx = modulated_signal.*1.25/max(modulated_signal);

% Received signal
rx = [zeros(1,200) tx];
t_received = 0:1/Fs:(length(rx)-1)/Fs;

% Computing equalization filter: rx and tx swapped!
 coeff = channel(rx,tx,20);

% Equalization filtering
 rx_corrected = filter(coeff, 1,  rx);

 figure
freqz(coeff)
title('Equalization filter')

figure
plot(t_received,rx,t_received,rx_corrected)
xlabel('Time (s)');
ylabel('Amplitude');
grid on;
legend('initial','corrected')

Shouldn't the equalization suggested here remove the phase delay and have the equalized and transmitted signal aligned?

I actually don't see how the following diagram can be possibly used for phase correction only as it will have an unavoidable amplitude correction as well. An all pass filter by definition will pass all magnitude without modification and only change the phase. The Time Delay block will similarly pass all magnitude without modification with a linear phase versus frequency. The resulting amplitude at the output will be the vector sum of the linear phase component with the all pass component. Since both will have the same magnitude then the sum will be dictated by the phase (for example if there are any frequencies where the phase for each are related by $pi$ radians there will be complete cancellation. The only way the resulting diagram can maintain a constant magnitude is if the all pass filter similarly has a linear phase (which is a trivial all-pass and further would not modify phase to perform phase correction).
To do phase correction you would use the all pass filter directly (this is its purpose) which as explained above is a filter that modifies only the phase response.


How does decimation work without losing samples? For example, changing the frequency from 32.768Mhz to 2kHz there is a scale down of 16384.
From the algorithms I've seen you have to discard 16384-1 samples and only count the 16384th sample. However, that many samples being lost will lose important parts of the information.
How does decimation work then?
EDIT
A signal is sampled at 2kHz and then modulated with a carrier wave of signal frequency 8.192MHz which is sampled at 32.768MHz. I am trying to revert it back to 2kHz, by decimation. I assume pass it through a LPF and  apply decimation however, the points sampled at 32.768MHz are then lost?

Your understanding is correct. You need to make sure the information content fits into the output rate, and that implies your signal needs to be band-limited to 1/16384. Otherwise you'll lose information.
Making sure that's the case with a filter is called anti-alias filtering.


Phase plot in Bode Diagrams I'm studying the transfer function of an aircraft that relates the pitch angle, $\theta$, with the horizontal stabilizer angle, $\delta$.
$$ H(s) = \dfrac{\theta (s)}{\delta (s)}$$

I have computed the frequency response by using Matlab, and I plotted the Bode Diagram: amplitude plot and phase plot.

I understand what is the amplitude curve indicating and I'm able to comment the result obtained. But when it comes to studying the phase curve (second figure above) I don't know what is the meaning and how it is related with the other curve.

So, my question is: what can I say from the phase plot by looking at the Bode Diagram of the studied transfer function, $H(s)$?

The Bode diagram shows you the properties of the transfer function $H(s)$ for $s=j\omega$, i.e. the complex frequency response. Since $H(j\omega)$ is generally a complex function it has a magnitude and a phase:

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

The plot shows the principal value of the phase in the interval $[-\pi,\pi)$, that's why you see a jump (of $2\pi$), which is no actual phase jump but it just comes from mapping the phase back to the interval $[-\pi,\pi)$.

The magnitude of $H(j\omega)$ tells you the relative amplification/attenuation of the different frequency components, whereas the phase shows you the relative alignment of the different components.

For example, a sinusoid $x(t)=A\sin(\omega_0t+\theta)$ will appear at the output of the system as

$$y(t)=A|H(j\omega_0)|\sin(\omega_0t+\theta+\phi(\omega_0))$$

where $\phi(\omega_0)$ is the system's phase response evaluated at the input frequency $\omega_0$.


Effect of sampling rate on BER simulation I am trying to build a bit error rate simulation for a digital baseband signal
(specifically using Manchester coding). The simulator generates a digital waveform from random symbols at a sample rate of at least $f_s = 2/T_l$, where $T_l$ is the chip duration (the duration of the "on" time).

The transmitted signal is mixed with AWGN, based on the desired S/N ratio. The detector correlates each symbol period from the received signal with an ideal 0 waveform or 1 waveform, and chooses the best. Either "ideal" is 1V during the on time, and 0V during the off time. These ideal waveforms are sampled at the same rate as the input.

What I am seeing is that the BER improves when the sampling rate is increased, for a given SNR (see plot). This was unexpected, since I believed that the error rate would be independent of the sampling rate.

Am I missing something in this BER simulation? Shouldn't the error rate be independent of the sampling rate?

Another hypothesis: as you increase the sample rate you also increase the bandwidth. Depending on how you generate the noise, this probably means that the spectral density of your noise decreases. If the spectral content of your signal stays constant, that means that the in-band SNR goes up.


What happens at the edges in fft2 filtering? I am doing a fairly simple homework assignment in MATLAB that is to filter an image in the frequency domain. So my code looks something like this:

output = ifft2(fft2(image, x, y).*fft2(filter, x, y));

After testing some images, this code filters properly, but has weird edge behavior. Out of curiosity, what exactly happens at the edges using this method? (0, symmetric, replicate, circular)?

Quick answer: FFT2 assumes periodic boundary conditions. This means that performing the FFT2 of an image is equivalent to have your original image periodically repeated in the horizontal and vertical directions, compute the DFT, and then only retain the central region. The weird artifacts you are seeing are probably the result of juxtaposing opposite edges that are not related to each other.

There are ways to work around this by, e.g., smoothing these discontinuities by using MATLAB’s edgetaper function before filtering.


Conceptual Question on Least Squares Estimation Method What is meant  by Least square plus Gaussian method of estimating the parameters of an unknown system? I am aware of least squares  in general but some works refer it as least square+gaussian method. What is the role of Gaussian here? Is it for some kind of assumption of the distribution of the residual errors? please clarify

Gauss did a two-way proof.  In rough terms, he not only showed that least squares gives the best fit to data containing a standard normal distribution of errors, but that if least squares is always giving the best fit, that implies that the distribution is indeed standard normal (also now called Gaussian).

The correlation coefficient, often computed in conjunction with a least square fit, is usually assumed to be parameterizing an underlying Gaussian distribution.  If the distribution is otherwise, it may well infer nonsense.


Is there any way to align temporally two signals in spectrogram? I am keen to know that If I have two signals say :
Farend_signal (y[n]) and Nearend_signal (x[n])
and x[n] also contains an echo of y[n] with some delay.
Is there any technique to align these two signals temporally is Spectrogram?
For example, any correlation technique to find the temporal displacement between signals after making a spectrogram(STFT) of two signals?

You can find a large number of  "model-based echo suppression" or "adaptive acoustic echo cancellation" (for the keywords). They often combine a sparsifying transformation (wavelets, time-frequency decompositions) and some form of fitting or adaptation of the model to the time and frequency localization where the echo dwells. Underlying hypotheszs are:

the model somehow looks like the echo,
the similarity will increase in the transformed domain,
this representation will help the morphing, or make is more efficient (faster, etc.)

For voice, you find a look of literature, overviews, etc. You may also look at neighboring domains: vibration-based non-destructive testing or monitoring, ultrasounds. I would mention seimic processing, where echoes are sometimes named multiples (waves bounding several times between layers). Here are two examples I have been involved in. A first one quite fast which uses 1-tap adaptive filters (hint: a complex filter), the second  more involved with optimization constraints. They use time-scale wavelets, but could be readily adapted to short-term Fourier transforms or multi-input/multi-output filter banks.

Adaptive multiple subtraction with wavelet-based complex unary Wiener filters, 2012  (preprint)
A Primal-Dual Proximal Algorithm for Sparse Template-Based Adaptive Filtering: Application to Seismic Multiple Removal, 2014 (preprint)


Manual filter or digital filter first for EEG? I am getting conflicting advise regarding how to clean my EEG data:

1) Manually remove artefacts first and then apply digital filters

OR

2) Apply digital filters first and then manually remove artefacts

The reason given for 1) is because artefacts are more visible and avoids accidentally accepting artefacts as EEG data.

The reason given for 2) is because it avoids accidentally rejecting EEG data.

Both reasoning makes sense to me and I am a little confused as to what I should do: 1) or 2). Anyone can advise? Thank you very much.

When it comes to preprocessing of EEG data, the first step is to filter the signal. Filtering is done such that it preserves information across different bands of the EEG data (alpha, beta, gamma, delta). Typically, many researchers use 0.5 - 50 Hz band in the first step (removes DC components). Note that the filter should be zero-phase such that no delay or phase changes are introduced in the EEG data on using the filter. One simple way is to use filtfilt command in MATLAB.

In addition, a notch filter is also designed to remove 60 Hz. You might be wondering why we need a notch filter at 60 Hz if we are filtering the original signal using a band-pass filter of 0.5 - 50 Hz. To answer this, we need to first understand the structure of the band-pass filter. If the band-pass filter is of lower order, then the transition bands are not very steep (or sharp). This in turn introduces some part of the 60 Hz interference into the band-pass filtered signal (which typically dominates in EEG data). So, to suppress the strong 60 Hz interference, a notch filter is also used.

After filtering the EEG data using a band-pass and a notch filter, you might want to address different types of artifacts that show up in the EEG data. Some of the them include EKG artifact, movement artifact, sweat artifact, etc. Each artifact removal requires a sophisticated approach such that the EEG data is preserved but the artifact is removed completely.


Why can't music recognition services like Shazam or SoundHound recognize music if its speed is changed by more than 5 percent? Music services like Shazam and SoundHound take an audio sample from your microphone and compare it to its database to recognize music. This algorithm works very well under high noise, unequal frequency bins, and even a simple hum or whistle to the music. However, if you change the song's pitch or speed by more than 5%, the algorithm cannot detect the song at all. A 5% change is very small and not even that noticeable, yet it completely fools the music detection algorithm. Why is this? Is it due to the different speed/pitch missing the targets in the audio fingerprint?

this is a guess from someone who does audio, but has never ever done recognition nor watermarking algorithms for someone in the industry (but i have a little idea how they may be done).

my guess is that they extract some features from the STFT of the music.  associated with each feature is an envelope that represents the quantity of that feature at the center window time of the STFT.

the collection of those envelopes can be thought of as a point in a metric space.  if you have a metric space, you can define a difference function

$$ \| \mathbf{X} - \mathbf{X}_n \| $$

between those of any pair.

then they compare the test envelope, $\mathbf{X}$, with each template, $\mathbf{X}_n$, of a set that come from the original targets (the copyrighted songs they're fishing for).  there is a norm taken of the difference between these points in a metric space, between the collection of envelopes taken from the test sound and the template.  if the norm is below a given threshold and if the norm is the minimum compared to all of the templates, then the song of that best matching template is identified as matching the input or test sound.

now if the pitch is wrong it won't match as good.  but they could make their difference algorithm insensitive to a constant offset to pitch.  it's like running the difference signal through a DC blocking filter before calculating the norm of it.  but they would have to do that, otherwise a constant pitch difference between test tune and target will always be an error.

and they could match a tempo difference, but that is much harder.  if memory and computational effort is not a problem (and they often are not for non-realtime applications running on a GHz-speed computer or server) then there could be a template for each of many speeds.  and see which one matches best.  but if these other tempo templates do not exist, and you submit a tune that does not exactly match a template of it at the original speed, then at some point in the tune or song, the pitches will not match at all.  that will make a dumb algorithm think they're different tunes.


I B P video frame reordering I am studying H264 video these days. I know video encoder will do frame reorder for GOP with B pictures.I am not sure if my understand is correct!

GOP input to encoder (it's raw data but I just map it as GOP for better understand).

I(1) B(2) B(3) P(4) B(5) B(6) P(7) B(8) B(9) P(10) B(11) B(12) I(13)

GOP output from encoder.

I(1) P(4) B(2) B(3) P(7) B(5) B(6) P(10) B(8) B(9) I(13) B(11) B(12)

GOP input to decoder.(the same as output from encoder)

I(1) P(4) B(2) B(3) P(7) B(5) B(6) P(10) B(8) B(9) I(13) B(11) B(12)

GOP output from decode. (the same as input to encoder)

I(1) B(2) B(3) P(4) B(5) B(6) P(7) B(8) B(9) P(10) B(11) B(12) I(13)

I think you have in general the correct idea. But you have to remember that the encode order (and therefore the order in which the pictures will appear in the bitstream and be decoded), depends on how the references are.

if you have: I(1) B(2) B(3) P(4) B(5) B(6) P(7). etc...

You could have a structure where in the group [B(2) B(3) P(4)], the first B(2) has a reference to the past I frame, but a future reference to the following B(3) frame. In that case, the encode order start the same as having I(1) P(4), but the next cannot be B(2), because B(2) needs to wait for B(3) to be encoded. And this could change from repetition to repetition.

So in general the encode order of frames depends on how the referencing structure is. Since the standards do not specify this and allow freedom of choice when encoding, the GOP structure is in general totally arbitrary, so every specific case has to be analyzed.


Image Interpolation Using the Yule Walker Equations I have been studying about the Yule-Walker equations for prediction of a time series data from knowledge of past values of the series. Is there any way I can use the same in an image to exploit the correlation between image pixels in order to perform image interpolation,, or to increase the resolution of the image?

The main problem is Auto Regressive (AR) / IIR models are based on the idea of something before and after.
While for time series it is well defined, how will you address that in image for a pixel in the middle of the image?

I addressed some ideas on that in my answer to the question Denoising an Image Using Kalman Filter.

Another idea is to work across rows and then across columns.
Then you can build an AR model based on the samples with unit "Time Parameter" and then ask what should happen for $ n / 2 $ or so.
You can even take it farther, and in order to use the 2D information take $ k $ rows / columns at a time and build a vectorized AR model.
Then it means that for each pixel you'll get $ k $ estimations which you could average or chose any other idea to fuse.

Actually it could be nice idea to have a try...


Linear-phase lowpass filters: frequency responses alignment These two filters have very similar frequency responses, but those responses are not identical except at a single frequency $\omega_0$. I want to find $\omega_0$.

One way I've seen that done is like this:

  $$\begin{align}
Y_1(z)&=X(z)+2X(z)z^{-1}+X(z)z^{-2}\\
Y_2(z)&=X(z)+X(z)z^{-1}
\end{align}$$

  if $Y_1(z)=Y_2(z)$

  $$
\begin{align}
\Longrightarrow 1+2z^{-1}+z^{-2}&=1+z^{-1}\\
\Longrightarrow z^{-1}+z^{-2}&=0\\
\Longrightarrow z\left(z+1\right)&=0\\
\Longrightarrow z=e^{j\omega}=-1\\
\Longrightarrow \omega=\pi\\
\omega_0=\frac{f_s}2
\end{align}
$$

But, to be perfectly frank that is quite confusing to me I don't quite know what's going on there, would someone please break it down for me?

Another way I've seen it done is like this:

  $$
H_1\left(\omega_0\right)= H_2\left(\omega_0\right)\rightarrow 1+2e^{-j\omega_0}+e^{-j2\omega_0}=1+e^{-j\omega_0}.
$$

  With algebraic simplification the above equality can be written as:

  $$
e^{-j\omega_0}=-e^{-j2\omega_0}.
$$
  Here's the tricky part,...we multiply both sides of the above equation by $e^{j\omega_0}$ to give us:
  $$
1=-e^{-j\omega_0}.
$$

I guess they're equivalent, but, in what ways specifically are they identical and in what ways are they divergent?

If someone could help me to see I would be so grateful.

Both methods are equivalent. Look at it like this: first write down the frequency responses (as already indicated in your question)

$$H_1(e^{j\omega})=1+2e^{-j\omega}+e^{-2j\omega}\\
H_2(e^{j\omega})=1+e^{-j\omega}$$

Now want to find $\omega_0$ such that

$$H_1(e^{j\omega_0})=H_2(e^{j\omega_0})$$

which gives

$$1+2e^{-j\omega_0}+e^{-2j\omega_0}=1+e^{-j\omega_0}$$

or, equivalently,

$$e^{-j\omega_0}(1+e^{-j\omega_0})=0\tag{1}$$

Equation (1) can only be satisfied if the term in parentheses vanishes. So you get

$$e^{-j\omega_0}=-1\Longrightarrow\omega_0=\pi$$

And since $\omega_0=2\pi f_0/f_s$ (where $f_s$ is the sampling frequency) you finally get

$$f_0=f_s/2$$


Why the unilateral Laplace transform? Why is the Laplace transform commonly taught as the unilateral Laplace transform?

I mean, for the Fourier transform, we commonly have the bilateral transform... if the signal is 0 for $t<0$, then it turns into a unilateral Fourier transform. Why not have this same convention for Laplace transform? Why specifically introduce the unilateral version?

The widespread use of the unilateral Laplace transform reflects the fact that in practice we often deal with causal systems and signals that have a defined starting time (usually chosen as $t_0=0$).

The Fourier transform is mainly used for analyzing ideal signals and systems, such as ideal filters (e.g., low pass, high pass, etc.) and ideal signals such as perfect sinusoids. In these cases we have to deal with non-causal systems with impulse responses that extend from $-\infty$ to $\infty$. The same is of course true for sinusoidal signals or complex exponentials. Note that neither ideal filters nor the signals mentioned above can be treated by the Laplace transform.

One of the most important features of the unilateral Laplace transform is that it can be used to elegantly solve differential equations with initial conditions. The initial conditions are taken into account by the well-known differentiation property of the unilateral Laplace transform:

$$\mathcal{L}\{f'(t)\}=sF(s)-f(0^-)\tag{1}$$

where $f(t)$ is a differentiable function, and $F(s)$ is its (unilateral) Laplace transform. The Fourier transform doesn't have an equivalent to $(1)$ for taking initial conditions into account. If the Fourier transform is to be used for solving a differential equation with non-zero initial conditions, then the initial conditions need to be modeled as additional sources.


Define comb function in terms of exponential function I have one signal with $N$ samples and I want to downsample it by factor of $M$. Essentially my downsampler is comb function which takes every M$^{th}$ sample while others are set to zero, i.e.,

$$
\check{Y}_k=\sum_{n=0}^{N-1}\sum_{m=0}^{\frac{N}{M}-1}x[n]\delta(n-mM)e^{-j\frac{2\pi nk}{N}}\tag{1}
$$
where $x[n]$ is original signal, $\check{Y}_k$ is downsampled signal, and $n$ are time samples .

How is comb function $\delta(n-mM)$ define in terms of exponential function ? i.e.

$$\sum_{n=0}^{N-1}\sum_{m=0}^{\frac{N}{M}-1}\delta(n-mM)=\sum_{n=0}^{N-1}\sum_{m=0}^{\frac{N}{M}-1}e^{??}$$

Also what is the definition of comb function in frequency domain?

The frequency domain representation of the discrete time downsampling operation (i.e.  $X_M(e^{j\omega})$ , the $DTFT\{x[Mn]\}$, in terms of $X(e^{j\omega})$, the $DTFT\{x[n]\}$ ) can be most easily observed by the following approach. We will begin from the end, $x[Mn]$ and proceed to its pre-compressed form by first deriving the frequency domain relation between the compressed signal $x[Mn]$ and its expanded form (pre-compressed) $x_s[n]$ and then we shall show the spectral relation between an input signal and its sampled form $x_s[n]$ (which is also the expanded signal at the input of the compressor)

First consider a periodic pulse train ,  $p_M[n] \triangleq \sum_{k=-\infty}^{\infty}{ \delta[n-kM]}$ whose period is $M$ and defined as :

$$p_M[n]=
\begin{cases}
1 ~~~~, ~~~ \text{n = 0, M, 2M, ..., kM , ...} \\
0 ~~~~, ~~~ \text{otherwise} \\
\end{cases} \tag{1}$$

Then the signal $x[n]$  will be first multiplied by this perodic pulse train $p_M[n]$, to produce $x_s[n]$ which will be selecting the required samples as follows:

$$x_s[n]=x[n]\cdot p_M[n] =
\begin{cases}
x[n] ~~~~, ~~~ \text{n = 0, M, 2M, ..., kM , ...} \\
0 ~~~~, ~~~ \text{otherwise} \tag{2}\\
\end{cases}$$

Now, we shall take every $M^{th}$ sample of $x_s[n]$ and discard the rest, which is called as the compression operation as graphically represented by:

$$ x_s[n] \longrightarrow \boxed{\downarrow M} \longrightarrow x_s[Mn]=x[Mn]=x_M[n] \tag{3}$$

Note that eventhough the conversion between $x[n]$ and $x[Mn]$ is not invertable, the conversion between $x_s[n]$ and $x_s[Mn]$ is so, as should be clear:

$$ x_s[Mn]=x_{sM}[n] \longrightarrow \boxed{\uparrow M} \longrightarrow x_s[n] \tag{4}$$

i.e. the inverse is the expander operator and we shall exploit it for easily showing the spectral relation between $x_s[n]$ and $x_M[n]$ as it can be shown that:

$$X_{s}(e^{j\omega}) = X_{sM}(e^{j\omega M}) \Longrightarrow  X_{sM}(e^{j\omega}) = X_{s}(e^{j\omega/ M}) \tag{5}$$

Based on this we can argue that (noting the equation (3)), $$ \boxed{ DTFT\{x[Mn]\} \triangleq X_M(e^{j\omega}) = X_{s}(e^{j\omega/ M}) } \tag{6}$$ where $X_{s}(e^{j\omega})$ is the DTFT of the sampled signal $x_s[n]$.

What comes next is the relation between the DTFTs of the signals $x_s[n]$ and $x[n]$ to complete the picture, which will be a straight forward observation if we consider representing the periodic impulse train by another form namelmy by its DFS (Discrete Fourier Series) representations (as is possible for any periodic discrete time signal with period M) :

$$ p_M[n] \triangleq \sum_{k=-\infty}^{\infty}{ \delta[n-kM]} = \frac{1}{M}  \sum_{k=0}^{M-1}{e^{j \frac{2\pi}{M}kn} } \tag{7}$$

Then we shall express $x_s[n]$ as :
$$ x_s[n] \triangleq x[n]\cdot p_M[n] = x[n] \big( \frac{1}{M}  \sum_{k=0}^{M-1}{e^{j \frac{2\pi}{M}kn} }\big) = \frac{1}{M}  \sum_{k=0}^{M-1}{x[n]e^{j \frac{2\pi}{M}kn} } \tag{8}$$

What comes next is using the linearity and modulation properties of DTFT operator and argue that:
$$DTFT\{x_s[n]\} = DTFT\{\frac{1}{M}  \sum_{k=0}^{M-1}{x[n]e^{j \frac{2\pi}{M}kn} } \} = \frac{1}{M}  \sum_{k=0}^{M-1}{ DTFT \{ x[n]e^{j \frac{2\pi}{M}kn} \} } = $$ $$ \boxed{X_s(e^{j\omega})=\frac{1}{M}  \sum_{k=0}^{M-1}{ X(e^{j \big(\omega -\frac{2\pi}{M}k \big)} ) } } \tag{9}$$

The final step is the mergal of the steps (6) and (9) to produce the result:
$$ \boxed{ X_M(e^{j\omega}) = X_s(e^{j\omega/M})=\frac{1}{M}  \sum_{k=0}^{M-1}{ X(e^{j \big(\frac{\omega}{M} -\frac{2\pi}{M}k \big)} ) } = \frac{1}{M}  \sum_{k=0}^{M-1}{ X(e^{j \big(\frac{\omega-2\pi k}{M} \big)} ) } } \tag{10} $$

I hope it was a clear way of expressing the spectrum of downsampled signal $x_M[n]=x[Mn]$ to that of the original signal $x[n]$

Also note that a more complete understanding of the spectral relation would only be possible if you fully explore the relation by a graphical plot as well, which would clearly demonstrate the resulting frequency stretching towards the $\omega  =\pi$ boundary, by a factor of $M$ after which a consideration of aliasing becomes clear...


Upsampling DFT at frequency domain Suppose we have a vector of length N of DFT, which every odd index (in frequency) equals 0:

$$X[k] = 0\quad \forall k:\, k\mod2\equiv1$$

What does that mean for the series in the time domain?

It can mean the following:
Consider any signal $x_{N}[n]$ of length $N$, and form the length $2N$ signal $x_{2N}[n]$ by repeating a copy of $x_{N}[n]$, then compute $2N$-Point DFT of this new signal $x_{2N}[n]$ whose every odd bin will be zero. i.e. $$X_{2N}[k] = 0 ~~,~~ \text{for}~~ k=1,3,5...,2N-1$$

whose demonstration can be performed with the following Octave/Matlab code:

 N = 16;
 x = randn(1,N);         % form an arbitrary signal x[n]
 x2 = [x x];             % repeat it once
 X2 = fft(x2,2*N);       % take DFT of length 2N
 figure,stem(abs(X2));   % display the result

Therefore it can mean that if you are performing an $M$ point DFT on your original signal $x[n]$, the signal $x[n]$ could actually be periodic with $\frac M2$ for $M$ even.

A more mathematical derivation of the result can be seen as follows:

Consider a signal $x_N[n]$ of length $N$ and its repetition $x_{2N}[n]$ of length $2N$ such that $x_{2N}[n]=x_N[n \mod N]$ for $n=0,1,\ldots,2N-1$.

Let's denote $N$-Point DFT of $x_N[n]$ as $X_N[k]$ for $k=0,1,\ldots,N-1$ and $2N$-Point DFT of $x_{2N}[n]$ as $X_{2N}[k]$ for $k=0,1,\ldots,2N-1$:

\begin{align}
&{\text{Start with defining $2N$-Point DFT of $x_{2N}[n]\ldots$}}\\
X_{2N}[k] &= \sum_{n=0}^{2N-1} {x_{2N}[n] e^{-j\frac{2\pi}{2N}nk}} \\ \\
&{\text{first split the sum and then substitude $n=n-N$ in the $2^{nd}$}}\\
&= \sum_{n=0}^{N-1} {x_{2N}[n] e^{-j\frac{2\pi}{N}n(k/2)}} + \sum_{n=N}^{2N-1} {x_{2N}[n] e^{-j\frac{2\pi}{2N}nk}} \\ \\
&{\text{recognise $x_{2N}[n]=x_N[n]$, $x_{2N}[n+N]=x_N[n]$ }}\\
&= \sum_{n=0}^{N-1} {x_{2N}[n] e^{-j\frac{2\pi}{N}n(k/2)}} + \sum_{n=0}^{N-1} {x_{2N}[n+N] e^{-j\frac{2\pi}{2N}(n+N)k}} \\ \\
&{\text{Expand the 2nd sum's multiplier factor}}\\
&= \sum_{n=0}^{N-1} {x_N[n] e^{-j\frac{2\pi}{N}n(k/2)}} + \sum_{n=0}^{N-1} {x_N[n] e^{-j\frac{2\pi}{N}(n+N)(k/2)}} \\ \\
&{\text{Replace $e^{-2j\frac\pi N}$ with $W_N$ for simplicity }} \\
&= \sum_{n=0}^{N-1} {x_N[n] W_N^{nk/2}} + W_N^{Nk/2}\sum_{n=0}^{N-1} {x_N[n] W_N^{nk/2}} \\ \\
&{\text{Recognise the sums as $N$-Point DFT $X[k]$ of $x[n]$, at $\frac k2$}}\\
&= X_N[k/2] + e^{-j\pi k} X_N[k/2] \\ \\
&= X_N[k/2] \cdot \big( 1 + (-1)^k \big)  \\ \\
&= \begin{cases}
2 X_N[k/2] , &\scriptstyle{\text{k=0,2,4,...,2N-2 }}\\
0  , &\scriptstyle{\text{k=1,3,5,...,2N-1 }}\\
\end{cases}
\end{align}

Shows that for such a defined signal , its odd DFT bins will be zero...