Removing a sinusoidal artifact from a set of movie frames I am doing some post-hoc analysis of a dataset consisting of a series of movie frames that are contaminated by a strongly periodic artifact. I would like to remove this artifact from my frames.

For ease of plotting I've just reshaped my array M of pixel values to [nframes, npixels], then averaged over all pixel values to give me a 1D vector m. Here's what this signal looks like in the time domain. You can see the oscillation quite clearly in the zoomed inset.

I then made a periodogram by taking Fm = rfft(m), and plotted abs(Fm)**2 against frequency. I see a very sharp peak at ~1.5Hz:

As well as the temporal periodicity, there also seems to be a weaker spatial component to this artifact, since at the exact peak frequency value there seems to be a smooth variation in phase across the x-axis of my frames, so that pixels on the right tend to lag pixels on the left:

As a brute force approach, I've tried just filtering each pixel in the time domain using a notch filter centred on 1.5Hz. I used an order 4 Butterworth filter with critical frequencies 1.46 and 1.52Hz (I'm not well versed in filter design, so I'm sure there may be more appropriate choices).

Here's what the mean pixel signal looks like after filtering:

And the corresponding periodogram:

The notch filter does a reasonably good job of reducing the artifact, but since it basically looks like a pure stationary sinusoid I can't help but think that I could do better than just attenuating that part of frequency space.

My initial (very naive) idea was to do something like:

Get the frequency, phase and amplitude of the oscillation from the Fourier spectrum for each pixel in the movie
Reconstruct the oscillation in the time domain
Subtract it from the movie frames

I realise this isn't something people usually do, since interference usually isn't so spectrally pure and temporally stationary, but I wonder if it might make sense in my case?

Data

Full 16bit TIFF stack (~2GB uncompressed)

Spatially decimated 8bit version (~35MB uncompressed)

Your proposed solution - calculating a sinusoid in the time domain based on the peak in the FFT, then subtracting it - should work, but there's an easier way to do essentially the same thing: modify that peak value in the FFT, then take the inverse transform.

So, for your rasterized video M[nframes, npixels] you find the frequency bin holding the artifact, then systematically flatten it (e.g., set its magnitude to the average of its neighbors) for every pixel:

import numpy as np
nframes, npixels = np.shape(M)
# Identify the bin containing the sinusoidal artifact
# Use the average intensity for each image
m = np.mean(M, axis=1)
# Calculate the FFT
Fm = np.fft.rfft(m)
# Find the largest bin away from the low-frequency region
lowfreq = 100  # or something
badbin = lowfreq + np.argmax(Fm[lowfreq:]**2)

# Now adjust the amplitude of that bin in the FFT of each pixel
for pixel in range(npixels):
   Fpix = np.fft.rfft(M[:, pixel])
   # Scale magnitude of artifact bin to be the mean of its neighbors
   Fpix[badbin] *= np.mean(np.absolute(Fpix[[badbin-1, badbin+1]]))/np.absolute(Fpix[badbin])
   # Rewrite the time sequence of that pixel
   M[:,pixel] = np.fft.irfft(Fpix)

This should work if the artifact is exactly constant amplitude and frequency, and its frequency falls right on to a submultiple of the sequence length (i.e., the sinusoids represented by the FFT).  In general, you might want to flatten out one or two bins either side of badbin to deal with a slightly broader set of narrowband corruptions, e.g.

# ...

   # Scale magnitude of artifact binS to be the mean of neighbors
   spread = 3  # flatten bins from (badbin - (spread-1)) to (badbin + (spread-1))
   # target value for new bins
   targetmag = np.mean(np.absolute(Fpix[[badbin-spread, badbin+spread]]))
   bins = range(badbin - (spread-1), badbin + spread)
   Fpix[bins] *= targetmag/np.abs(Fpix[bins])
   # ...

If you want to constrain the component removed from each pixel to have the same frequency and phase of the artifact detected in the mean intensity, you could remove just the projection of the badbin magnitude onto that phase, e.g.

badbinphase = np.angle(Fm[badbin])
# ...

   Ncomponent = np.abs(Fpix[badbin])*np.cos(np.angle(Fpix[badbin]) - badbinphase)
   Fpix[badbin] -= Ncomponent * np.exp(0+1j * badbinphase)
   # ...

Note that the resulting component at badbin will now always be 90\deg phase shifted (orthogonal) from the global badbinphase in every pixel - any signal component at exactly that frequency and phase cannot be separated from the artifact.


Cross-correlation of time compressed 1D signal If I ave 2 otherwise identical 1D signals that are phase shifted from each other then cross-correlation is a perfect way to identify the time lag between the 2 signals.

Now supposing I double the speed of one of those signals (eg sampled at every other sample).  The cross correlation will fail miserably at this point.

However, in image we can account for this kind of scaling (or even rotation) by converting the images to log polar coordinates and then performing the [phase] correlation.  This has the added advantage of being able to identify the angle of rotation and amount of scaling, I believe.

Is there a similar, simple transform for 1D signals to allow for identification of the time lag as well as identifying the amount of time compression of the signal?

I guess you could use log polar coordinates directly on the 1D signal but it would, instantly, square the memory requirements.  Any other thoughts?

With known time scaling, you can resample to a common rate and then cross-correlate.  It is natural to think you might need to go to the higher of the two rates, but that does not increase information beyond the lower rate, it only interpolates the information. This is analogous to zero-padding an FFT to interpolate the resolution of a peak.

To determine unknown time scaling, cepstral processing may provide some insights.  I'm afraid it is not my area, but it sounds related to what you are describing.

Essentially, you've got two dimensions of uncertainty (start time, rate) across which you'd like to search to find the best correlation.

You can reduce it to one dimension by considering the magnitude of the two signals' spectral data. The lower rate data should then look like a stretched version of the higher rate data.


How to recognise different symbols from their slope ? like (,^) etc i am trying to build a program in which user input a symbol with mouse and program recognise the symbol by saving mouse coordinates in array ,for this i have followed following steps:

Firstly check if line is straight or not by straight line equation.
I am able to recognise symbols like vertical,horizontal lines, forward and backward slash with the help of the slope of line.
but for >,<,^ i don't know what to do, for starters i have solved it like this, for example if < is to be recognised first symbol should be backward and second to be forward slash but this conflict with ^.

please let me know if there's another method for solving this problem.

My way of doing it would be:

Extract LBP feature vectors for the reference symbols are store it.
Now extract LBP feature vector of the test symbol.
Compare it with the list of available reference feature vectors and choose the one with minimum difference.

This method can also be extended to other symbols.


Relation between erfc and the BER Antipodal Signalling has a bit error probability of $P_b = \frac{1}{2} \mathrm{erfc}\big(\sqrt{\frac{E_s}{N_0}}\big)$ and on-off-signalling has a bit error probability of $\frac{1}{2} \mathrm{erfc}\big(\sqrt{\frac{E_s}{2 N_0}}\big)$. Why does this correspond to a higher SNR of 3 dB which we need for on-off-signalling to achieve the same bit error probability as with antipodal signalling?

I know that the factor 2 corresponds to 3 dB. But the complementary error function does not decrease by a factor of 2 but by a factor of $\sqrt{2}$

Assuming your bit error probability formulas are correct, for antipodal signalling

$$P_{b, aS} = \frac{1}{2} \mathrm{erfc}\left(\sqrt{\frac{E_{s,aS}}{N_0}}\right)$$

and for on-off signalling

$$P_{b, oS} = \frac{1}{2} \mathrm{erfc}\left(\sqrt{\frac{E_{s,oS}}{2N_0}}\right)$$

If these two error probability are identical, i.e. $P_{b, aS} = P_{b, oS}$,

$$\frac{1}{2} \mathrm{erfc}\left(\sqrt{\frac{E_{s,aS}}{N_0}}\right) = \frac{1}{2} \mathrm{erfc}\left(\sqrt{\frac{E_{s,oS}}{2N_0}}\right) \tag{1}$$

Because both $\mathrm{erfc(.)}$ and $\sqrt{.}$ are bijective functions, $(1)$ implies $\frac{E_{s,oS}}{N_0} = 2\frac{E_{s,aS}}{N_0}$ in linear scale and 3dB difference in log scale.


Palmprint Identification - Why do we align the images before we use the Fourier Transform? I was reading the paper PALMPRINT IDENTIFICATION BY FOURIER TRANSFORM by WENXIN LI, DAVID ZHANG and ZHUOQUN XU about identifying persons based on an image of their palm.

The first step in their identification algorithm consists of aligning and pre-processing the image since they spent a few pages on that front, this gave me the impression, that this step is kinda important.
What I didn't understand however: Why do we do this? The essential part in my opinion is the Fourier Transform afterwards, that is used to extract/generate the feature. But from my understanding, this should even work without this aligning/pre-processing step.
I have found a property of the DFT that states that rotating the spatial domain contents rotates the frequency domain contents, can be read here.
In order to extract features, some rotation might be necessary, at some point. In my opinion however, this does not justify to think up such an elaborate method of aligning/preprocessing the image, there has to be more behind that.
Therefore, I am most likely overlooking/missing something. Thanks in advance.

From the article:

  [...]only a sub area of a palmprint image is used in feature extraction.

The alignment by a shift and other image transformations is done in the full hand image which is then cropped to a sub image containing only the features of interest. To analyze the sub image they Fourier transform it. They are not shifting the sub image alone but the full image so the shift operation does not have a direct equivalent operation in the frequency domain of just the sub image.


Compute minimum phase version of a FIR I am working with HRIR filters, in particular I am trying to interpolate them.
One commod method in the literature to perform interpolation of HRIR is to use the minimum-phase decomposition and interpolate separately the minimum-phase part of the filter and then the all pass part.
You can find more here: https://www.researchgate.net/publication/277879431_On_the_Minimum-Phase_Nature_of_Head-Related_Transfer_Functions
What I am trying to do is to convert my arbitrary HRIR FIR into a minimum phase one. In particular, I am working in Python and I have been using this function.
However while performing the conversion from the original FIR into the minimum-phase one I get the following error:

RuntimeWarning: h does not appear to by symmetric, conversion may fail

With the following code:
min_phase_HRIR_0 = minimum_phase(hrtf[0], 'hilbert')

Where hrtf[0] is a 256 long FIR extracted from a HRIR database (in particular the HUTUBS database).
If I plot the magnitude  of my original FIR and the minimum-phase one I get the following plot:

Which is obviously wrong since I am expecting the min phase and the original FIR to have the same magnitude spectrum.
I guess there is something wrong with the Scipy function. Is there anything for Python that converts a FIR into a minimum-phase version?

In the end I used the Hilbert transform to compute the minimum-phase phase of my FIR.
I have used the following formula:

Where H is the spectrum of the desired filter.
Which has been used in:

Individualisation d’indices acoustiques pour la synthèse
binaurale
Phase Unwrapping for Spherical
Interpolation of Head-Related
Transfer Functions
Transformation characteristics of the external human ear

However I am not sure of the results I got.
A follow up question discussing the results and implementation can be found here: minimum-phase phase via Hilbert transform returned values


Loss in converting YUV to PNG I have a YUV 420NV12 image. I'm using ffmpeg to convert the YUV file to PNG and then back again to YUV.

ffmpeg -s 640x480 -pix_fmt nv12 -i testold.yuv test-640x480 test.png
ffmpeg -i test.png -s 640x480 -pix_fmt nv12 testnew.yuv

When I compare the two YUV files, there are few differences. I expected PNG to be a lossless compression. What is the reason for difference in two files?

The PNG format is lossless for RGB24 image data.  However the conversion from YUV to RGB24 is not lossless, as the two formats quantize the color space differently.

To see this, the following animated gif was made by applying your two ffmpeg operations 200 times back-and-forth and collating the resulting 200 images into the gif.

By contrast, the following gif was generated using yuv444p as the YUV format instead of nv12. This format does not do subsampling / as much quantization.


Root Locus, Transfer Functions and Unit Step Response? Consider the generic feedback loop, and the transfer function $G(s)$ shown by the following root locus plot.

Where $\mathbf{x}$ denotes the open-loop poles and $\square$ denotes the closed loop poles.
I want to determine if this root-locus produces any of the following output responses to a unit step reference signal:

Attempt:
I attempted to build the open loop transfer function:
$$G_{OL} = \frac{K}{s^2+2s+2}$$
Find the value $K$ using the FVT for the step response and graph a)
$$\lim_{s\to 0} = s\cdot\frac{1}{s}\frac{K}{s^2+2s+2} = 1$$
which gives $K = 2$, plug this back in and see if the step response matches any of the figures. It is possible none of the figures match, but I think I am doing something incorrect.
Root locus shows all the possibilities of K, but I am also given the closed loop poles. Can I use those to get K? Is the unit step response in relation to the Closed Loop Transfer Function? or the open loop? I can see that building the connection between the graphs is important, I seem to be missing something though. Thanks!
Edit/Update
I believe the corresponding Step Response Graph is C) but I'm still uncertain on if I followed the right process: Finding the closed-loop transfer function I get
$$G_{CL} = \frac{K}{s^2+2s + 2 + K}$$
Applying FVT to graph C I get $K = 3$, which makes
$$G_{CL} = \frac{3}{s^2+2s + 5}$$
and plotting the unit step response with Matlab gives C). What if I didn't have access to Matlab? Is there an easier way? Is it sufficient to match the closed loop poles after I determine $K$?

The closed loop poles are the roots of the polynomial
$$D(s)=s^2+2s+2+K\tag{1}$$
and, according to the root locus plot, they are $s_{1,2}=-1\pm 2j$. Consequently, we get
$$2+K=|1+2j|^2=5\quad\Longrightarrow\quad K=3\tag{2}$$
With $K=3$ we obtain $$H(0)=\frac{K}{2+K}=\frac35\tag{3}$$
which leaves step response $C$ as the only option.


Book recommendations on DSP, statistical SP and signal detection I am looking for a good reference on DSP. I am a little bit more than beginner and have some experience in DPS. I think I would have to start with rather easier books where a concept of DSP is explained based on practical examples, intuitively.

Eventually I'd take more advanced books to have solid understanding of DSP, including statistical signal analysis such as power spectra, cross-spectra, coherence, autocorrelation, and cross-correlation and ultimately I'd be interested in signal, sound, noise, or speech detection.

For beginners, people recommend

Rick Lyons' Understanding DSP
Sanjit K. Mitra, Digital Signal Processing: A Computer-Based Approach
The Scientist and Engineer's Guide to Digital Signal Processing - free online book

All recommendations of advanced books are by Alan V. Oppenheim's e.g.:

Alan V. Oppenheim, Ronald W. Schafer, John R. Buck, Discrete-Time Signal Processing, Prentice-Hall Signal Processing Series.
Alan V. Oppenheim, Alan S. Willsky, with S. Hamid, Signals and Systems, Prentice Hall

There was one more recommendation:

John G. Proakis, Dimitris K Manolakis, Digital Signal Processing: Principles, Algorithms and Applications, Prentice Hall

These books are usually quite expensive which is fine but I would not like to buy something I'd be disappointed with. If you could recommend or say whether the books I listed are worth the price, I'd be grateful :)

Thanks.

EDIT.

I have decided to buy and go through "Understanding Digital Signal Processing" by Rick Lyons. I am still waiting for it but have had a chance to read through first two chapters via Kindle App. I must say I am already happy and don't regret this purchase. I'd recommend this book for all starters.

The Proakis and Manolakis book is good if you're looking for one book.

If you're looking for depth about statistical signal processing, I recommend the series of three by Steven M. Kay:

Fundamentals of Statistical Signal Processing Volume I: Estimation Theory
Fundamentals of Statistical Signal Processing Volume II: Detection Theory
Fundamentals of Statistical Signal Processing Volume III: Practical Algorithm Development

If you want some more advanced statistical signal processing books, then the ones by Harry L. van Trees might be more your style:

Detection, Estimation, and Modulation Theory (Part I)
Nonlinear Modulation Theory (Detection, Estimation, and Modulation Theory, Part II)
Detection, Estimation, and Modulation Theory, Radar-Sonar Signal Processing and Gaussian Signals in Noise (Part III)
Detection, Estimation, and Modulation Theory, Optimum Array Processing (Part IV)


Which anti aliasing filter algorithm for efficient audio resampling? I'm building a real-time audio resampler (think pitch bend) that needs to have a several different performance vs quality configuration options. I understand that I'll need to apply a low-pass filter before resampling in order to avoid aliasing.

The requirements for the first filter are this:

very efficient
better than nothing

For this filter, I don't mind if I have a very slow roll-off, or if there's significant phase distortion. I just want it to be as fast as possible, and to attenuate frequencies above nyquist to some degree.

The requirements for the second filter are:

reasonably efficient
significant reduction in frequencies above nyquist

I'm guessing that I'll want to go with one of the classic algorithms I hear about for the second filter (chebyshev, butterworth), and potentially something like a recursive moving average filter for the first.

I realize my question may so vague as to be useless ("very efficient" and "reasonably efficient" are obviously relative terms!) I realize that filter selection and design is an art that one can spend a lifetime on. I'm quite inexperienced with DSP and am just looking for a couple algorithms I can get my hands dirty with. I'm hoping someone with audio experience can pip up and say "here are my general purpose, go-to algorithms for efficient lowpass filtering".

For a first cut, I would consider a windowed Sinc interpolation kernel for both resamplers.  Very narrow window with a small polyphase table for the fast one.  Much wider window with more phases in a larger, possibly interpolated lookup table for the slower higher quality filtered interpolator/resampler.  Lots of parameters (including choice of window) that can be adjusted for both quality and performance to meet widely varying requirements.


How Do You Save / Convert an Image to Be in the Frequency Domain? Quick Note: I don't believe this is a subjective question... however if it is I will gladly modify to be more objective

In ImageJ, if I take the FFT of an image, how can I save this new frequency domain image? The goal being that I would like to be able to open it up and take the IFFT.

When I try this, imageJ says you need to have a frequency domain image. This leaves me thinking that whatever information it has about the frequency domain is lost once I save.

My end goal is to take images that I have taken the FFT of from a python script that I am debugging and load them up in imageJ and take the IFFT.

A trick would be to use the TIFF format which supports many channels.
For each channel in your original image create 2 channels, one for the magnitude and the other for the phase.

TIFF supports 32 bit hence data will be well preserved.


When modelling noise,how do they decide whether additive or multiplicative? How does one take a decision of whether noise is additive or multiplicative?

It will depend on what you are modelling in your system. Consider an active radar system. There are both additive noise components and multiplicative noise components. What distinguishes them in this case:

If I increase the transmit pulse power by a factor of $K$ and the SNR improves by a factor of $K$ then the noise is additive
If I increase the transmit pulse power by a factor of $K$ and the SNR does not change then the noise is multiplicative

To give some examples to clarify:

Thermal noise at the input of the ADC is generally additive - noise level does not generally change with power of the received signal, so SNR will change with transmitted/received power.
Clipping in the ADC (distortion) is non-linear - it is more multiplicative than additive. Increasing the receive power won't affect the SNR (assuming clipping still occurs).
Phase noise in the ADC (timing errors) is multiplicative
Speckle noise - Presence of a lot of random scatters in radar resolution cell is usually modeled as multiplicative - increasing the transmit power does not reduce the speckle effect.

Another common example is in active Sonar. In near range there is a lot of reverberation (volume reverberation), but as the transmitted pulse travels away from the sonar the reverberation drops and the system performance becomes limited by the ambient noise. In these two cases the sonar is operating in reverberation limited mode and noise limited modes respectively.

Reverberation - Multiplicative - Increasing transmit power offers no improvement in SNR.
Noise Limited - Additive noise - Increasing the transmit power improves the SNR.

This assumes a very simplistic model of underwater sound propagation. You will have other wave propagation effects mostly due to the sound velocity profile (Reverberation at surface/bottom, sound channels, convergence zones, etc).

So within a single system you can have both. In what mode of operation will determine which type of noise dominates. It should also be noted that the noise may not be entirely multiplicative or entirely additive. You can have noise which contains elements of both. In this situation you will get some SNR improvement but it may not be linearly related to the increase in transmit power.


Polyphase components Spectrum formula in Discrete time I am trying to understand the meaning behind some of the components of the formula for polyphase components in discrete time. Note that $\Omega$ indicates discrete in this notation and the formula is given by:

Polyphase components is dividing the spectrum into blocks of size $L$ and each block divided into $k$ parts. If we remove the exponentials $e$ we have the formula for downsampling which makes sense since we are getting rid of the other polyphase components. The first exponential $e^{jk \Omega/L}$ seems to be the shift to move the index $k$ of the block of size $L$. However, I do not know what is the purpose of the second exponential $e^{-2\pi jkp/L}$

Thank you

Let's start with your equation.

$$ P_k\left( e^{j\Omega} \right) = \frac{1}{L} e^{jk\Omega/L} \sum_{p=0}^{L-1} e^{-2 \pi j k p / L} H\left( e^{j(\Omega-2\pi  p)/L} \right) $$

Do a little rearranging.

$$ P_k\left( e^{j\Omega} \right) =  e^{j\Omega\frac{k}{L}} \frac{1}{L} \left[ \sum_{p=0}^{L-1} H\left( e^{j\Omega\frac{1}{L} }  e^{-j2\pi \frac{p}{L}} \right) \left( e^{j2\pi/L} \right)^{-kp}   \right] $$

It can be made simpler to understand by a few variable substitutions.

$$ y = e^{j\Omega} $$

$y$ is a complex value on the complex unit circle at $\Omega$ radians around the circumference.

$$ h[p] = H\left( y^{\frac{1}{L}} e^{-j2\pi \frac{p}{L}} \right) $$

$h[]$ is a sequence of values sampled from the $H()$ function.  $y^{\frac{1}{L}}$ is a complex value on the unit circle at $1/L$th the angle of $y$.  The $e^{-j2\pi \frac{p}{L}}$ is a clockwise walk around the unit circle in $L$ even steps.  Conceptually:

$$ e^{-j2\pi \frac{p}{L}} = \left( e^{j2\pi } \right)^{-\frac{p}{L}} = 1^{-\frac{p}{L}} $$

Therefore $h[]$ has a period of $L$.  If $p$ were to go below zero or above $L-1$, it would be a repeat pattern.

With these substitutions, your equation becomes this:

$$ P_k\left( y \right) =  y^{\frac{k}{L}} \frac{1}{L} \left[ \sum_{p=0}^{L-1} h[p] \left( e^{j2\pi/L} \right)^{-pk} \right] $$

The part in brackets is the definition of the DFT.

$$ X[k] = \sum_{n=0}^{N-1} x[n] \left(e^{i2\pi/N}\right)^{-nk} $$

It's arranged the same way as in my answer here Fractional powers of complex numbers (DSPrelated computation)

The $\frac{1}{L}$ is a normalizing factor making the magnitudes of the DFT calculation independent of $L$.

$ y^{\frac{k}{L}} $ is a complex value $ k/L $ of the way to $y$ along the unit circle which means the results of the DFT bin are rotated that far.

So, a long answer to your question.  The purpose of the second exponential $ e^{-2\pi j k p/L} $ is to execute the DFT transform.  It is part of the definition and the reason it works the way it does.  Lots of stuff available on that if you are unfamiliar.


Question regarding using cardinality in pattern recognition So I've read this article for a project im working on:

Ortiz-Cataln: Cardinality as a highly descriptive
feature in myoelectric pattern
recognition for decoding motor
volition, Chalmers Uni, published in "Frontiers in Neuroscience", 2015

They use cardinality in pattern recognition and claim that it is "highly descriptive".

I can't understand how they have used it or why it would be of any use, if my understanding of cardinality as simply a measure of the size of dataset is correct. Can anyone explain to me what they mean and how they have used cardinality?

There's a lot of red flags in that paper that shouldn't have passed review. In fact, I'm a bit surprised – Mr. Ortiz-Catalan is a well-published researcher in a research group called "signals and systems" at a very prestigious university, and yet there's sentences like:

  Altering the dimension of the unit used for sampling (ADC
  resolution) to a high precision unit (for example, a double) would
  alienate the discrimination power of cardinality.

There's no ADC that give double precision floating point values. (Aside from "double precision" being a relative term, only defined within one programming language or CPU architecture). At any rate, the "D" in ADC stands for "digital", and that means "discrete in time and value, so

  This is because
  every sample value would be unique, so cardinality will always be
  equal to the number of samples in the time window (maximum
  value of cardinality).

is class-A hogwash, because sampling leads to discrete samples, and thus, the maximum number of unique values you can get is the cardinality of the alphabet of your ADC.

The same applies to:

  A simple operation like dividing the EMG sampled values
  by the gain of the amplifiers, in an attempt to use the original
  EMG amplitude in further computations, could increase the unit
  precision automatically by the processing software; that is, from
  the original acquisition values (for example, 14 bits) to a double,
  in the case of MATLAB (Massachusetts, United States).

For a fixed gain, you get as many unique values out of that division as you put in – division by a constant is a function (school math), for the same input value it should give the same output value; the data type doesn't matter here (it might, in a corner case, mean that the number of output values is smaller than the number of input values. Never more. And generally: Information can't be generated during signal processing – they cite Information Theory as basis for their research, so they must know the Data Processing Theorem, and gaining values would be equivalent to generating information).

I have to credit the author for the "I cite a software tool as coming from Massachusetts" joke – I liked that, it was cute :)

Having cleared up that ADCs only give discrete values, let's move forward.

Every ADC I encountered in EMG applications so far is used with a proper anti-aliasing filter up front. Needing to do that is mathematical/physical fact – if you don't band-limit the analog signal before sampling it, you get artifacts, namely aliasing in your signal, and also, especially for "peaky" signals such as EMG tends to be, instead of noticing every peak, you'd, in the best case, just get a non-1 probability of observation.

So here we go, feeding in band-limited signals into an ADC with a finite set of values.

Now, I mainly work in communication technology, so I'm much more familiar with the mathematical tools meant for analysis of harmonic signals, but we can still work with this pretty well:

Nyquist's sampling theorem says that, if the original signal is sufficiently band-limited, we can reconstruct the original signal from a periodical observation of instantaneous values. That periodicality is the sampling rate, and in the real-valued case, that needs to be just a bit higher than twice the highest frequency in the signal (i.e. the signal's bandwidth).

Nyquist doesn't require any synchronization between the signal and the sampling.

Thus, let's do the following thought experiment:

Sample a signal with bandwidth $B$, at times $t_i = i T_s$, $T_s < \frac1{2B}$.
Following Nyquist, we now have a complete representation of the signal
We find that if we configured our sampling system optimally, the extreme values of the sampled signals are at least close to the extreme values that the ADC can produce – otherwise, we'd be wasting dynamic range and could very well better represent the signal by amplifying it a bit before, reducing the quantization errors. This should have tipped us of that cardinality is extremely sensitive to the physical voltage amplitude – and thus a very hard-to-control meausure, considering that depends on a multitude of factors like electrode placement, contact quality, temperature of the system...
Now, the information how many different values where taken by the signal is identical to how many were NOT taken – since we know how many values our ADC can produce.
The observational windows were 200ms @ 2000Hz sampling rate, i.e. 400 samples. Hence, there were exactly 400 different cardinalities that could be realized here: 1 unique value (constant) or up to 400 different values. We can, with a high degree of certainty, say that certain cardinalities are more likely than others – for example, the constant signal is very unlikely to happen, and 400 different values isn't much more likely, either, since the signal is discrete and supposed to be correlated.
The paper finds that decreasing the sampling rate to $\frac12$ of the orignal 2 kHz only reduces the quality by 0.2 % – that is surprising, as half of the signal goes missing! We must pick one of three hypotheses:

A) the signal digitized was sufficiently bandlimited in the first place for the lower sampling rate, and in that case oversampling was done, which should reduce the meaningfulness of value set cardinality (as oversampled signals will just take more intermediate signals, thus shifting the density functions for that measure closer together higher up in of the range of possible cardinalities). The interesting part is here that additive noise will effectively lead to dithering – ie. an increase of cardinality, regardless of the original signal.
B) the signal was undersampled when it was sampled with 1 kHz, by a factor of $u = \frac{2B}{f_\text{sample}}$.
C) Most of the energy in the signal was in the lower half (or, in fact, the lower fourth, as 500 Hz sampling rate still seems to work nearly as good – which at some point should make one wonder what we're really looking at), but not quite all – and thus, the fast changes that (almost) always lead to unique values are rare and correlate highly with the classification.

Intuitively, C) works quite nicely: anything that changed by more than 1 quantization step between two successive samples at the 1 kHz rate will almost always have one intermediate value (excluding strong noise for now) sampled at 2 kHz. But: Knowing typical EMG spectra, there's less enery at higher frequencies (i.e. "faster jumps"), and thus, it's very unlikely that strong peaks and intermediate values from these will be duplicate within an observation window (especially if that is only 400 samples long).
Thus, the counting of unique values is but an aliasing detector – or, if you want so, a detector for energy in the upper half of the spectrum. A high-pass filter!

I personally have my problem with papers that offer $p$-values but don't even mention how they calculated those. Maybe you could nicely ask to get their algorithms/software. This should be no problem – the paper boldly boasts "open access" insignia, and that ideologically implies that they're proud to share not only their results but also their methodology. I do believe that the paper is on to something – it seems to have demonstrated a classifier that outperforms others. However, these others might have been badly parameterized, and might not perform that much worse, especially on properly preprocessed data.


What are some practical examples of complex signals? I am aware that complex discrete-time signals are of the form:
$$
x(n) = a(n) \cdot e^{j\omega n T}
$$
However, I have not yet worked with complex signals, only real ones. What are some practical examples of complex signals? Answers with plenty of detail are appreciated.

Here is one aspect that often creates confusion: In many cases, complex signals are used solely as a mathematical convenience, not because the complex representation is required by anything physical.
A good example is circuit analysis in electronics: all physical quantities (voltages, currents, impedances) are real and you can do the analysis with these quantities and the associated physical equations directly.
However, it's a lot easier if you convert the real quantities to complex ones, run the analysis in the complex domain and convert the results back to real. Hence, that's the way everyone does it.
So you don't HAVE to use complex signals, it just gets a lot easier if you do.


You need to do partial fraction expansion.
Most of the time, the result will be 3 parallel sections with the same poles as the original filter but only with a real zero in each parallel section. The first section is just a single real scale factor.
$$H(z) = r0_0 + \frac{r0_1 + r1_1 \cdot z^{-1}}{a0_1 + a1_1 \cdot z^{-1} + a2_1 \cdot z^{-2}}  + \frac{r0_2 + r1_2 \cdot z^{-1}}{a0_2 + a1_2 \cdot z^{-1} + a2_2 \cdot z^{-2}} $$
Note that the "a" coefficients are the same as in the cascaded form, i.e. you only need to calculate the "r" coefficients.
This is all relatively straight forward if your poles are unique and complex conjugated pairs. Repeated poles and single real poles require special treatment.
The Matlab function residuez() can typically help with this although it can have numerical problems if the poles are too close to the unit circle. If you need to do it by hand take a look at https://ccrma.stanford.edu/~jos/filters/Partial_Fraction_Expansion.html
Since you only have two sections, the following approach should work

Write both versions (parallel and serial) as a single fraction
Discard the denominators, they are the same
Both numerators are a 4th order polynomial in $z^{-1}$, i.e. they have 5 coefficients.
Each power of $z^{-1}$ needs to have the same coefficient in both numerators.
That gives you five linear equations with 5 unknowns that you can solve for $r0_0, r0_1, r1_1, r0_2, r1_2$


Which image filter can be applied to remove gridded pattern from corrupt jpegs? I have around 1,400 jpegs that have been corrupted somehow and have lost the backup images. They all seem to have the same gridded pattern of lines over each (i.e. the gridding does not shift from image to image.

Here's what one of these images looks like:

Are there any image filtering techniques in Matlab particulary or otherwise that will remove or smooth this gridding pattern?

You can use a standard inpainting algorithm. These algorithms replace marked pixels in an image with the pixel values that surround these marked pixels. The challenge here is to detect the grid (my tests seem to show that it is not a completely regular grid). So, I came up with this solution:

from PIL import Image
import requests
from io import BytesIO
import cv2

url = "http://i.stack.imgur.com/Ahrnl.jpg"
response = requests.get(url)
img = Image.open(BytesIO(response.content))

plt.imshow(img)
A = np.array(img)
A2 = A.copy()
A_gray = cv2.cvtColor(A, cv2.COLOR_RGB2GRAY)

# Do some rough edge detection to find the grid
sX = cv2.Sobel(A_gray, cv2.CV_64F, 1, 0, ksize=3)
sY = cv2.Sobel(A_gray, cv2.CV_64F, 0, 1, ksize=3)
sX[sX<0] = 0
sY[sY<0] = 0

plt.subplot(221)
plt.imshow(sX)

plt.subplot(222)
plt.imshow(sY)

plt.subplot(223)
# the sum operation projects the edges to the X or Y-axis.
# The 0.2 damps the high peaks a little
eX = (sX**.2).sum(axis=0)
eX = np.roll(eX, -1) # correct for the 1-pixel offset due to Sobel filtering
plt.plot(eX)

plt.subplot(224)
eY = (sY**.2).sum(axis=1)
eY = np.roll(eY, -1)
plt.plot(eY)

mask = np.zeros(A2.shape[:2], dtype=np.uint8)
mask[eY>480,:] = 1
mask[:, eX>390] = 1

A2[mask.astype(bool),:] = 255
plt.figure()
plt.subplot(221)
plt.imshow(A)

plt.subplot(222)
plt.imshow((A2))

restored = cv2.inpaint(A, mask, 1, cv2.INPAINT_NS)
plt.subplot(223)
plt.imshow(restored)

The program output is as follows:

To detect the grid I did a quick-and-dirty solution. It can be improved a lot, but it shows the initial idea. The general flow is:

detect the grid
create a mask that describes which pixels are corrupted by the grid
inpaint the corrupted pixels.

For inpainting I used OpenCV inpaint operation. For detecting the grid, I performed edge detection in X and Y direction using a Sobel filter. Then I add all edge values in the X-direction and Y-direction to find peaks, where the grid lines are. Then, I choose the highest peaks as the coordinates where the grid lines are estimated. It's not working perfect (e.g. strong edges in the image are falsely detected as grid lines), but it shows the idea. It can be improved by e.g. Hough transformation to find lines, kicking out very strong edges etc.

Alternatively, if the grid is really the same for all images, then you can perform the grid detection jointly for all images, which would yield a much better accuracy (just do the technique above, but before choosing the peaks, sum up the results from all pictures). In more detail, you would calculate eX for all images and add all these eX together into a single vector. This vector will have a much clearer peak structure and the thresholding can be done easier.


describing a block diagram // How do I describe multiplying a signal through multiple branches formally? Let's say I made this block diagram and I want to explain it:

FYI: $x$ is a signal and each $y$ box is a matrix

I want to say that:

The signal $x$ is multiplied by each matrix $y$ in the different branches independently to produce
$f_0=xy_0$

$f_1=xy_1$

...

$f_m=xy_m$

I want to ask the experts in the signal processing world if this is formal enough? or am I sounding weird?

Edit 1:

What if I said:

The signal $x$ is multiplied by a set of matrices $y_k$ where $k = 0,1,..., m$

does it sound better? or would it imply that $x$ is multiplied like: $x*y_0*y_1*...*y_k$

which is NOT what I intend to say

What I would do is:

Add $m$ multiplier blocks, using a circle with a $\times$ inside
Each multiplier would have inputs $X$ and $Y_k$
Each multiplier would have output $f_k$
Be consistent in the use of bold typeface, uppercase and lowercase

Then I would say something like "The system has input $X$ and outputs $f_k=XY_k$, for $k = 0,\ldots,m$. In the diagram, $\times$ stands for signal-matrix multiplication, which is defined as follows...".


FIR Filter implementation without padding Is it possible to implement FIR filtering action without padding the input and coefficients?

i.e. Let's say if the input and filter coefficients are of size 4, then the output will be 7 samples. So, while implementing, we generally add 3 more zeros to both input and filter coefficients making them equal to output size.

But, if the input and filter coefficients are of size 1024, then the output will be of 2047 samples. So, now, we need to add 1023 zeros to both input and filter coefficients. This is inefficient, right?

So, I just want to know is there any other way to implement FIR Filtering without padding?

The below code gives the idea I was talking about.

int main()
{
   int x[4],h[4],y[7]={0};

   int i,j;

   for(i=0;i<4;i++)
   {
        x[i] = i+1;
        h[i] = i+1;
   }

   for(int i=0;i<7;i++)
   {
      for(int j=0;j<4;j++)
      {
         y[i+j] = y[i+j]+x[i]*h[j]; //filtered signal of length M+N-1
      }
   }

   for(i=0;i<7;i++)
     printf("%d\n", y[i]);

}

Output: 1,4,10,20,25,24,16 - expected

But obtained results - 1,4,10,20,garbage value, garbage value, garbage value

As I said in a comment, no padding is necessary. A simple FIR implementation might look like (it is assumed that the filter kernel h[] has already been implemented)

float y[N+M-1];
//set y[] to zero
for(int i=0;i<N;i++){ //N - signal length

 for(int j=0;j<M;j++){ //M - filter kernel length (usually much shorter than N

   y[i+j] = y[i+j]+x[i]*h[j]; //filtered signal of length M+N-1
   }
}

h[] represents the filter kernel and is usually implemented as a window-sync filter of length M (in your case it's just the original signal).

For your use case, you'd have

int x[7],h[7],y[13]={0};

       int i,j;

       for(i=0;i<7;i++)
       {
         if(i<4)
         {
            x[i] = i+1;
            h[i] = i+1;
         }
         if(i>=4)
         {
            x[i] = 0;
            h[i] = 0;
         }
       }

       for(int i=0;i<7;i++)
       {
          for(int j=0;j<7;j++)
          {
             y[i+j] = y[i+j]+x[i]*h[j]; //filtered signal of length M+N-1
          }
       }

Note that since you're convolving a signal with itself (x*h), you don't even need h[] but I kept it for clarity.


1D signal to 2D spectrogram I struggling at transforming a 1D breathing signal into 2D spectrograms,
power spectral density as a function of time and frequency.
Does anyone have an idea about that?

If you're working in the MATLAB environment, the STFT of a 1D signal can be obtained by using the Spectrogram() command. Refer here

And as for PSD, you can use the FFT to obtain the same. Refer here


What is the Sound of a Solution to the Wave Equation? What is the sound of a solution $y(x,t)$ to the Wave Equation

$\frac{\partial^2y}{\partial t^2}=c^2\,\frac{\partial^2y}{\partial x^2}\enspace$?

To be clear, I don't look for a solution of the Wave Equation, I have a concrete solution $y(x,t)$. I only want to describe the sound of this concrete solution.

In my case I want to describe the sound of a vibrating string, so the ends of the string are fixed, i.e. $y(0,0)=y(L,0)=0$ where $L$ is the length.

With a solution $y(x,t)$, I get perfectly the shape of the string between $x=0$ and $x=L$ for every time $t$.

I can show that $y(x,t)$ must be also $2L$-periodic according to $x$, i.e. $y(\lambda+2L,t)=y(\lambda,t)$, and $(2L/c)$-periodic according to $t$, i.e. $y(x,\lambda+2L/c)=y(x,\lambda)$.

To describe the sound, I would like to

find a frequency spectrum (or magnitude/amplitude spectrum) or
produce a sound file with a series of samples $x(i)$ to hear the vibrating string.

For 1. I don't know how to do that. Should I develop the Fourier-Series according to $x$ or to $t$ (a Fourier-Series exists, because I've showed the periodicity) or do complety other things? For a Fourier-Series according to $x$, the Fourier-Coefficients $a_k,b_k$ depends on $t$ and for a Fourier-Series according to $t$, the Fourier-Coefficients depends on $x$.

Any idea?

Edit:

Here is an example for a string that get plucked at $x_p$ with a displacement of $h$. Its shape is a triangle as follows:

$y_0(x)=
\begin{cases}
x(h/x_p)&\text{if}\enspace 0\leq x < x_p\\[1ex]
h\left(1-\frac{x-x_p}{L-x_p}\right)&\text{if}\enspace x_p \leq x \leq L
\end{cases}$

The initial velocity is 0 so finally the string is described as:
$y(x,t)=y(x,t)=\frac{1}{2}\left(y_0(x+ct)+y_0(x-ct)\right)$

In the following picture, you can see the process of the string starting at $t=0$

Since you already have the solution to the wave equation, what you need to do now is define a "pick-up" point somewhere along the length of the string (some $x$) and obtain its trace in time, i.e. look at what that specific $x$ value is "doing" across all the $t$ instances that you have. In other words, draw a straight line down your vibrating string diagram and take the vertical displacement of your string at each time point as the amplitude of your sound signal at that time point.

In reality, when you hear an electric instrument you hear exactly that, you hear what the vibrating string is doing exactly above the pick-up coil. Furthermore, electric guitars (and other instruments) have more than one pick-up points along the length of a string and a switch that is used to change the way these pick-up points are combined to offer different qualities of sound. For example, the coil closest to the neck (closest to the middle of the distance between the neck and bridge) produces sound with more low frequencies, while the coil closest to the bridge (the fixed point at the lower end of the instrument) produces sound with more high frequencies. The placement of the coils is a very important issue. Imagine for example what would happen if you placed the coil at a nodal point on the string (when it is excited to produce a standing wave).

Similarly, in acoustic instruments, the sound that reaches the ear is the projection of the vibrating instrument (notice now we are talking about the whole instrument) towards the direction of the ear or microphone. So, in effect, you are still applying this 'slice of vibrating element in time' concept. By the way, this is why mic placement is important when you record acoustic instruments (see for example: http://www.soundonsound.com/sos/1995_articles/sep95/mikinginstruments.html)

Hope this helps.


Wiener Khinchin theorem : struggle in the derivation I am reading Introduction to quantum noise, measurement and amplification, and I need to understand the Wiener Khinchin theorem: how to derive it. I also need to understand some context around this theorem (why some object are defined the way they are).
The theorem is derived on the page 55 of this document. I will explain to you what I understand from the starting point of this theorem and where I am stuck.

I assume a stationary stochastic gaussian process. I call $R_{XX}(\tau)=\langle X(\tau) X(0) \rangle$ the autocorrelation function of a stochastic variable $X$. I assume the signals to be real.
We define the power spectral density of noise as:
$$ S_{XX}(\omega) \equiv \lim_{T \to +\infty} \langle |X_T(\omega)|^2 \rangle $$
Where $X_T(\omega)$ is the windowed Fourier transforme of $X$ on a period $T$. Basically we define the signal $X_T(t)$ that is equal to $X(t)$ for $t \in [-T/2,T/2]$ and that vanishes elsewhere. The Fourier transform of this signal can thus be defined as:
$$ X_T(\omega) = \frac{1}{\sqrt{T}} \int_{-T/2}^{+T/2} dt X_T(t) e^{i \omega t} $$
To be consistant we must define the inverse Fourier transform as:
$$ X_T(t) = \frac{\sqrt{T}}{2 \pi} \int_{-\infty}^{+\infty} d \omega X_T(\omega) e^{-i \omega t} $$
First questions:

Why do we take such "uncommon" definition for the Fourier transform
(usually it is $\frac{1}{\sqrt{2 \pi}}$ for both sides
(direct+reciprocal), or $\frac{1}{2 \pi}$ in front of one of the
sides. Is there a theoretical motivation for this in this context ?
Why do we need to window the function ? Is it to avoid technical difficulties with divergent integrals ?

Now, I go to the derivation and express where I am stuck. The goal is to show that $R_{XX}(\omega)$, the Fourier transform of the autocorrelation function corresponds to the noise power spectral density $S_{XX}(\omega)$. I define $S^T_{XX}(\omega)$ the quantity before the limit of $T \rightarrow +\infty$. I know how we can derive computing the inverse Fourier transform of $S_{XX}(\omega)$ directly as suggested in a link in the comment. However I am stuck on deriving it from the other side (I have $S_{XX}(\omega)$ on the left hand side and I make appear $T.F(R_{XX})$ on the right handside). It might look as a boring question but in the document I follow they always follow this way of proving equalities, and I really don't get the manipulation they do. To understand how they do I need to understand this "other way" of deriving the theorem.
$$ S^T_{XX}(\omega)=\frac{1}{T} \int_{-T/2}^{+T/2} dt \int_{-T/2}^{+T/2} dt' \langle X_T(t) X_T(t') \rangle e^{i \omega (t-t')}$$
Now, I use the fact my process is stationnary. And I recognize the autocorrelation function:
I start from the definition (before taking the limit):
$$ S^T_{XX}(\omega)=\frac{1}{T} \int_{-T/2}^{+T/2} dt \int_{-T/2}^{+T/2} dt' R_{X_T X_T}(t-t') e^{i \omega (t-t')}$$
I do the change of variables: $(t,t') \rightarrow (\tau=t, \tau'=t-t')$. It gives me:
$$ S^T_{XX}(\omega)=\frac{1}{T} \int_{-T/2}^{+T/2} d \tau \int_{t-T/2}^{t+T/2} d \tau' R_{X_T X_T}(\tau) e^{i \omega \tau}$$
At this point, I have absolutely no idea how we could proceed.
We can also try another method (which I think is the way they do in the document but I really don't understand their justification). We can do the change of variable $(t,t') \rightarrow (v=t+t',u=t-t')$. The Jacobian gives $1/2$ and the boundaries are the one I write below:
$$S^T_{XX}(\omega)=\frac{1}{T} \int_{-T}^{T} dv \int_{-B(v)}^{+B(v)} \frac{du}{2} R_{X_T X_T}(u) e^{i \omega u}$$
$B(v)=v+T$ for $-T \leq v \leq 0$ and $B(v)=t-v$ for $0 \leq v \leq T$. It gives:
$$S^T_{XX}(\omega)=\frac{1}{T} \int_{0}^{T} dv \int_{-(T-v)}^{+(T-v)} \frac{du}{2} R_{X_T X_T}(u) e^{i \omega u}+\frac{1}{T} \int_{-T}^{0} dv \int_{-(T+v)}^{+(T+v)} \frac{du}{2} R_{X_T X_T}(u) e^{i \omega u}$$
And same problem here: I really don't see how we could get to the result with this way of deriving (which I think is roughly the path taken is the document I linked)
Thus, second question: How can prove the theorem with this method ?

Let me state at the beginning that the details that make this rigorous do not bring any extra understanding of the statistical behavior of signals, so the desire to use the Dirac delta distribution is completely understandable. If I could have figured out how to use it here rigorously, I would have.

I am going to repeatedly skip proof of the validity of "bringing expectation inside the integral". Note that computing an expectation is yet another integral (and over a finite measure space), so it really is just another change in the order of integration. To do so rigorously, appeal to Tonelli's Theorem and then Fubini's Theorem.

Let $x$ be a continuous-time wide-sense stationary stochastic process.
\begin{eqnarray}
r_{xx}(\tau) &=& \mathsf{E}\left[x(t)x(t+\tau)\right]\\
X_T(\omega) &=& \frac{1}{\sqrt{T}}\int_{-T/2}^{T/2}x(t)e^{-i\omega t}dt
\end{eqnarray}
We assume that the integral defining $X_T(\omega)$ converges for almost all realizations of the process $x$.
\begin{equation}
\begin{split}
\mathsf{E}\left[X_{T}(\omega)\overline{X_{T}(\omega)}\right] &=~
\mathsf{E}\left[\frac{1}{\sqrt{T}}\int_{-T/2}^{T/2}x(t)e^{-i\omega t}dt\overline{\frac{1}{\sqrt{T}}\int_{-T/2}^{T/2}x(t')e^{-i\omega t'}dt'}\right]\\
&=~ \mathsf{E}\left[\frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}x(t)x(t')e^{-i\omega(t-t')}dtdt'\right]\\
&=~ \frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}\mathsf{E}\left[x(t)x(t')\right]e^{-i\omega(t-t')}dtdt'\\
&=~ \frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}r_{xx}(t-t')e^{-i\omega(t-t')}dtdt'\\
&=~ \frac{1}{T}\int_{-T/2}^{T/2}\int_{-T/2}^{T/2}\color{red}{r_{xx}(t-t')e^{-i\omega(t-t')}}dtdt'
\end{split}
\end{equation}

Note the presence of $t-t'$ as the argument of factors in the integrand. Our final goal is to have the Fourier transform of $r_{xx}$ on the right-hand side. Perhaps we can convert this to an iterated integral in $t-t'$ and a different variable. Let
\begin{eqnarray}
p &=& t - t',\\
q &=& t + t'.
\end{eqnarray}
We want to change this to an iterated integral over $p$ and $q$. We are parametrizing the $(t,t')$ points with $p$ and $q$, and the Jacobian determinant of the map from $(p,q)$ to $(t,t')$ is $\frac{1}{2}$, so $dt~dt' = \frac{1}{2}dp~dq$.

The messy details that remain are describing the region of integration in terms of $p$ and $q$. The set over which we integrate is the diamond with corners at $(0,T)$, $(T,0)$, $(0,-T)$, and $(-T,0)$ in the $(p,q)$-plane. For $-T\leq p \leq 0$, we integrate $dq$ from $q = -T-p$ up to $q = T+p$. Hence the integral over $-T \leq p \leq 0$ is
\begin{equation}
\begin{split}
&~\frac{1}{T}\int_{-T}^{0}\left[\int_{-T-p}^{T+p}dq\right]r_{xx}(p)e^{-i\omega p}\frac{1}{2}dp\\
~=&~\frac{1}{T}\int_{-T}^{0}\left[T+p-(-T-p)\right]r_{xx}(p)e^{-i\omega p}\frac{1}{2}dp\\
~=&~ \frac{1}{T}\int_{-T}^{0}\left(2T + 2p\right)r_{xx}(p)e^{-i\omega p}\frac{1}{2}dp\\
~=&~ \int_{-T}^{0}\left(1 + \frac{p}{T}\right)r_{xx}(p)e^{-i\omega p}dp
\end{split}
\end{equation}

Now I add a restriction that a more seasoned professional might not need: I assume that
\begin{equation}
\int_{-\infty}^{\infty}\left|p~r_{xx}(p)\right|dp < \infty.
\end{equation}
Then
\begin{equation}
\begin{split}
&~\lim_{T\to\infty}\int_{-T}^{0}\left(1 + \frac{p}{T}\right)r_{xx}(p)e^{-i\omega p}dp\\
&=~\lim_{T\to\infty}\int_{-T}^{0}r_{xx}(p)e^{-i\omega p}dp
+ \lim_{T\to\infty}\frac{1}{T}\underbrace{\color{red}{\int_{-T}^{0}p~r_{xx}(p)e^{-i\omega p}dp}}_{\color{red}{\textrm{bounded as $T\to\infty$}}}\\
&=~ \int_{-\infty}^{0}r_{xx}(p)e^{-i\omega p}dp
\end{split}
\end{equation}

A similar calculation for $0\leq p \leq T$ (for which $-T+p\leq q\leq T-p$) yields
\begin{equation}
\int_{0}^{\infty}r_{xx}(p)e^{-i\omega p}dp.
\end{equation}
Adding these shows that
\begin{equation}
\lim_{T\to\infty}\mathsf{E}\left[\left|X_T(\omega)\right|^2\right]
~=~ \int_{-\infty}^{\infty}r_{xx}(p)e^{-i\omega p}dp.
\end{equation}


Embedded Image Processing System that sends photos to Smartphone via Bluetooth I would like to create an embedded system that takes a Full HD photo by pressing a button,processes and compresses it(for ex. to jpg) and then sends it via Bluetooth to the smartphone - where it's saved in the photo gallery.

I have already started to read "Digital Image processing" by Gonzalez, Woods, and I was thinking if, for my specific purpose, there isn't a better book with which to start.

I am interested to gain understanding about the overall process and the number of components needed for my task and how would they connect with each other.

Do you know any reference material(books/articles/experiments etc.) with which I could start?

I don't believe that I will find everything in one single book, so you can suggest different sources of information for each step needed in the process.

Thanks.

This question is more of a hardware setup. Books like Gonzalez will be more relevant when you have images as pixels and process them to transform to more desirable images.

Here are the few tips:

Setup a hardware like Raspberry PI - a small general purpose computer
Add hardware I/O interface like Camera Module
You can use standard USB Dongles and make it work on your pi.
Write your application which will now collect pics and send across to the app on your smart phone!

Thats' it, really!


Is it possible to reconstruct an image using its magnitude only? Let us take any image like Lena or cameraman. I know that when Fourier transform is applied on image,spatially global features like amplitude and phase are obtained locally. My question is whether it is possible to reconstruct an image from its magnitude only keeping phase value zero or different than actual value? Can anybody explain?

Re-construction using just the magnitude (and an assumed phase of zero) will only work for exactly symmetric images.  The phase is needed to make the top look different from the bottom and left look different from the right side of the image, etc.

That's because all the cosine waveform DFT basis vectors are exactly symmetric around the center of an FFT or DFT window.  Thus, unshifted, they can only reconstruct a symmetric result.  You could shift the phase of each basis vector randomly (or 90 degrees to make them all sinewaves) but, without knowing the phase shifts of the original image's DFT (or guessing right), you would just get another garbage result.


How do I determine whether I have pattern-free noise? For microscopy, we're frequently testing cameras. Since my applications involve very low signal-to-noise ratio, it becomes important that the noise is free of correlations and patterns, because local correlation is all that really distinguishes signal from background.

To test the noise, I usually acquire a series of ~100 dark frames, i.e. frames where no external light hits the camera, determine the fixed camera pattern by time-averaging, and subtract that from the series.

I've observed patterns in the noise by simply taking the standard deviation for each pixel through time and looking at the resulting image (where e.g. different rows/columns of the camera had different noise standard deviations), and by doing row and column-wise cross-correlation (where I noticed for some interleaved camera that the noise was correlated between every other row).

The first of these tests is qualitative only, and the second only gives me (relatively) global correlations. Are there better (and faster?) ways to determine whether there is any correlation or dynamic pattern in the noise of the camera?

If I were you, I would simply take a large number of no-signal, where you aren't measuring any real signal. Run it through whatever identification algorithms you have, and see if you see anything. If you do, then you need to worry about correlations.

I think what you might be missing is that correlation does not always mean a false detection, especially if you have a robust algorithm to this kind of noise.


How does a sound card translate amplitude values to sound when playing WAV/PCM/FLAC files? I'm reading the documentation for WAV files (and importing them in C++)It seems like 44100 amplitude (this can change, and is specified by the WAV header) values get sent to the sound card every second. What I don't understand is that these amplitudes can be negative. I'm confused as to what these amplitude values represent. How can they be negative and positive? What happens when you send an amplitude value to a sound card to produce a sound that you can hear?

They represent fluctuations from atmospheric pressure.

Positive values = increases in pressure = speaker cone moving outward
Negative values = decreases in pressure = speaker cone moving inward
Zero = no change in pressure = speaker cone in its relaxed center position = 1 atm

Since we can't hear a constant pressure, the constant pressure is averaged out and ignored, and only the relative vibrations are important.


Why does the excitation signal appear, separated, at high quefrencies in the cepstrum? So, I've just begun a speech and language processing course and have found the explanation of the process of getting the cepstrum of a signal and its properties a little confusing. The following is a description of my current understanding and an explanation of confusion it's causing me:

start with the speech signal. we can think of it as the formant signal convolved with a the excitation signal which is a dirac comb (approximately).
take the FFT, giving the spectrum of the excitations multiplied with the spectrum of the formants. the FFT of the excitation signal is another dirac comb with period 1/T
take logs. so the 2 signals above are now added
inverse fourier transform - the two signals from 1 should now be combined in addition (the FT transform is linear)

So if those 4 step are right, then why does the excitation appear at a particular region in the quefrency domain? it should emerge as a Dirac comb, added to the formant impulse response shouldn't it?

When using an FFT, an evenly spaced sequence of events in one domain usually produces a strong component in the other domain at a location related to the spacing of the events in the first domain.

A voiced speech signal usually includes a lot of harmonics which are evenly spaced in the frequency domain.  These evenly spaced events in the frequency domain will thus produce a component in the quefrency domain whose location is related to the spacing of the harmonics.  And the spacing of the harmonics in the frequency domain is related to the spacing of the exciter events in the time domain (the period of the fundamental).


The frequency spectrum of a static non-linearity driven by colored noise Does there exist any straight-forward theory that can explain what happens to colored noise that passes through a static non-linearity? That is, if you have colored noise
$$v = H u \, , \; u \sim N(0,\sigma)$$
generated by driving an LTI filter $H$ by white noise, and pass it though a static non-linear function $y = f(v)$, e.g. $f(v) = v^2$; what is the spectrum of $y$?

If your nonlinearity can be expressed as a polynomial (i.e., in terms of addition and multiplication), you can make use of:

The linearity of the Fourier transform, i.e., if $f$ and $g$ are (benign) functions, $a$ and $b$ are numbers and $ℱ$ denotes the Fourier transform, then:

$$ℱ(a·f+b·g) = a·ℱ(f) + b·ℱ(g)$$
The convolution theorem, which states that multiplication in the time domain corresponds to convolution in the frequency domain (and vice versa):

$$ℱ(f·g) = ℱ(f) ∗ ℱ(g)$$

So, in your specific example, the spectrum of $y$ is the spectrum of $v$ convoluted with itself:

$$ℱ(y) = ℱ(v·v) = ℱ(v) ∗ ℱ(v)$$


What are the typical causes of clicks and pops in an audio file? If a device generates clicks and pops at start and end of audio file playback, even when those audio files are completely zeroed out and dc-offset corrected, what could've caused them?

Clicking in an audio stream is typically indicative of some kind of discontinuity in the time-domain waveform. I'm not sure what you mean by "zeroed out and dc-offset corrected," but these sort of artifacts can be created via unsophisticated audio editing methods. If you're suggesting, though, that an all-zero (i.e. pure silence) audio file will cause this behavior, then there's probably some little implementation detail that's causing it. Examples:

Audible effects caused by turn-on/turn-off of the output audio amplifier. This is relatively common.
Significant differences between the initial/final sample values and the output level forced to the DAC when "idle" (I can't think of a sane reason why the designer wouldn't make the DAC output zero when not in use, but it's possible to do something else).


Wiener filtering/deconvolution for non-stationary noise Consider a stationary discrete-time random process $x[k]$ which undergoes low pass filtering by a filter with impulse response $h[k]$ and is subject to additive, temporally uncorrelated noise $n[k]$ of time-varying power $\sigma^2_n[k]$:
$$
y[k] = \sum^\infty_{m=0} h[m] x[k-m] + n[k].\qquad
$$
The time-variation in noise power is described by:
$$
 E[n^2[k]]=\sigma^2_n[k]= Ae^{k/B} \qquad B>>1
$$
such that the noise power increases (slowly) with time.
Specification of length $L$ FIR filter coefficients of a Wiener deconvolver seeking to estimate $x[k]$ from $y[k]$ is straight-forward for stationary noise.
Things are more complicated for non-stationary noise.
Strictly speaking non-stationary noise requires a different filter for every sample instant $k$ which is prohibitively complex.  In my case, the noise changes slowly enough for it to be modeled as locally stationary over significantly more than $L$ samples.  I could always just update the filter coefficients every $NL$ samples to reflect a new noise power assumption, but I'm concerned about abrupt transient responses due to the sudden jump in deconvolver coefficients.  Something more frequent but more gradual sounds better.
Is anyone aware of a smart recursive form of the Wiener filter for cases like this (especially where the noise power itself follows a recursion: $\sigma^2_n[k]=e^{1/B}\sigma^2_n[k-1]$)?

The Wiener Filter can written as a Least Squares problem (See How Is the Formula for the Wiener Deconvolution Derived?).
In you case, Since noise is basically with different STD per sample you may use the simple extension of Weighted Least Squares on top of the terms for the Wiener Filter.
Once you do that you'll have the optimal solution.


What is the difference between photon shot noise and salt and pepper noise? I am currently working on an application which deals with the restoration of fluorescent beads examined using widefield microscopy. There are two major sources of noises that are being encountered in the image: photon shot noise and salt and pepper noise. I have read in a paper that, both of them can be filtered using a median filter. So, is the essential difference between them only their sources?

photon shot noise follows a Poisson Distribution, which is essentially a shifted Normal Distribution. It is an electrical noise usually caused by charge surges. Salt and Pepper Noise is more due to digital conversion, and is usually modelled using a Fat-tailed distribution (which is essentially a shifted Normal Distribution). What to note is for any unaltered image, the order of magnitude of each type of noise, with photon shot noise usually much lower than salt and pepper noise (unless specifically created)

True both noise can be treated using median filter, and often people use salt and pepper noise to represent shot noise.


How to remove a tone/interference from a signal, without using classical linear filtering I am wondering if a methodology(ies) exist for essentially 'filtering' out a tone or interference, WITHOUT resorting to classical linear filtering techniques. (band pass filter, etc for example).

The reason for this question is because inherent in any classical linear filtering operation, I am at the end of the day linearly combining samples among themselves, and while this might do well to remove a tone of another frequency that I am trying to remove, it also has the effect of 'smearing' my signal in time. This smearing effect causes information such as the 'sharpness' of the signal to be lost, and in this specific case, the 'sharpness' carries information such as start of signal energy, etc, that I use later.

What method(s) exist (do they exist?) for this goal?

I should also add, please assume omnidirectional case, (I only have one sensor), but if assuming a case with multiple sensors lends itself to a solution, feel free to use that assumption as well.

Thanks

EDIT#1:

I must say that for the tonal case,  I think it might be straight forward as many of you have indicated, vis-a-vis estimating the tonal parameters, reconstructing, and subtracting. The actual problem though, is a little more complicated...

The unwanted signals can range from anything like a group of whales having a 'conversation' to a boat's motor overhead. In such cases no known robust model exists, all I know is that my signal of interest is among all those other energies - they are seperable in frequency, but I need a way to remove them as well, without filtering because that removes my sharpness.

FYI, I think this falls into the realm of non-linear filtering or other techniques from different areas? I should also add that I have a lot of computing power.

A sharp edge or transient may require a very wide-band spectrum to represent it.  Thus any classical or frequency domain filter that completely removes a spectral band could unsharpen or "smear" a sharp transient, by removing a necessary component of that transient.

There are two non-linear methods based on analysis/resynthesis techniques that might work.  Either estimate your tone or interference signal precisely enough to synthesize it and subtract it.  Or characterize your desired signal well enough to resynthesize it without the noise.  You will need a very accurate model for either of these.  Both of these will only remove or add the appropriate fraction of energy in any spectral band, instead off all or none, or some random percentage (relative to the actual S/N) given by a filter response curve.


Why is my filter unstable and self-oscillating in this case? I am using resonant bandpasses to simulate modes of guitar/piano strings for audio synthesis. While attempting to introduce coupling effects I have encountered a problem.

I will illustrate with some simple code:

exciterFinal = exciter - (output * scalar_0to1);

resoBP.setInput(exciterFinal);
resoBP.setFilterCutoff(frequency);
resoBP.setFilterQ(Q);

output = resoBP.getOutput();

In other words, I am subtracting the prior sample's output multiplied by a scalar from 0 to 1 from the usual input of the bandpass to simulate a damping effect.

The most extreme effect would be with a scalar of 1 which should immediately damp the bandpass completely.

This works fine at frequencies below around ~9.4 kHz with normal damping behavior at 48 kHz sample rate. But once you go above that frequency with a scalar of 1, the filters tend to immediately self oscillate to overload. I tried with a State Variable Filter and a Butterworth Filter and had the same result with both.

Increasing the sample rate to 96 kHz allows it to work up to about 17.7 kHz. Above that point it still oscillates. So sampling seems to be the limiting factor. But I cannot afford the processing power to go up this high or higher to get the full spectrum of sound.

What is the reason for this behavior and why does it only happen at frequencies above roughly 19% of the sample rate? Is there any way to fix it?

Update: It just seems the filters become unstable when they are too aggressively damped this way at high frequency and low sample rate. My solution is to damp the filters instead by modifying their Q values when I need aggressive damping and be more gentle when I need to use the subtraction-from-input method. Good enough I guess.

The question is not very specific so this will also be just a general answer.

The stability of a composite linear-time-invariant (LTI) system composed of smaller LTI systems cannot be deduced from the stability of the component systems if the composition introduces feedback. Instead, you should test the stability of the full composite system. The system is stable if the poles of the transfer function lie inside the $z$-plane unit circle. So write out the transfer function of the component systems and also the transfer function of the composite system in terms of the component transfer functions. Find the poles of the transfer function of the composite system. Write an inequality that is the stability condition, and you may be able to solve from that an inequality for the acceptable values of a system parameter of interest, in terms other parameters.


Bandpass filter for audio wav file How do I apply a bandpass filter on an audio file which is in wav format?
And write back the filtered data in a new wav file for further analysis.

Consider looking at this StackOverflow answer which provides the full code for creating a Butterworth bandpass filter. In the sample code the answerer provides, the filter is applied to a manually constructed simplistic signal x. If you substitute this signal with your own NumPy array, it should achieve the intended effect.

WAV files can be quickly made into NumPy arrays, consider using the librosa module. librosa.core.load will be what you use. To write out the filtered signal, simply use librosa.output.write_wav.


What is Imaginary in Fourier transform? How to plot graph of $e^{-t}$ in frequency domain. What would be the axis? If its Fourier transform is $1 /(1+j\omega)$, then how can we plot imaginary on frequency domain (amplitude vs frequency graph) .

First of all, in order to have a Fourier transform, the original signal has to be multiplied by the unit step function:

$$x(t)=e^{-t}u(t)$$

Giving indeed the transform:

$$X(\omega)=\frac{1}{j\omega+1}$$

The way to plot a complex function on the frequency domain is by finding both its amplitude and its phase, and drawing one graph for each.

For the amplitude case, find $|X(\omega)|$, which is:

$$X(\omega)=\frac{1}{\sqrt{\omega^2+1}}$$


Time domain to frequency domain conversion of audio signals to extract 1/3 octave frequencies I am writing first time in this forum and I am not expert in programming and FFT. We have developed an android app (Noise Tracker) for noise measurement using smartphones. It displays noise levels in Leq sound pressure levels.I want to implement A & C weighting. As per my knowledge it is desired to convert received signals from time domain to frequency domain and extract the noise levels for each 1/3 octave frequency band.it looks
Can anyone suggest me how to proceed or there is a ready available program that can be used as such.
your help can solve the problem

Thanks in Advance

As per my knowledge it is desired to convert received signals from time domain to frequency domain and extract the noise levels for each 1/3 octave frequency band.

No. This is needlessly complicated. A-weighting and C-weighting can be implemented through simple FIR or IIR filtering directly in the time domain. There are public domain filter designs avialble. For example check out https://www.mathworks.com/matlabcentral/fileexchange/69-octave

  It displays noise levels in Leq sound pressure levels.

Your biggest problem will be to deal with microphone calibration. Both frequency response and overall sensitivity of Android microphones are all over the place and vary a great deal from model to model, can change between Android version and sometimes are actually actively modified by each carrier. Without some sort of calibration, you can display a number but it won't be Leq and it will mostly be meaningless.

iOS might be easier: Apple does a much better job in controlling microphone consistency in design and supply chain. For a meaningful Leq estimate you still need calibration, but it's going to be less severe and there is a much smaller number of models/iOS versions to deal with.


Simulating spatial perception - with cosine wave based sound - in MATLAB I'm trying to create cosine wave, one for the right channel and another for the left. The one on the left is to be moved in phase to simulate change to the spatial perception of the tone. A delay in time to the left side makes us think the sound comes from the right, and since time is converted to phase in a cosine I am working with the phase.
I wish to use a $500\textrm{ Hz}$ cosine and a $2\textrm{ mSec}$ signal. MT ultimate goal here is to have 9 sounds to simulate the locations around a person's head.  1-4 being left, 5 the center and 6-9 the right.
This a snippet of the code responsible for this:

t = [ 0 : 1 : 2000 ];   % Time Samples
f = 500;    % Input Signal Frequency 500Hz
fs = 44100; % Sampling Frequency

x =  cos(2 * pi * f / fs * t);  % Generate Sine Wave

y =  cos(2 * pi * f / fs * t - DELAY);  % Generate Sine Wave

Y = [x,y];

audiowrite(['file' num2str(j) '.wav'], Y, fs);

The issue I'm having is that no matter what values I try for DELAY the resulting file sounds just like any other, I cannot hear any change.
What am I doing wrong?

Try this

DELAY = 0.5e-3; % half a millisecond for starters.
t = [ 0 : 1 : 2000 ]';   % Time Samples, COLUMN vector
f = 500;    % Input Signal Frequency 500Hz
fs = 44100; % Sampling Frequency
x =  cos(2 * pi * f / fs * t);  % Generate Sine Wave
y =  cos(2 * pi * f * ( t/fs - DELAY));  % Generate Sine Wave
Y = [x,y];

A bunch of problems with your code

Your variable t was a row vector so you generated a mono output signal
The units of your DELAY variable are unclear. I have written it in a way so that it is in actual seconds.


What is the difference between time complexity and sample complexity? What is the the difference between time complexity and sample complexity in the context of algorithm development?

I'm currently studying the sparse Fourier transform and this idea of sample complexity keeps coming up. The majority of sample complexity literature is in the context of machine learning which is not what I'm interested in.

Some of those notions can seem quite vague, and may even depend on the field.

In algorithmics, different types of complexities come at play. One could be related to actual implementations (machine or computational complexity). Another to data access (sample, training set complexity).

Computational complexity is generally machine-, transmission-, power dependent. You have for instance time-complexity (the time need to run) or space-complexity (the memory required).

Data complexity are generally less machine dependent, and more conceptual. Sample complexity can be an estimate of the number of operations required, depending on the number of samples. In the Landau notation, an $O(N^3)$ algorithm grows as the cube of the number of samples $N$. It is a gross measure of scalabilty. Here, one is not concerned with data precision, storage, etc. If one can bound the execution time of operations, one could relate sample complexity to time complexity, but complicated computer architectures make this conversion doubtful: caching, pipelining, parallelizing,  swapping alter a lot the expectations.
In learning tasks, "data complexity" may refer to the quantity of data to be learned, on the training sets, to achieve certain prediction goals. One classical measure is the  Kolmogorov complexity.

A related SE question: Computational complexity vs other complexities.


Does the auto-correlation function of stationary random process always converge? The auto-correlation function of the stationary random process only depends on the time difference $\tau$.

http://web.ntpu.edu.tw/~yshan/chapter6_han.pdf

64th slide of this lecture note mentions that for a zero-mean process, the autocorrelation converges to zero as $\tau$ goes to infinity.

I want to know why.

No it does not necessarily.

For example the following discrete-time, WSS random process

$$x[n] = A \sin(\omega_0 n + \phi) $$

which is called the random phase sinusoid, where $A$ and $\omega_0 \neq 0$ are fixed values and $\phi$ is a random variable uniformly distributed in $\phi \in [-\pi,\pi)$ has an auto-correlation function of the form

$$ r_{xx}[k] = \frac{A^2}{2} \cos(\omega_0 k) $$

which does not go to zero as $k$ goes to infinity; $ \lim_{k \to \infty} r_{xx}[k] \neq 0$.

Similarly for a continuous-time process, the same can be shown.

Note, however, that as MattL indicated in his answer as well, the information contained within a random process is mostly included in its innovations part, this is also expressed in Wold decomposition theorem that any random process can be broken into two parts as a predictable periodic part and a regular unpredictable part (which is the innovations part), then if a WSS random process only includes a regular part but no predictable part, then its covariance (or correlation if zero mean) sequence shall go to zero as the lag goes to infinity if it has a uniformly convergent power spectral density; i.e., the DTFT of auto-correlation sequence (ACS) should converge, and this requires the sum of abs value of ACS be finite and this also requires that $\lim_{k \to \infty} r_{xx}[k] = 0$ for a zero-mean, WSS process containing no predictable (periodic) parts in it.

Furthermore, when the concept of ergodicity is introduced, then one of the necessary conditions for a WSS random process to be ergodic in the mean is that its auto-covariance (or equivalently auto-correlation for a zero mean process) should go to zero as $k$ goes to infinity. May be that was stated in that document.


Split audio signal into frames, apply window function and merge again First off: I'm really new in the signal processing department and MatLab as well, so please bear with me here.

What I want to achieve is to split my audio signal into smaller 25ms frames, where the frames overlap with 50%, apply a window function, thus the overlapping (probably Hamming or Hann), and then add it together again (write out a new sound file), so that I basically get the original signal back. Someone mentioned overlap-add, but I'm not quite sure if convolution is correct here, as mentioned, I'm new to all this.

So I have a sound file that I read like this:

function output = readSoundFile(filename, channels)

if isscalar(channels)
    chN=channels;
    channels=ones(1,chN);
else
    chN=length(channels);
end

fid=fopen(filename,'r','n');
[output, elementCount] = fread(fid,[chN,inf], '*int16');
fclose(fid);

for i=chN:-1:1
    if channels(i) == 0
        output(i,:)=[];
    end
end

end

In another MatLab file, I use the readSoundFile function:

fileUrl = 'whatever/file.format';
sound = readSoundFile(file);
soundsc(double(sound), 16000); % sound is sampled with 16kHz. Don't know, if this information is needed

How do I split the resulting data into 25ms frames with an overlap of 50%, apply a window function, and write out a new file?

Should it help, I can convert the file into a wav file, if wavread/audioread would be any easier.

Von Hann (but not Hamming) windows with a 50% overlap sum to unity (except possibly at the very start or end of a file).  If you do any processing, rather than just analysis, in the frequency domain, you may find overlap add or overlap save with a rectangular window to be more appropriate.


What is spectral whitening? What is meant by "spectral whitening" in DSP?

What effect does spectral whitening have when used in image processing?  (visually or otherwise...)

Where might spectral whitening be useful in audio processing or analysis?  What would a spectrally whitened audio signal sound like?

What is meant by "spectral whitening" in DSP?

Spectral whitening is usually an attempt to make the spectrum of the signal "more uniform".  One reason this might be a good thing to do is that it can have the effect of making the autocorrelation of the signal "narrower" (and closer to a Dirac delta, for discrete-time signals). This can help localize in time.

  What effect does spectral whitening have when used in image processing? (visually or otherwise...)

It's generally not pretty.  Most images are "low pass" (most of the information is in the low frequency part of the spectrum).  One simplistic approach to whitening in images is to do a column-wise (or row-wise) difference (i.e. diff in matlab).

This will mean negative pixel values, which generally do not map to anything sensible with standard images.

This example shows how prewhitening can improve localization in image processing template matching. The picture from that link is:

  Where might spectral whitening be useful in audio processing or analysis?

If you are trying to localize (in time) the onset of a sound, then it's possible that spectral whitening can improve this.  It's also possible that it can reduce (disimprove) the SNR.

  What would a spectrally whitened audio signal sound like?

For audio of speech or music, it will tend to bring in more higher frequencies.


Time signal FFT different for whole signal than half signal on high frequencies I have a time signal of u from 0 to 1000 time units. If I chop the signal in two halves and compare the frequency spectrum of both it looks similar, as it should (the signal should be statistically steady from where I start recording it). However if I compare it to the frequency spectrum of the full signal I get a different slope in the high-frequency region.

I do not understand this. If I use the whole signal there should be a change on the low-frequency region rather than the high-frequency one right?

This is the code I use compute the frequency spectrum. I think it's okay, but maybe I'm doing something wrong.

def time_spectra(t, u):
    import numpy as np
    from scipy.interpolate import interp1d
    # Re-sample u on a evenly spaced time series (constant dt)
    u_function = interp1d(t, u, kind='cubic')
    t_min, t_max = np.min(t), np.max(t)
    dt = (t_max-t_min)/len(t)
    t_regular = np.arange(t_min, t_max, dt)[:-1] # Skip last one because can be problematic if > than actual t_max
    u_regular = u_function(t_regular)

    # Compute fft and associated frequencies
    uk = np.abs(np.fft.fft(u_regular)) / len(u_regular)
    freqs = np.fft.fftfreq(u_regular.size, d=dt)

    # Take only positive frequencies
    freqs = freqs[freqs > 0]
    uk = uk[:len(freqs)]

    return freqs, uk

The primary issue is due to using a rectangular window with the floor being the result of spectral leakage.

This is demonstrated in the plot below showing the frequency response (kernel) for the two rectangular windows used; one 500 samples long as red and the other 1000 samples long as blue. The corresponding rectangular window is effectively multiplied by your signal in time, thus convolves with your signal in frequency, resulting in spectral leakage as the spectrum given in the original question.  Thus you are seeing the smearing (convolution) of the actual spectrum with this kernal in each case. Using a window with a wider dynamic range (lower sidelobes) such as the Blackman  would result in a better view of the actual spectrum to the extent that spectral components now are getting buried under the spectral leakage.

For further details on windowing and window choices, please see this post: Why would one use a Hann or Bartlett window?

You will also see a 3 dB difference between two windows if the signal energy is distributed wider than the resolution bandwidth in each case, since the resolution bandwidth for the rectangular window is $1/T$ where T is the length of your waveform. (So when you double the length from 500 to 1000, the resolution bandwidth is half resulting in half the power in each bin (unless a spectral tone occupies one bin in which case you would not see a change).


How to find the output signal of a filter using state space matrices? I have a filter.
It has two poles and two zeros.

I found the state space equations and the matrices A, B, C, and D

Now. I have 9 samples that I need to process with my filter.
How do I use A,B,C,D matrices to process the samples?

We don't know anything about your system, but general rule is:

$$\begin{array}
&y[n] &= \mathbf{C}x[n] + \mathbf{D}u[n]\\
x[n+1] &=\mathbf{A}x[n] + \mathbf{B}u[n] \\
\end{array}$$

Where $u[n]$ is the input signal value, $x[n]$ is state vector value at a given point in time, and $y[n]$ is obviously the output.

What you need to do is to:

Initialize the initial state $x[0]$ to some arbitrary value (usually vector of zeros).
Run the calculation by consequently calculating values of $y[n]$ and $x[n+1]$ to be used in the next step.

Speaking Matlab'ish (not tested, cause I am in travel, but given all informations its extremely easy to figure out what is going on):

b = [0 2 3];
a = [1 0.4 1];
[A,B,C,D] = tf2ss(b,a);

x_0 = zeros(size(A,1), 1); % Initial state
N = 9;
u = randn(1, N); % Some input signal
y = zeros(N, 1); % Output signal

x = u_0;
for n=1:N-1
  y(n) = C*x + D*u(n); % System output
  x = A*x + B*u(n);    % Transition
end


What's the state-of-the-art algorithms in image edge detection? There are several state-of-the-art algorithms in image denoising including Non-Local Means, BM3D and so on.

I am wondering whether there is some well-known state-of-the-art algorithms in image edge detection (with reproducible code) ?

Thanks a lot for your answering.

Despite its age, Canny Edge Detection is still a state of the art filter. The results produced by this algorithm make for it always being included in image editing software. Solid and descriptive edges that are often overly represented by other filters. It lacks the simplicity of, say, the Boolean Edge Detection, included in the paper, "Edge Detection Techniques: Evaluations and Comparisons" by Ehsan Nadernejad, Sara Sharifzadeh, Hamid Hassanpour.

The quality of Edge Detection and Segmentation filters are difficult to appraise, and algorithm speed is usually the most important consideration when implementing them.

For code and discussion, you should take your research to codeproject.com, and lookup Image Edge Detection.


Question about zeroes of simple 2nd order FIR filter I have a question about the zeroes in a simple second order FIR filter.  I'll list out the major beat points of calculating where the zeroes are before asking the question (I'm following "Designing Audio Effect Plug-Ins in C++", section 5.17).

Starting with a simple second order FIR filter difference function:

$y(n)=a_0x(n)+a_1x(n-1)+a_2x(n-2)$

We can z transform it to:

$Y(z) = X(z)(a_0+a_1z^{-1}+a_2z^{-2})$

Which gives us the transfer function:

$H(z) = a_0+a_1z^{-1}+a_2z^{-2}$

We can factor out $a_0$ to make it a gain parameter, and replace $a_1$ with $\alpha_1$, and $a_2$ with $\alpha_2$, where $\alpha_1 = \frac{a_1}{a_0}$ and $\alpha_2 = \frac{a_2}{a_0}$:

$H(z) = a_0(1+\alpha_1z^{-1}+\alpha_2z^{-2})$

To find the zeros, we can ignore $a_0$ and just look at the part inside the parentheses.

$0=1+\alpha_1z^{-1}+\alpha_2z^{-2}$

We know that the zeroes need to be complex conjugates, so can change it to this:

$0=(1-Z_1z^{-1})(1-Z_2z^{-1})$

where $Z_1=Re^{j\theta}$ and $Z_2=Re^{-j\theta}$

Using Euler's equation, we can expand and do algebra to get this:

$0=1-2R\cos(\theta)z^{-1}+R^2z^{-2}$

Comparing that to the earlier equation:

$0=1+\alpha_1z^{-1}+\alpha_2z^{-2}$

We can see that this equation has a zero where:

$R=\sqrt{\alpha_2}$

$\theta=cos^{-1}(\frac{\alpha_1}{-2R})$

The other zero has the same radius $R$, but has a negative theta ($-\theta$).

Ok so my question is this. Using values $a_0=1$, $a_1=-1.27$, $a_2=0.81$, we get $R=0.9$ and $\theta=\frac{\pi}{4}$.

That seems to be correct and works well, but let's say that $a_1=-3.6$. We get the same radius $R$, but to calculate $\theta$ we have:

$\theta=cos^{-1}(2.0)$

... which is not a valid thing.  cosine never has a value of 2, so acos(2.0) is not defined.

Where did i go wrong or how is this situation handled?

The reason for this is that the maths centres on a false assumption:

  We know that the zeroes need to be complex conjugates

Do the zeros need to be complex conjugates? Can they not be 2 real zeros?

I sympathise here because I also read this book and the maths portion ended up utterly confusing me because of very odd methods.

A more general method would be at this point:

$0=1+\alpha_1z^{-1}+\alpha_2z^{-2}$

Rearrange thus:

$0=\frac{z^{2}+\alpha_1z^{1}+\alpha_2}{z^{2}}$

Here on the top we have a quadratic equation, which you can use any quadratic solver for.

For your first case, I arrive at 0.635 $\pm$ 0.637i, which is a pair of complex numbers denoting the conjugate zeros - correct.

For your second, I arrive at two real zeros, one at 3.36 and one at 0.24 - so it is most definitely able to be computed!

And finally, you will have noted the $z^{2}$ term lying on the bottom, which indicates that there is a square pole at 0+0i. If these zeros were poles then there would similarly be a square zero at 0+0i.

All of these results can be verified with the zplane function in MATLAB.

It took me some time with this book and checking results with the zplane function in matlab to realise that his methods of working really don't work on a general basis, and frankly overcomplicate things. It has however been awhile since I've looked at it so maybe I similarly misinterpreted his maths


How to calculate step response for $y''(t) - y(t) = x'(t) - x(t)$ in time domain? How to calculate step response for $y''(t) - y(t) = x'(t) - x(t)$ in time domain? So, without Laplace or Fourier Transforms.

This is what I tried:

The Homogeneous solution of the differential equation would be: $ y_\text{homogeneous} = A*e^t + B*e^{-t}  $.

$u(t)$ is the unit step function.

$$ u(t) \triangleq
\begin{cases}
1, & \text{if }t>0 \\
0, & \text{if }t<0
\end{cases} $$
When I try to fill in the stepfunction for $x: x(t)=u(t)$. So $x'(t) = u'(t) = \delta(t)$

And I try a step response $y(t) = y_\text{homogeneous}(T)*u(t) =  (A*e^t + B*e^{-t})*u(t) $

$$\frac{\partial ^2 (  \left(A e^t+B e^{-t} \right)u(t) )}{\partial t^2}-\left(A e^t+B e^{-t}\right) u(t)=\delta(t)-u(t)$$

So, I will get:

$2 (A-B) \delta(t)-A \delta(t)+A \left( \delta ' (t) \right )+B \delta(t)+B \left( \delta ' (t) \right)=\delta(t)-u(t)$

But there's no $u(t)$ at the left hand side.

It is indeed possible to convert the system to an equivalent first order system, but let's first directly solve the given problem. The system's step response satisfies

$$y''(t)-y(t)=\delta(t)-u(t)\tag{1}$$

This implies that we can compute the step response of the given system by considering another system

$$y''(t)-y(t)=x(t)\tag{2}$$

and compute its impulse response ($x(t)=\delta(t)$) and its step response ($x(t)=u(t)$), and obtain the step response $(1)$ as the difference of the impulse response and the step response of the new system $(2)$.

For the step response we need to solve

$$y''(t)-y(t)=u(t)\tag{3}$$

Assuming the system is initially at rest, we get initial conditions $y(0)=y'(0)=0$. The homogeneous solution of $(3)$ is $y_h(t)=Ae^{t}+Be^{-t}$, $t>0$, whereas a particular solution is $y_p(t)=-1$, $t>0$. The complete solution is then

$$y(t)=y_h(t)+y_p(t)=Ae^{t}+Be^{-t}-1,\quad t>0\tag{4}$$

The constants $A$ and $B$ are determined by the initial conditions:

$$\left.\begin{align}y(0)&=A+B-1=0\\y'(0)&=A-B=0\end{align}\right\}\Rightarrow A=B=\frac12\tag{5}$$

Consequently, the step response of system $(2)$ is

$$y(t)=\frac12\left(e^t+e^{-t}\right)-1,\quad t>0\tag{6}$$

Its impulse response satisfies

$$y''(t)-y(t)=\delta(t)\tag{7}$$

and it can be obtained by differentiating the step response $(6)$:

$$y(t)=\frac12\left(e^t-e^{-t}\right),\quad t>0\tag{8}$$

Finally, according to $(1)$ and $(2)$, the step response of the original system is the difference between the impulse response and the step response given in $(8)$ and $(6)$:

$$y(t)=\frac12\left(e^t-e^{-t}\right)-\frac12\left(e^t+e^{-t}\right)+1=1-e^{-t},\quad t>0\tag{9}$$

Another approach would be to see that the original system can also be represented by a first-order system. Note that

$$y''(t)-y(t)=(y(t)+y'(t))'-(y(t)+y'(t))\tag{10}$$

So if $y''(t)-y(t)=x'(t)-x(t)$, it follows from $(10)$ that

$$y(t)+y'(t)=x(t)\tag{11}$$

With $x(t)=u(t)$ you obtain from $(11)$ the same step response as in $(9)$:

$$\begin{align}y_h(t)&=Ce^{-t}, &t>0\\
y_p(t)&=1, &t>0
\end{align}\\y(t)=y_h(t)+y_p(t)=Ce^{-t}+1,\quad t>0$$

From $y(0)=0$ we get $C=-1$, and the step response becomes

$$y(t)=1-e^{-t},\quad t>0$$

just like in $(9)$.

As a final check it is useful to see what happens when you solve the problem using the Laplace transform. The original difference equation transforms to

$$Y(s)(s^2-1)=X(s)(s-1)$$

which gives for the transfer function

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

This shows that you have a pole-zero cancellation, and the system is actually a first-order system, just as described by Eq. $(11)$. For $x(t)=u(t)$ you get $X(s)=1/s$, and, consequently, for the Laplace transform of the step response

$$Y(s)=\frac{1}{s}\frac{1}{s+1}=\frac{1}{s}-\frac{1}{s+1}$$

which corresponds to the time domain function

$$y(t)=u(t)-e^{-t}u(t)=(1-e^{-t})u(t)$$

which is the same as $(9)$.


A Fast Power Grid Frequency Estimation Approach Using Frequency Shift Filtering I have read the paper A Fast Power Grid Frequency Estimation Approach Using Frequency Shift Filtering.
I want to prove it using a numerical example, I have problem how to calculate the hp(n) values,
Let's suppose that I have a sample signal which I have sampled with an ADC and for simplicity I will use Matlab to make the signal. (the signal has a 50.1HZ frequency) and I want to estimate it with this paper approach,
fs = 12800;                    % Sampling frequency (samples per second)
dt = 1/fs;                   % seconds per sample
StopTime = 0.02;             % 20 mili seconds or a 1 full cycle of 50Hz AC
t = (0:dt:StopTime-dt)';     % seconds
F = 50.1;                    % Sine wave frequency (hertz)
X = sin(2*pi*F*t);           % Generate sine wave

Now How should I calculate the hp(n) values using MATLAB?
Finally I should convolve the result with this sample sin wave to get 50.1 Hz estimation.

Basically they cascade moving-average filters, or CAF.
a Moving-average filter has the following coefficients
$h_{av}(m) = 1/M, m = 0, 1.. M-1$
In Matlab it would look like this
h = ones(1,M)./M;
y = filter(h,1,x);

If you want to cascade three moving-average filters
h = ones(1,M)./M;
y1 = filter(h,1,x);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);

M should be equal to the number of samples per period. So if your frequency is approximately 50 Hz and your sampling frequency is 12800 Hz then M = 256;.
Edit : I modified your script to make it work. There were 2 problems, first of all you only processed 256 samples, your couldn't measure anything significant because the transient response of your 3 cascaded filters is 768 samples long. Secondly, you filtered the reference and not your complex signal...
clear
clc
%ADC sin wave, it represent the data that I have captured with the ADC
%and want to estimate the frequency
ADCsin = dsp.SineWave(1,50.1);
ADCsin.ComplexOutput = false;
ADCsin.SampleRate  = 12800;
ADCsin.SamplesPerFrame = 12800;
ADC = ADCsin();
plot(ADC);
grid on
hold on

%ref sin wave with 50Hz complex value
Refsin = dsp.SineWave(1,50);
Refsin.ComplexOutput = true;
Refsin.SampleRate  = 12800;
Refsin.SamplesPerFrame = 12800;
ref = Refsin();
plot(imag(ref))

%calculating  xs
xs = ref .* ADC;

%calculate the CAF
M = 256;
h = ones(1,M)./M;
y1 = filter(h,1,xs);
y2 = filter(h,1,y1);
y3 = filter(h,1,y2);

%calculating the final result
Fest = 50 - (angle(y3(end)) - angle(y3(end-1))) * 12800/(2*pi)