FFT "Brick Wall" Low Pass Filter Not Smoothing Output I'll try to be as brief as possible. I have very limited signal processing experience. I have a data set that consists of essentially a big old sine function (with deviations from ideality), and lots of high-frequency noise (what I don't want). I desiced on a low pass filter because it doesn't change amplitude much or significantly shift peak/band locations.

I played around with QTIplot (scientific plotting soft-ware) and found I can smooth the dataset quite well by simply doing a FFT Low pass filter with a  1Hz cut-off (seem image below).

So I looked at the source code in QTI plot to get an idea of how they implemented it. It turns out it's a brick wall algorithm. Their code goes like this -

//(C) 2007 - 2008 by Ion Vasilief, Knut Franke

double df = 1.0/(double)(d_n*(x[1]-x[0]));//Sampling frequency
gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(d_n);
gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(d_n);
gsl_fft_real_transform (y, 1, d_n, real, work);//FFT forward
gsl_fft_real_wavetable_free (real);
...
//Brick wall
for (int i = 0; i < d_n && ((i+1)/2)*df < d_low_freq; i++)
    y[i] = 0;
break;
gsl_fft_halfcomplex_wavetable *hc = gsl_fft_halfcomplex_wavetable_alloc(d_n);
gsl_fft_halfcomplex_inverse (y, 1, d_n, hc, work);
gsl_fft_halfcomplex_wavetable_free (hc);
gsl_fft_real_workspace_free (work);

So I borrowed from their code but translated it into JAVA and decided on using the JTransforms FFT library. The JTransforms FFT library leaves the first two index positions for storage of the transform and the rest of the indices are alternating Re, Im components ("Bins").

    int size = recievedData.size();
    DoubleFFT_1D fftDo = new DoubleFFT_1D(recievedData.size());
    double fft[] = new double[size*2];
    fftDo.realForwardFull(fft);
    double df = 1.0/(((double)size)*(0.027413951));
    for(int i = 2; i < 2*size && (((double)i+1.00)/2.00)*df < cutOffFreq; i++){
        fft[i] = 0.0;
    }
    fftDo.complexInverse(fft, false);

Now when I run my code I see no real noise reduction. Sometimes I see the minima of my signals increase to some weird baseline, if I allow the FFTinverse to scale my code sometimes I see all blocks. Ultimately it seems the code does little to nothing to my signal, while the QTIPlot code is exactly what I want.

So a couple questions:
- What am I doing wrong in my code? Or is it the library differences (GSL vs JTransforms)?
- Is there a better way to do this kind of filtering? If so can you take me to a resource to learn how it's implemented. I read about using masks everywhere but could not find a single example of them that wasn't "hard-coded", IE: my data may be dynamic. Seems most people who did masks like to show pictures but not detail how to implement them even in lecture powerpoints from universities.

Also this is not for homework or anything I'm just playing around with an electronics project.

Sorry for the post I just solved it.

The problem was I was using a full transform. The half-transform utilizes a symmetry condition and for some reason works infinitely better!


Methods in Singing Synthesis I want to know more about the latest attempts to create synthesizers for singing. I would ideally like to know more about both the most widely respected algorithms/models as well as mainstream software or programs that implement them (whether free or proprietary). Can someone offer me resources that would give me a decent grasp of the above?

I chose to post this question on Signal Processing Stack Exchange since I know such software relies to some degree on signal processing type methods for analyzing or creating waves.

I'm very new to this area of research so any suggestions for where to start would be welcomed.

Thanks.

Well, there is sampling of singing voices.  it's not the latest.  one thing that you perhaps need to worry about is that, unless you get a sample root for each note going up the MIDI note scale, then you will likely be pitch shifting the notes up or down a little.  it turns out that this causes a noticeable distortion sometimes called "munchkinization" or the "Alvin the Chipmunk" effect.  to avoid that, each cycle of the quasi-periodic note needs to be isolated into little "wavelets" (not to be confused by the wavelet transform) or "grains" or "FOFs".  then you would output the grains or FOFs at a rate of the fundamental frequency of the note without stretching or scrunching the grain (this requires overlapping and adding).  there's an old paper by Keith Lent that discusses pitch shifting vocal notes without shifting formants, and another paper by me that analyzes the algorithm a little more rigorously.  the latter, i can send a copy.

but for synthesis, you need not perform the real-time pitch detection that you would have to for pitch-shifting the voice live.  what you can do is sample some notes from these sung voices and perform analysis on them off-line and generate a file of these grains or FOFs and store them in your synthesizer.


fast 2D cross correlation algorithm I want to match a small template to a larger image, that the distance between the template and the subimage with the same size of the template is minimized. It can be solved directly or by applying 2D cross correlation, and both methods have an O(n^4) time complexity. Are there any method to simplify the algorithm by either reducing the complexity or applying pre-processing?

If your template or kernel is small, then straight convolution might be the fastest approach. There's a crossover point when performing convolution in the frequency domain is faster than straight time/spatial domain convolution and it can be hardware dependent, but usually when the kernel (template) approaches 1/4-1/2 the size of the image frequency domain convolution is faster.

If you have a multicore machine you can split the larger image into say 4 quadrants and run the spatial convolution on 4 threads and fuse the results together. That should actually be pretty speedy.


Detection of wobbling sound with rising notes In an audio file I would like to find calls of a specific bird specie, one that makes wobbling sound with rising notes. Can you recommend what would be the best audio feature to calculate in time or frequency domain to capture its essence?

Among others I tried MFCC, but the results were not promising (or perhaps my implementation poor) .

[Edit]

In this post I explained briefly the problem and included exemplary audio file with spectrogram: Example.

Here is another example. All marked areas contain kiwi (famous flightless bird from New Zealand) calls I am trying to identify. Mark that the original audio is very noisy; presented spectrogram is a result of applying high-pass filter and spectral subtraction on areas identified as noise-only.

I would start with following set of parameters:

MFCC's (I that know you tried it, but stay with me) without static energy (1'st coefficient)
Some descriptors from MPEG-7, like: Spectral Flatness, Harmonicity, Fundamental frequency, Spectral Spread, etc. You can find out more about them here: click or in this great book: MPEG-7 Audio and Beyond

You should also include $\Delta$'s and $\Delta\Delta$'s of MFCC's and Spectral Flatness. They should really improve accuracy. You can also try to calculate these for other parameters. Having that, you can start with Gaussian Mixture Model or Support Vector Machine (I think that GMM will be suitable for that task). After training your model simply validate it on testing set (seems obvious but I've seen many projects where people tested their Machine Learning algorithms on training set with added noise or, what's even worse, same training data).

After this, you should perform Feature Optimisation. This is important part, as some of the features are irrelevant (i.e. Spectral Flatness for few first filter banks, while bird's chirp has higher frequencies). Nice publications with many algorithms for that purpose can be found here: click,   click.
Good luck!


Time resolution of the Short Time Fourier Transform (STFT) I'm struggling to figure out how the time points of an STFT are calculated, and I can't find a definitive answer.  Let's say I have a 4Hz stationary signal and I'm going to use a 64 second window with 3 second overlap. So that's a 256 point window and a 12 point overlap.

Assuming I start at time=0, take the first 64 seconds, and perform the FFT/Power Spectrum Density/etc... Can I then say that is the value at t=32? Is the next window, after the 3 second slide localized at t=35, and so on?

If so, and I really wanted to start at t=0, would I then effectively start at t=-32, fill the first 128 points with zeros and take the first 128 points from my signal, thus centering on t=0?

There is no single "time instant" associated with a short-time Fourier transform. As you noted, if you perform a DFT on data collected from $t = 0$ to $t = 64$, then there isn't a single point in time that you can associate with the output from that DFT; it is a function of every sample in its time interval.

For this reason, there is no standard convention for how you might denote the time axis associated with successive STFTs, if you were tiling them together into a spectrogram, for example. You'll just need to pick a convention that is meaningful for your application. You've already identified a couple reasonable ones (using the beginning or center of the DFT's window as its "time instant").


On what it means for a system to be minimum phase From my undestanding I get that when a system's poles and zeros are all inside the unit circle then it means that it's minimum phase. But I don't get how the location of zeros and poles are related to the phase being "minimum" or whatever the "minimum phase" notation implies. So basically I'm asking what does it mean for a system to be minimum phase othen than the zeros and poles part and why was the minimum phase as a name was chosen for this kind of systems.

The minimum phase system will have the minimum group delay for a given magnitude response. The phase response will have the least excursion over the frequency domain (due to all zeros being inside the unit circle). Since group delay is the negative derivative of phase with respect to frequency, a minimum phase system also has the minimum group delay.
The minimum phase excursion for a minimum phase system is well demonstrated with the graphic below, knowing that the complex frequency response (magnitude and phase) of an FIR filter can be determined graphically through phasors rotating at a rate according to the coefficient number (the first coefficient rotates zero times, the second coefficient rotates once, etc as the frequency axis is swept from $0$ to $2\pi$). Both filters shown below will have an identical magnitude response, but the one with the larger coefficient appearing first is minimum phase while the one with the coefficients reversed (the reverse filter) is maximum phase. We can see with the minimum phase system how the total phase deviation is minimized, while in the maximum phase case the phase deviation will be maximized, but in both cases having the same magnitude.
The transfer function of the first filter is $1-0.5z^{-1}$ or $(z-0.5)/z$ with all zeroes inside the unit circle as one zero at $z=0.5$ (minimum phase!). While the transfer function of the second filter is $0.5-z^{-1}$ or $(z-2)/z$ with all zeroes outside the unit circle with one zero at $z=2$ (maximum phase!). Notice too how in the minimum phase case, with the dominant tap placed earlier, how most of the signal will be at the output of the filter sooner (least delay).
When all the zeros are inside the unit circle, the frequency response as equivalently mapped out as in these graphics by rotating phasors will never enclose the origin, thus the phase excursion never exceeds $2\pi$.

For further details and how this relates to Cauchy's Argument Principle please also see: Can I set a constraint on the first tap of an FIR filter such that its inverse is stable?


Reduce noise of sensor by combining two I've had the basics of signal in university but little real-life practice. I've got an acceleration sensor with a digital output signal (signed 16 bit value). This signal is a little noisy. The standard deviation is about 20 (referenced to the full 2¹⁶-1 scale of the sensor).

On the hardware I have, there are two of these sensors present. Is it somehow possible to merge the signals of both sensors, in order to reduce the noise?

The reason I'm asking this is because the noise of the two sensors should be uncorrelated (at least in theory) whereas the acceleration signal is correlated. Is there a way to reduce the noise?

Assuming that the sensors share the same characteristics, have the same timing (acceleration signals are aligned), the model with $y_1= x + n_1$ and $y_2= x + n_2$, $n_1$ and $n_2$ being uncorrelated noises of the same power, averaging them is a way to reduce the noise. The theory that asymptotically, averaging $N$ sensors reduce the variance by a factor of $N$, hence the noise amplitude by $\sqrt{N}$. For small $N=2$, and with quantization, one should not expect too much (since the asymptotic formula is not valid), but this can be tried, you may have nice surprise. If noise powers are different, this method can be improved to take this into account.

Further, this problem can be addressed as a form of source separation, for which many methods have been developed (independent component analysis, Non-negative matrix factorization, low-rank decomposition to name a few). The choice will depend on how much more information you have, and your algorithmic constraints.


print a date and time in command window MATLAB I use matlab R2015a, when I run my script I want print the date and the time of exection

I tried this code

Dtimes=datetime('now','TimeZone','local','Format','d-MMM-y HH:mm:ss Z');
% the first test
disp(' Date and time is ' Dtimes ); % But not work
% second test
disp(sprintf('date and time is %s  ',Dtimes)) %also not work

But with this command is worked :

disp(sprintf('date and time is %s  ',datestr(now)));

What is the problem?

Based on Mathworks' documentation here, datetime returns an array, not a string. Therefore, you cannot use disp to display a string Date and time is together with that returned array.

Your second method is also not working because again Dtimes is an array and you cannot use %s to print it.

Your last method is working because datestr returns a string see documentation.

Please let me know if this explanation is not clear/enough.


How to simulate analog VU meter response to audio? I'd like to simulate the movement of an analog VU meter dial as shown below and often found on stereos and other music monitoring equipment.

I just learned through this article by Rod Elliot (ESP) and Wikipedia that the response of a VU meter is standardized and designed to measure the perceived loudness of audio signals.

  A VU meter is designed to have a relatively slow response. It is driven from a full-wave averaging circuit defined to reach 99% full-scale deflection in 300ms and overshoot not less than 1% and not more than 1.5%. Since a VU meter is optimised for perceived loudness it is not a good indicator of peak (transient) performance. Nominal sensitivity for 0VU is 1.228V RMS, and the impedance is 3.9k.

I'm playing with music visualization and data processing and a measure of perceived loudness is precisely what I need. I also like the analog nature of the response.

I'd like to implement this in an algorithm in Python (for an entire music file) and later perhaps in Processing (live processing in batches). The underlying implementation will probably be very similar.

How can I generate the response from raw audio data to match the specifications of a VU meter?

if you have input signal as vector of floating point numbers, try following code. It won't behave as real VU meter, but by modifying two constants you should get something very similar. Visualize the value of C in logarithmical scale

// memory element (can be imagined as capacitor)
c = 0
kCharge = 0.1
kDischarge = 0.001

for each sample x:
{
  // rectify the signal
  x1 = abs(x)

  if x1 > c then
  {
    // charge
    c = c * (1-kCharge) + x1*kCharge
  } else
  {
    // discharge
    c = c * (1-kDischarge)
  }
}


Inverse discrete Fourier transform If anyone can help solving this exercise I'll be grateful. It's urgent.
(I've added my answer, but I think it's wrong)

No it's not correct.

Derive it like this :

1-) $2N$-point DFT of $x[n]$ is: $X_{2N}[k] = \sum_{n=0}^{2N-1} x[n] e^{-j \frac{2\pi}{2N} k n }$

2-) $Y[k] = X_{2N}[2k+1]$ be the odd-indexed samples of $X_{2N}[k]$

3-) We are looking for $y[n] = \text{N-point IDFT}\{ Y[k]\}$.

4-) Elaborate on step-2 and step-1 to see that $Y[k] = \text{N-point DFT}\{ x[n] e^{-j \frac{\pi}{N}n}\}$

5-) Then from steps 3 and 4 we get :

$$y[n] = \text{N-point IDFT}\{ \text{N-point DFT}\{ x[n] e^{-j \frac{\pi}{N}n}\} \}  $$

$$ y[n] = x[n] e^{-j \frac{\pi}{N}n}  $$


Is there a algorithm for detecting bursts in signal? Suppose we have a signal like that and we want to know when a burst have happened. Rigth now our approach is pretty straightforward. We measure std of signal at some time window and then with each new signal bin we measure it's difference with respect to previous and if it's much more than std, then we have found our burst. But maybe this is known problem and there is a library for that, or the way to do it more efficient.

signal time plot

You want to build a burst detector.

First of all please notice that you can't do magic so you will always have false alarms and miss detections.

We'll suppose that you have a 1-D signal and that what you call a burst is a short peak in signal intensity, we'll also suppose that your burst is very close to 1 point in duration.
If your "burst" is in fact a short deterministic signal whose shape is known, you'll want to first use a matched filter, if your burst is a short random signal whose second order moment is known, you'll want to first use a stochastic matched filter. If you are neither in these cases, you'll need to use more advanced and/or more had-hoc algorithm.

Now we have the 1 point long burst, or signal, we want to detect, it is burried in some random noise.

Usualy people go for ad hoc approach, he main one being the one you stated.
A more rigorous one is to use the neyman pearson lemma, that will allow you to calculate a thresold, but you need to know both the signal and noise probability densities.

The main problem is to estimate the signal and noise pdf, or equivalently under gaussian hypothesis, their mean and covariance. The main difficulties arrising from the fact that both aren't stationary and from the presence of signal biasing your estimation of noise's pdf, mean, variance etc ..., you'll want to go for robust estimator like median and median absolute deviation.


Deconvolution Using Response to an Heavy Side I'm measuring a "charge" signal in function of time from an amplifier.
Here is a measured signal (x-axis is the time in some arbitrary units, y-axis is the charge in ADU):

I would like to get the "true", original signal, by deconvolving the instrument's response from the measured signal.

Here the response of the instrument to an Heaviside step function:

I guess this should help me characterize the instrument's response, and do the deconvolution.
I read about SVD and Bayesian unfolding, but these don't seem to do the trick (I need to built a transfer matrix first).
Could you suggest me any tool/algorithm to do this?

Thanks in advance for the help.

If we can assume no noise (Or the SNR is very high) you can get the response by applying the inverse filter in frequency domain.

Lets say $ y [n] $ are the signal samples.
Given $ x [n] $ the samples of the ideal signal you can apply on both of them the DFT to get $ Y [k] $ and $ X [k] $.

The response is given by the Inverse DFT of the division $ \frac{Y[k]}{X[k]} $.

If there's some noise you need to regulate the result.
Easy choice would be the Wiener Filter or more specifically Wiener Deconvolution.


A filter for for a realtime system? Apologies if this is too basic.

I need to implement a band-pass filter, for a realtime system.
Not sure what to use. FIR? IIR? Something else ?

Any tips would be highly appreciated.

For an EEG signal, I suppose you need to keep the shape of the signal. Therefore the filter should have linear phase, so the reasonable way to go is a FIR filter.

If shape (transients) don't matter much for your application, then an IIR filter would be preferable since it will give less end-to-end delay in general.

Note that suprressing 0 - 1 Hz, inevitably gives a group delay of at least 1 second, depending on how sharp the transition should be, so your real-time requirements may conflict with the start frequency.

You may try some filter design tools (e.g. from matlab) to find a reasonable compromise between frequency response and group delay.


calculate SNR of two signals with power spectral density You can see my system's diagram as shown below

At output,two different input signals are sited in each ($\frac{H}{3}$) position $(w_2=w_1+ \frac{H}{3})$ and if output power spectral density equation is included different sentences such as

Dc term reaction of input noise to itself.
reaction of inputsignal(w1 or w2) to itself.
reaction of input signals(w1 and w2) to themselves
reaction of input signal(w1 or w2) and input noise
reaction of input signals(w1 and w2) and input noise to themselves.

How can obtain input and output SNR?

Since you only say your channel is nonlinear, but not exactly which math it follows, the only answer possible to give here is

1.-5.: By filling in the signals and noise in your device equation, and transforming the math until you get an expression for what you're looking for.

"nonlinear" doesn't really say anything about what's going on (aside from the transfer function not being linear; especially, questions like determination, time-variance, boundedness, memory... are left unanswered).

If you think about it, it can be shown that with the lack specification given, the statement above is the strictest statement possible.

Even incorporating the fact that you claim the power spectral density exists doesn't do us much good – it only tells us the output of the device is weak-sense stationary, and that might as well mean the thing is just generating random tones.


Modelling propagation of sound wave by particle simulation How might one model the propagation of a sound wave through a 3D environment with obstacles in it?

My initial thought would be to use freely moving particles that repel one another; with sufficiently many particles I see no reason why the simulation wouldn't work.

However, if any particle could collide with any other then this introduces a huge complexity hit.

So I'm wondering whether the particles could be effectively anchored on a sphere packed lattice (dodecahedron?), so each particle has 18 neighbours. (you could think of placing one tennis ball on the table, then six fitting around it, now 6 can be fitted on top, and by symmetry 6 could be fitted underneath. so the original ball has 18 neighbours.

And then shockwaves could propagate through this mesh, and at each iteration each particle just needs to calculate the resultant force of its 18 neighbours acting on it.

I can't see whether this model would work or not.  I suspect it is going to model sound propagating through a diamond rather than sound propagating through air, but isn't it going to be essentially the same?

Also, this isn't quite on topic for DSP, my apologies. Can anyone recommend a more suitable place to ask this question?

There is a huge body of research and commercial work on doing acoustic simulation. There are a two basic classes of simulation systems:

Finite element or finite boundary: these model the actual physics by dividing air and or/surfaces into small patches and modelling step by step the interaction between the patches by locally solving the wave equations given the boundary and initial conditions of the patches near by
Approximation methods: Ray tracing, particle tracing, mirror image method, etc. These are macroscopic models that use certain simplifying assumptions such as specular reflections, simplified diffusion models, no diffraction, etc.

This is a very complicated problem and it takes many person years to create a modelling system that actually produces useful results. There are commercial systems out there that you can buy. They vary greatly in complexity, capability, cost, and training required so it all depends on your specific application.


Unable to low pass filter the signal I am trying to use a basic example of a low pass filter as shown in the example:
https://de.mathworks.com/help/signal/ref/lowpass.html

However, I change the frequency:

fs = 20e6;
t = 0:1/fs:1;
x = [1 2]*sin(2*pi*[70e3 300e3]'.*t) + randn(size(t))/10;
lowpass(x,100e3,fs)

The following is the response that I get:

The location of the peaks are at 68 kHz and 300 kHz. I am not sure why 300 kHz is not filtered. Any guidance will be helpful.

See the documentation at https://www.mathworks.com/help/signal/ref/lowpass.html and in particular the "Steepness" parameter: You didn't specify a steepness, so it is just using the default value of 0.85 which is a transition width of 15% of $(f_{nyquist}-f_{pass})$, or 15% of (10MHz-100KHz) = 1.485 MHz. You see this directly with the shape of your post-filtered noise.

You could increase this to the minimum achievable that the lowpass function provides which is 1% (at the expense of filter complexity). This would result in a transition width of (10MHz-100KHz)/100 = 99 KHz which should result in significant attenuation of the 300 KHz signal.

Another approach would be to lower the sampling rate, which would be more efficient as long as there are not higher frequency signals present. This is often combined in a multi-rate approach where a more relaxed low pass filter is first used, then the sampling rate of the filtered signal is reduced (by simply taking every Nth sample after filtering where N depends on how much filtering is done) and then implementing the filter to your desired cut-off at the lower sampling rate.


Goertzel algorithm produces incorrect phase I've implemented Goertzel algorithm according to the Wikipedia page (https://en.wikipedia.org/wiki/Goertzel_algorithm) and another page (http://www.mstarlabs.com/dsp/goertzel/goertzel.html), which are consistent with each other. I am then testing it using a naive implementation for a single DFT frequency (https://en.wikipedia.org/wiki/Discrete_Fourier_transform). I am getting correct amplitude, but incorrect phase. Please help me find the mistake in my code, I've spent a couple of days on it now, and can't see the problem...

/**
 * Complex imaginary number i.
 */
template<typename T>
constexpr auto COMPLEX_i = std::complex<T> (
  static_cast<T> (0.0), static_cast<T> (1.0));

/**
 * Convert frequency in Hertz to Omega (radian per sample).
 */
template<typename T>
T
hertz_to_omega (T hertz, T sample_rate)
{
  return hertz / sample_rate * static_cast<T> (2.0) * M_PI;
}

/**
 * Simple implementation of the Goertzel algorithm
 * for detecting a given omega frequency from a range
 * of raw signal values.
 */
template<typename T, typename Range>
std::complex<T>
goertzel (const Range &range, T omega)
{
  const auto coeff = static_cast<T> (2.0) * std::cos (omega);

  T s1 = 0.0;
  T s2 = 0.0;
  for (const auto value : range)
  {
    const auto s0 = coeff * s1 - s2 + value;
    s2 = s1;
    s1 = s0;
  }

  return std::complex<T> (
    std::cos (omega) * s1 - s2,
    std::sin (omega) * s1) /
    static_cast<T> (range.size ());
}

BOOST_AUTO_TEST_CASE( test_goertzel )
{
  using T = double;

  const T FREQ = 543.2;
  const T SAMPLE_RATE = 5678.9;
  const T OMEGA = hertz_to_omega (FREQ, SAMPLE_RATE);
  const T START_RADIANS = M_PI / 3.0;
  const size_t LENGTH = 20000;

  std::complex<T> dft_sum (static_cast <T> (0.0), static_cast<T> (0.0));
  std::vector<T> data;
  {
    T radians = START_RADIANS;
    for (size_t sample_index = 0; sample_index < LENGTH; ++sample_index)
    {
      const T value = std::cos (radians);
      dft_sum += value * std::exp (- COMPLEX_i<T>
        * OMEGA * static_cast<T> (sample_index));
      data.push_back (value);
      radians += OMEGA;
    }
  }

  const auto result1 = dft_sum / static_cast<T> (LENGTH);
  const auto result2 = goertzel<T> (data, OMEGA);

  std::cout << "DFT: " << result1 << ", ABS: " << std::abs (result1) << std::endl;
  std::cout << "GOE: " << result2 << ", ABS: " << std::abs (result2) << std::endl;
}

Running the test above produces the following output:

Running 1 test case...
DFT: (0.250009,0.433004), ABS: 0.499997
GOE: (0.114572,0.486693), ABS: 0.499997

The magnitude is detected correctly (it should be 0.5), but the phases of the naive DFT implementation and the Goertzel algorithm output don't match. Please help me find the mistake in either my single-bin DFT or Goertzel implementation.

EDIT:

It seems that DFT calculation is detecting the phase correctly, as atan(0.433/0.25)/180 = 0.333 (and I specified the phase of pi / 3 for the generated signal on which I am testing), so its my implementation of Goertzel algorithm which is incorrect.

The short answer is that when Goertzel algorithm is used for omega, which doesn't correspond to of the finite set of frequencies computed by DFT (which correspond to integer k), then an additional correction factor should be applied. I've got an answer from Peter K after describing my question more clearly here: What's wrong with my Goertzel algorithm implementation?


What is the distance metric in dynamic time warping? I am reading the lecture slides from this site http://luthuli.cs.uiuc.edu/~daf/courses/CS-498-DAF-P/Lecture%2018%20-%20Time%20series,%20DTW.pdf

In slide 24, the lecture suggest using metric of the following

d(i, j) = ||f(i) - f(j)||

select top frequency in i and top frequency in j,
    find the Euclidean distance between the two

I am not exactly sure what do they mean by Euclidean distance between the two.

Did they meant to have top n frequency sorted in order and match the order for subtraction and find the euclidean distance?

For example

f1 = [1, 2, 3, 4]

f2 = [1, 100, 200 ,300]

d(1,2) = |1-1| + |100-2| + |200-3| + |300-4|

Any help would be appreciated

Thanks

The Euclidean distance between two vectors $\mathbf a$ and $\mathbf b$ is $\mathrm d(\mathbf a, \mathbf b) \equiv \sqrt {\sum_i (a_i - b_i)^2}$

In your example you forgot to square the differences. Remember the way you calculate hypotenuse of a triangle; it's the same idea.


Parallel form filter design I have a transfer function coefficients as following:
[b,a] = ellip(4,.2,40,[.41 .47]);

I'd like to factor this system into its PFE using [r,p,k] = residue(b,a). Then, I'd like to implement a filter which uses a parallel combination of second-order subsections. To do this, I need to combine each complex pole with its complex conjugate so that the overall second order subsections will have real valued coefficients. Hence; I selected four pairs of first order terms and recombined them into second order subsections using "residue" function. How can I store the coefficients of the resulting second-order subsections in the matrices c and d so that each row corresponds to one second order system?

I tried to do this but I don't know whether it's correct or not?

[c1 d1]=residue([r(1) r(2)],[p(1) p(2)],0);
[c2 d2]=residue([r(3) r(4)],[p(3) p(4)],0);
[c3 d3]=residue([r(5) r(6)],[p(5) p(6)],0);
[c4 d4]=residue([r(7) r(8)],[p(7) p(8)],k);

[c d]=[c1 c2 c3 c4 d1 d2 d3 d4]

Can anybody please give me an idea? Any help would be appreciated

A few points first

Since you are using digital filters, you should be using residuez() not residue()
Matlab's residuez() is not particularly well implemented: in case of higher order or steep filters, you can get a lot of numerically noise or results that are just totally wrong. Always check the result against your original target.

Back to the question: each residue comes as a conjugate complex pair (at least if you are lucky and don't have to deal with real or multiple poles). Just add these to form a single real second order section.

$$H(z) = \frac{r}{1-p \cdot z^{-1}} + \frac{r^{*}}{1-p^* \cdot z^{-1}} $$
$$= \frac{r \cdot (1-p^* \cdot z^{-1})+r^* \cdot (1-p^* \cdot z^{-1})}{(1-p \cdot z^{-1}) \cdot (1-p^* \cdot z^{-1})}$$

$$= \frac{(r + r^*) - (r \cdot p^* + r^* \cdot p  ) \cdot z^{-1} }{1 -(p+p^*)\cdot z^{-1}+ p \cdot p^* \cdot z^{-2}}$$

The last equation  gives you directly the real valued biquad coefficients (if I did the math correctly,  that is).


Program to display amplitude of sine signal at each interval and display separately in MATLAB I am facing a problem on how to separate amplitude of sine signal and compare with other sine signal having same frequency but different amplitude. How can one write a code for this?

it's actually easy, you only need to find the maximum value of your sine, which is in turn happens to be your amplitude, then you can compare them with if an if statement.

this gives you the amplitude

max(signal);

use if statement to compare the two values

if max(signal1) == max(signal2)
%% your condition comes here
end

Hope this helps.

cheers.


Why two images with same cluster value have different representation? k-mean cluster classifies pixel randomly. I wrote this code to repeat clustering until both images meet proper cluster.

 while 1
  idx_fix = kmeans(double(fix(:)),4);
  idx_float = kmeans(double(float(:)),4);
  idx_fix = reshape(idx_fix,size(fix));
  idx_float  = reshape(idx_float,size(float));
  f1 = idx_fix(127,303);
  f2 = idx_float(114,205);
  if (f1 == f2)
      disp(idx_fix(127,303));
      disp(idx_float(114,205));
      if (f1 == f2)%when these pixles same that mean this is my target classification.
          break;
      end;
  end
  end
  mat1 = uint8(255*mat2gray(idx_fix));
  mat2 = uint8(255*mat2gray(idx_float));
  disp(['mat1 = ' num2str(mat1(127,303))]);
  disp(['mat2 = ' num2str(mat2(114,205))]);% the value still correct after coversion.
  imshow(mat1); title('Fix image class');
  hold on
  plot(127,303,'or');
  hold off;
  figure;
  imshow(mat2); title('Float image class');
  hold on
  plot(114,205,'or');
  hold off;

  similarity_chart(idx_fix,idx_float);

The output image is not as expected.
I point the compared pixel by a red circle, however, the area around those circles should have the same color... why this happened?

The original images -> https://1drv.ms/f/s!Al4zZ4gTzuGXhOtUnqo7Oiui1HMpSg

Frankly speaking I did not quite understand your question, but let me write my comments to this kind of problems. This is the signal processing forum not only MatLab, so I will not discuss the code. All calculations I did in Mathematica.

First of all you cluster pixel intensity, meaning that you clustering 1D array and spatial correlation don't make sense and it is better to deal with the image histogram.

First let me plot your image histogram.

As you can see the left spike corresponds to the background. It is not a big problem, but it's better to eliminate the background from the clustering problem.

Using various filters I create the mask for your image,

and check that the mask is correct.

Let's plot the image histogram without the background (please note, I plot the smoothed histogram).

Now you can start clustering of your histogram. If you know a priory the number of clusters, you can do Kmean clustering, for example. Kmean simply finds the local minimum of the variation of your data from clusters. You may have many local minimums and you have to vary initial conditions to find the global one.

Let me give you a simple example.

Let's consider a mixture of 4 Gaussian distributions

where mean values are {10, 15, 25, 35}, standard deviations are {1, 1, 1, 1} and amplitudes {1/4, 1/4, 1/4, 1/4}. You can write the problem of clustering of the mixture on two cluster as minimization problem of function of two variables (because you know the distribution analytically). You can plot this function

where the x-variable corresponds to the boundary between the first and the second classes, and y-variable between the second and the third classes.

As you can see on the plot, there are two minimums, where the top-right minimum is the global one. Let's plot these minimums on the mixture PDF

where red lines show the global minimum and green lines the local one. Even for such a simple problem Kmean may converge to the local solution, and you should do the best to be sure that the method converges the global one.

If you have 1D data, finding the global minimum can be done by a simple brute-force. You can split the histogram on all possible classes, calculate the variation and find the global solution. I took your histogram and brute-force it to separate on three classes. Please, find the result below

where the best combination is

{{1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, {11, 12, 13, 14, 15, 16, 17}, {18,
  19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
  36, 37, 38, 39, 40, 41, 42, 43, 44, 45}}

and the histogram is

intensity={0., 0.0204082, 0.0408163, 0.0612245, 0.0816327, 0.102041, 0.122449,
0.142857, 0.163265, 0.183673, 0.204082, 0.22449, 0.244898, 0.265306,
0.285714, 0.306122, 0.326531, 0.346939, 0.367347, 0.387755, 0.408163,
0.428571, 0.44898, 0.469388, 0.489796, 0.510204, 0.530612, 0.55102,
0.571429, 0.591837, 0.612245, 0.632653, 0.653061, 0.673469, 0.693878,
0.714286, 0.734694, 0.755102, 0.77551, 0.795918, 0.836735, 0.857143,
0.897959, 0.938776, 0.959184}

bins={0.00766421, 0.0344291, 0.0280162, 0.0230938, 0.0201404, 0.0210605,
0.0252008, 0.0263049, 0.0271238, 0.0306568, 0.0330398, 0.0265901,
0.0287339, 0.0391123, 0.0618841, 0.0936635, 0.091078, 0.0935254,
0.120474, 0.0829354, 0.0292859, 0.0114089, 0.00848308, 0.00669813,
0.00557564, 0.0050144, 0.00423234, 0.00348708, 0.00262221,
0.00204256, 0.00185855, 0.0012237, 0.000901672, 0.000726858,
0.000634851, 0.000414033, 0.000266821, 0.000184015, 0.0000736059,
0.0000828066, 9.20073*10^-6, 9.20073*10^-6, 9.20073*10^-6,
9.20073*10^-6, 0.0000184015}

If you don't know the number of classes, you may use clustering which can detect it (DBSCAN, Mean shift, Hierarchical clustering, variational Gaussian mixture algorithm, spectral clustering algorithm
 just to mention a few). For your data the best looking result gives variational Gaussian mixture algorithm, at lest it gives 3 clusters. Other method gives more than 10 clusters or less than two. Please find below result of clustering with the variational Gaussian mixture algorithm.

For my view, the result is still not good enough. If you histogram is the mixture of highly-overlapped distributions, this type of clustering may fail. In this scenario clustering of the histogram in scale-plane may help. I use idea from the  paper, with my modifications. The basic of the idea tells that you take your histogram, convolve it with the Gaussian kernel increasing the scale parameter (the standard deviation) and track all local minimums of the histogram (local minimums of the function). Then, you may find that some minimums survive longer than others and survivors minimums are the clusters boundary. You can find the plot which shows the positions of minimums depending on the scale parameters in the convolution for your histogram

Then, using Kmean clustering, you can separate long lines from the short ones (you can find details in the paper). The separation gives

And finally it shows, that there are three minimums, where the right minimum looks like the boundary artifact. These minimums are the boundaries of your classes. Please, find the plot the boundaries

which looks reasonably good. Please note I slightly modified the method from the paper, if you need more details, just let me now.


FFT with day-periodic data I have a collection of data: number of hours per day, over the course of several years.  I want to create an FFT of this data.  Since my sampling rate isn't seconds, but days, I think I have two choices when constructing the FFT.  One, I can specify that the sample rate Fs=1/86400 samples per second (86400 seconds in a day.)  But I think that's not right.  So maybe I can: two, say Fs=7 samples per week, then my resulting FFT has a scale, instead of cycles-per-second, of cycles-per-week.  Is this reasonable, or is there a better way to go about doing FFT for daily-periodic data?

Saying "7 samples per week" is perfectly reasonable.  The FFT doesn't care how the rate is defined, it only requires it to be constant.

When I fiddle with really low rates, I use the same concept.  Say, samples per minute or samples per year.

Mathematically, it doesn't matter.  You know yourself that the number 1/86400 samples per second is mathematically the same as 7 samples per week - you calculated them both, after all.

The only effect it will have is in how you translate the FFT bins into a frequency - and you have to do that anyway, regardless of how you represent the sampling rate.

For myself, I find it easier to think of the numbers and samples when represented as some whole number of samples per time unit.


Why is the Fourier transform so important? Everyone discusses the Fourier transform when discussing signal processing. Why is it so important to signal processing and what does it tell us about the signal?

Does it only apply to digital signal processing or does it apply to analog signals as well?

This is quite a broad question and it indeed is quite hard to pinpoint why exactly Fourier transforms are important in signal processing. The simplest, hand waving answer one can provide is that it is an extremely powerful mathematical tool that allows you to view your signals in a different domain, inside which several difficult problems become very simple to analyze.

Its ubiquity in nearly every field of engineering and physical sciences, all for different reasons, makes it all the more harder to narrow down a reason. I hope that looking at some of its properties which led to its widespread adoption along with some practical examples and a dash of history might help one to understand its importance.

History:

To understand the importance of the Fourier transform, it is important to step back a little and appreciate the power of the Fourier series put forth by Joseph Fourier. In a nut-shell, any periodic function $g(x)$ integrable on the domain $\mathcal{D}=[-\pi,\pi]$ can be written as an infinite sum of sines and cosines as

$$g(x)=\sum_{k=-\infty}^{\infty}\tau_k e^{\jmath k x}$$
$$\tau_k=\frac{1}{2\pi}\int_{\mathcal{D}}g(x)e^{-\jmath k x}\ dx$$

where $e^{\imath\theta}=\cos(\theta)+\jmath\sin(\theta)$. This idea that a function could be broken down into its constituent frequencies (i.e., into sines and cosines of all frequencies) was a powerful one and forms the backbone of the Fourier transform.

The Fourier transform:

The Fourier transform can be viewed as an extension of the above Fourier series to non-periodic functions. For completeness and for clarity, I'll define the Fourier transform here. If $x(t)$ is a continuous, integrable signal, then its Fourier transform, $X(f)$ is given by

$$X(f)=\int_{\mathbb{R}}x(t)e^{-\jmath 2\pi f t}\ dt,\quad \forall f\in\mathbb{R}$$

and the inverse transform is given by

$$x(t)=\int_{\mathbb{R}}X(f)e^{\jmath 2\pi f t}\ df,\quad \forall t\in\mathbb{R}$$

Importance in signal processing:

First and foremost, a Fourier transform of a signal tells you what frequencies are present in your signal and in what proportions.

  Example: Have you ever noticed that each of your phone's number buttons
  sounds different when you press during a call and that it sounds the same for every phone model? That's because they're each composed of two different sinusoids which can be used to uniquely identify the button. When you use your phone to punch in combinations to navigate a menu, the way that the other party knows what keys you pressed is by doing a Fourier transform of the input and looking at the frequencies present.

Apart from some very useful elementary properties which make the mathematics involved simple, some of the other reasons why it has such a widespread importance in signal processing are:

The magnitude square of the Fourier transform, $\vert X(f)\vert^2$ instantly tells us how much power the signal $x(t)$ has at a particular frequency $f$.
From Parseval's theorem (more generally Plancherel's theorem), we have
$$\int_\mathbb{R}\vert x(t)\vert^2\ dt = \int_\mathbb{R}\vert X(f)\vert^2\ df$$
which means that the total energy in a signal across all time is equal to the total energy in the transform across all frequencies. Thus, the transform is energy preserving.
Convolutions in the time domain are equivalent to multiplications in the frequency domain, i.e., given two signals $x(t)$ and $y(t)$, then if

$$z(t)=x(t)\star y(t)$$
where $\star$ denotes convolution, then the Fourier transform of $z(t)$ is merely

$$Z(f)=X(f)\cdot Y(f)$$

For discrete signals, with the development of efficient FFT algorithms, almost always, it is faster to implement a convolution operation in the frequency domain than in the time domain.
Similar to the convolution operation, cross-correlations are also easily implemented in the frequency domain as $Z(f)=X(f)^*Y(f)$, where $^*$ denotes complex conjugate.
By being able to split signals into their constituent frequencies, one can easily block out certain frequencies selectively by nullifying their contributions.

  Example: If you're a football (soccer) fan, you might've been
  annoyed at the constant drone of the vuvuzelas that pretty much
  drowned all the commentary during the 2010 world cup in South Africa.
  However, the vuvuzela has a constant pitch of ~235Hz which made it
  easy for broadcasters to implement a notch filter to cut-off the
  offending noise.[1]

A shifted (delayed) signal in the time domain manifests as a phase change in the frequency domain. While this falls under the elementary property category, this is a widely used property in practice, especially in imaging and tomography applications,

  Example: When a wave travels through a heterogenous medium, it
  slows down and speeds up according to changes in the speed of wave
  propagation in the medium. So by observing a change in phase from
  what's expected and what's measured, one can infer the excess time
  delay which in turn tells you how much the wave speed has changed in
  the medium. This is of course, a very simplified layman explanation, but
  forms the basis for tomography.

Derivatives of signals (nth derivatives too) can be easily calculated(see 106) using Fourier transforms.

Digital signal processing (DSP) vs. Analog signal processing (ASP)

The theory of Fourier transforms is applicable irrespective of whether the signal is continuous or discrete, as long as it is "nice" and absolutely integrable. So yes, ASP uses Fourier transforms as long as the signals satisfy this criterion. However, it is perhaps more common to talk about Laplace transforms, which is a generalized Fourier transform, in ASP. The Laplace transform is defined as

$$X(s)=\int_{0}^{\infty}x(t)e^{-st}\ dt,\quad \forall s\in\mathbb{C}$$

The advantage is that one is not necessarily confined to "nice signals" as in the Fourier transform, but the transform is valid only within a certain region of convergence. It is widely used in studying/analyzing/designing LC/RC/LCR circuits, which in turn are used in radios/electric guitars, wah-wah pedals, etc.

This is pretty much all I could think of right now, but do note that no amount of writing/explanation can fully capture the true importance of Fourier transforms in signal processing and in science/engineering


Inverse DFT to a real signal when only odd harmonics are present There is an old paper by L.R. Rabiner, "On the Use of Symmetry in FFT Computation" which describes (among other things) an optimized method to calculate the DFT of a real signal $x_n$, $n=0,\dots, N-1$ with the property that $x_{N/2+n} = -x_{n}, \ n=0,\dots, N/2-1$. The resulting transform $X=\mathrm{FFT}[x]$ only contains odd harmonics, that is $X_k=0$ when $k \,\mathrm{mod}\,2=0$, and $X_k=X_{N-k}^*$. The method requires an $N/4$ complex-valued FFT (or, equivalently, a $N/2$ real-to-complex FFT) with some pre- and post-processing.

What I wonder is if the corresponding inverse transform is described somewhere, that is obtaining the first half of $x_n$ given the first half of nonzero $X_k$. Unfortunately, Rabiner only provides an algorithm for the forward transform, and I cannot figure it out myself.

So, I've read the paper I cited more attentively, and it turned out that the author, while not having described the inverse algorithm explicitly, left enough clues in the derivation of the forward algorithm, after all. In case anyone needs it in the future, here it is.

Input: the nonzero (odd) harmonics $X_k$, $k=0,\dots,N/4-1$ (assuming here for simplicity that $N$ is divisible by 4).

Build the array $Y_k$, $k=0,\dots,N/4$: $Y_0 = -2\mathrm{Im}X_0$, $Y_{N/4} = 2\mathrm{Im}X_{N/4-1}$, for the rest $Y_k = i(X_{k+1}-X_k)$.
Calculate $y = \mathrm{IRFFT}(Y)$ (the standard inverse FFT to real signal, e.g. numpy.fft.irfft in Python or dft_c2r from FFTW). The resulting vector $y$ is real and has length $N/2$.
$x_0 = \frac{2}{N}\sum_{k=0}^{N/4-1} \mathrm{Re}X_k$, $x_n = y_n / (4 \sin(2\pi n))$, $n=1,\dots,N/2-1$. This gives the first half of the full vector $x$.

Of course, this method is not ideal, since the division by a sine can introduce instabilities. There may exist a method to deal with that (there is one for the FFT of an (anti)symmetric signal, which requires a similar post-processing).


How to mix two signals without changing the overall loudness? I have two audio signals that I want to mix at various mixing ratios. Initially, I simply went for something like

$y(t) = \alpha \cdot x_1(t) + (1-\alpha) \cdot x_2(t)$
where $\alpha$ is the ratio between zero and one and $x_1$ and $x_2$ are the two signals.

However, I have seen a few panning curves that did not mix signals linearly. They told me that a linear mix would change the overall loudness while mixing.

Is there something to that or is a linear mix the correct way to do it?

That depends somewhat on the application. In most cases a "constant energy" pan will be best. This can be expressed as $$y(t)=\sqrt{\alpha} \cdot x_{1}(t)+\sqrt{1-\alpha}\cdot x_{2}(t)$$ where $\alpha = .5$ is the point of equal energy.

If you are working with fixed point signals, such as wave files for examples, you may run into clipping problems. That could be a topic for a separate question.


How to plot the poles and zeros of this Bandpass Filter in MATLAB? So I'm trying to design a band pass filter in MATLAB (with a center frequency of 10kHz and a sampling frequency of 44kHz).

I have calculated the transfer function but I'm not sure how to enter this transfer function into MATLAB to generate a plot of the frequency and phase responses of the filter, and the poles and zeros.

It could also be that the transfer function is calculated wrong? I haven't seen the denominator have so many terms before so I'm not sure..

This is my first time posting on SE, so any help is appreciated.

PS: I kinda always thought that the poles of a band pass filter would be pi/2 always, so I guess I was wrong?

You've designed a second-order band pass filter, so the number of terms in the denominator is just fine (it's a second-order polynomial). Your design is OK, apart from the scaling, i.e., the maximum of the magnitude of the transfer function is not equal to $1$, which might be desirable.

Your transfer function has the form

$$H(z)=\frac{z^2-1}{(z-p)(z-p^*)}=\frac{z^2-1}{z^2-2zr\cos(\omega_p)+r^2}\tag{1}$$

where $p$ is a complex-valued pole:

$$p=re^{j\omega_p}\tag{2}$$

with radius $r$ and pole angle $\omega_p$. Note that in general the pole angle $\omega_p$ and the peak frequency $\omega_0$ (where the magnitude of the frequency response $H(e^{j\omega})$ achieves its maximum) are not equal. However, in your case they are almost identical because $\omega_p$ is relatively close to $\pi/2$. See this answer for more details on this matter.

If you want the frequency response peak value to be equal to $1$, you need to scale your transfer function by the factor $(1-r^2)/2$ (also explained in the answer quoted above).

The pole angle $\omega_p$ depends on the desired center frequency, so it can't always be equal to $\pi/2$ as you seem to have thought (even though in your calculations you chose it according to the desired center frequency).

The Matlab function freqz computes the complex frequency response of a discrete-time filter.
For plotting the poles and zeros you can use the function zplane.


DFT of real sinusoids - why sum over -$N/2$ to $N/2-1$ as opposed to $0$ to $N-1$? I'm going through a Coursera course on signal processing, and we're just introduced to DFTs.

We are told that if you have a complex sinusoidal signal
$x[n]$ where $n=0,1...\ N-1$, its DFT is given as

$$X[k]=\sum_{n=0}^{N-1} x[n] e^{-j2\pi kn/N}$$

And this makes sense to me. Because $x[n]$ has $N$ samples. You take each one, multiply it with the sample of a basis function, essentially taking the dot product of the signal with the basis function, and you get how much of that frequency is present in the input signal.

But we are told that if you have a real sinusoidal signal $y[n]$ where $n=0,1...\ N-1$, then the DFT is given as

$$X[k]=\sum_{n=-N/2}^{N/2-1} y[n] e^{-j2\pi kn/N}$$

My question: Why do we sum over the different values? What is the intuition and explanation behind this?

I'm not even sure what $y[-N/2]$ would mean.

Is the $-N/2^{th}$ sample the same as the $N/2^{th}$ sample because we assume the input signal is periodic? Even if so, why is taking the DFT summation over such a range beneficial? We could also do this for complex sinusoidal signals since they're periodic too right? What's the advantage?

The DFT (Discrete Fourier Transform) of an $N$-point sequence $\left\{x[n]\right\}$ is given as:
$$X[k]=\sum_{n=0}^{N-1} x[n] e^{-j2\pi kn/N}, \quad k=0, 1, \ldots, N-1\tag{1}$$

If we assume the sequence $\left\{x[n]\right\}$ and its DFT are complex, they can be expressed as follows (with the subscript $r$ and $i$ for real and imaginary part respectively):

\begin{align}
x[n] &= x_r[n]+jx_i[n], \quad 0\leq n\leq N-1\\
X[k] &= X_r[k]+jX_i[k], \quad 0\leq k\leq N-1
\end{align}

With a bit of algebra and using equation $(1)$, you get:

\begin{align}
X_r[k] &= \sum_{n=0}^{N-1}\left[ x_r[n]\cos\left(\frac{2\pi kn}{N}
\right)+x_i[n]\sin\left(\frac{2\pi kn}{N}
\right)\right]\\
X_i[k] &= -\sum_{n=0}^{N-1}\left[ x_r[n]\sin\left(\frac{2\pi kn}{N}
\right)-x_i[n]\cos\left(\frac{2\pi kn}{N}
\right)\right]
\end{align}
If the sequence $\left\{x[n]\right\}$ is real-valued, this means $ x[n]=x_r[n]$ as $x_i[n]=0$. The DFT will then be expressed as:

$$X[k]=\sum_{n=0}^{N-1} x[n]\left[\cos\left(\frac{2\pi kn}{N}
\right)-j\sin\left(\frac{2\pi kn}{N}
\right)\right]\tag{2}$$

And it can be shown from equation $(1)$ and $(2)$ that for these real-valued signals:

$$X(N-k)=X^*(k)=X(-k)\tag{3}$$

From the above results, it suffices to compute the DFT for half of the values of $k$ to get the full characteristic of the spectrum. From the $N$ values you compute less (roughly half) $k$ values. The boundaries of the summation are still the same, it still goes from $n=0$ to $N-1$ but your $k$'s shouldn't go up to $N-1$ as it will be a waste of resources.
$$X[k]=\sum_{n=0}^{N-1}y[n] e^{-j2\pi kn/N}$$

For $N$ even, $ 0\leq k \leq N/2$ suffices. With the results from equation $(3)$ you can get the remaining values from $N/2 + 1$ to $N-1$ or from $-N/2 + 1$ to $-1$.
You have $ 0\leq n \leq N$. So, $y[-N/2]$ does not exist and $y[N/2]$ is the value of your sequence $\left\{y[n]\right\}$ at $n=N/2$ (if $N$ is even).
However, for a real even sequence $\left\{y[n]\right\}$, the DFT $X[-N/2]$ equals $X[N/2]$ and this is the value at Nyquist.
Also note from equation $(3)$ that if $N$ is even, $X[0]$ and $X[N/2]$ must be real, and if $N$ is odd $X[0]$ is still real but you don't have a frequency bin at exactly $N/2$. Equation $(3)$ applies for real-valued sequences, what has been concluded from that equation doesn't apply for strictly complex-valued sequences. The advantage is that you've saved computations by computing only half of what you need to compute.


What is the difference between linear and non-linear filters? How a mean filter is called as linear filter and a median filter is called as non linear filter? I understand how a mean and median filter operates, but I was not able to relate with the term linear and non-linear. Please explain me with an example.

To understand the linearity property more easily.Let us consider the above diagram,here we have 2 sequences namely Xn and Yn. when we add both the sequence we get Xn+Yn whose amplitude value are represented with blue colour. when any system which satisfy this condition then it is called linear.
In case of mean filter, mean value for sequence Xn is 1+1+3/3=5/3.mean value for sequence Yn is 1+2+0/3=1.mean value for Xn+Yn is 2+3+3/3=8/3.

$ mean(Xn)+mean(Yn)= mean(Xn+Yn), \\ \\
 5/3    +    1   = 8/3 $

hence we called mean filter as linear filter.
In case of median filter, if we calculate median value for sequence Xn, we get 1 (arrange the sequence in ascending order and then find the middle value). similarly median value of sequence Yn is 1 .median value of sequence Xn+Yn is 3.

$median(Xn)+median(Yn)  \neq  median(Xn+Yn).\\ \\
1 + 1  \neq    3 $

hence we call median filter as non-linear filter


Is it possible to Denormalize a matrix of complex numbers? In this article, FFT and I-FFT are achieved as the following:

FFT steps

process starts with a Bitmap image.

Bitmap-image to Integer-image
Integer-image to Complex-image
Apply Forward-FFT to Complex-image (output: Complex-FFT)
Apply Shift to FFT (output: Shifted Complex-FFT)

I-FFT steps

process starts with the FFT (a complex matrix) of an image.

Apply Backward FFT to Complex-FFT (output: Complex-image)
Complex-image to Integer-image
Integer-image to Bitmap-image

Normalization

To display an FFT's Magnitude/Phase, FFT is Normalized and then that normalized matrix of complex numbers is converted to a Bitmap. Without Normalization, I see a complete 'Black' image output.

My Question

As I posted this question, I was thinking about applying an I-FFT in the following process:

Obtain the Bitmap-image of a magnitude plot of an FFT.
De-normalize the Bitmap-image to obtain the original FFT of the image.
Apply I-FFT to the complex matrix and obtain the Complex-image.
Convert the Complex-image to Integer-image.
Convert the Integer-image to the original Bitmap.

N.B. For some reason, they used,

(1) Bitmap $\Rightarrow $ Integer$\Rightarrow $Complex, and, (2) Complex$\Rightarrow $Integer$\Rightarrow $Bitmap

to obtain the conversion between a Bitmap and its Complex counterpart.

Obtain the Bitmap-image of a magnitude plot of an FFT.

  De-normalize the Bitmap-image to obtain the original FFT of the image.

This can't work. Preserving only the magnitude throws away half of the information contained in the complex values, and thus, you can't get the original FFT back, ever.


How do I implement IFFT to filter low frequencies from a group,after Freq detection by FFT? I use a combination of Visual Basic Programming interfaced with FAMOS which is a DSP software like MATLAB but much easier. I am trying to separate signals from a group, which requires filtering each signal. I used FFT and then IFFT on a 600s signal but the IFFT was longer than the original.
I don't understand if this changes the originality of the signal in terms of shape,which might affect the further processing of the signal.
After re-sampling the original 10min signal according to Shannon-Niquist theorem and using rectangular windowing, I detect the frequency using FFT.

After I apply IFFT to the signal, does this change the phase or amplitude even slightly that might vary from the original signal?
How do I filter one particular freq before IFFT after detecting it from FFT (Maybe I use peak detection)? Is there a general method of doing it?
DSP is totally new to me,and I am still attempting to learn.

Frequency domain processing (i.e FFT, manipulating signal, IFFT) is actually pretty complicated and requires careful choice and a thorough understanding of things like framing, step size, windowing, overlap, wrap-around, circular vs. linear etc. Simple examples for these type of algorithms are overlap-add or overlap-save for linear time invariant filtering.

If the only thing you want to do is getting rid of the low frequencies in your signal, filtering with a high pass filter is much easier and quicker.


Correct method for up-sampling cross-power spectrum for sub-pixel motion estimation Let's say I have two images, one reference($g$) and another shifted($f$). I tried to get subpixel accuracy of the motion. Here are my steps:

Firstly I would calculate the normalized cross-power spectrum $C(u,v)$, where $F$ and $G$ are fft2-ed.

$$C(u,v)=\dfrac{F(u,v)G^*(u,v)}{\lvert F(u,v)G^*(u,v) \rvert} $$

Next I would like to get $1/k$ subpixel accuracy, by padarray the $C(u,v)$ by $k$ factor, resulting in $C_k(u,v)$
ifft2 $C_k(u,v)$, resulting in $c_k(x,y)$
Find the arg max of the $c_k(x,y)$

Here's my code in Matlab:

reference = fft2(reference);
shifted = fft2(shifted);
C=shifted.*conj(reference)./abs(shifted.*conj(reference));
[m, n]=size(reference);
ck=ifft2(padarray(C,[k*m k*n],'both'));
[~,w] = max(real(ck(:)));
[dxk,dyk] = ind2sub(size(ck),w);

After run the code I can't get the accurate result. I hope I explained it clear enough(from what I understood from the attached paper). I really need help on this and much appreciate any help!

Referred paper: Here (start from section 3)

Preface

I know the OP asked for Matlab, but if you follow our discussion here in the Mathematica chat it seems any answer will be helpful. Therefore, I'm going to present the approach, visualizing each step, so it can be repeated and debugged in Matlab.

Answer

The approach itself is simple and as far as I can see there are only two things to keep in mind. But first things first. I'm using a downscaled version of the Lena image that can be accessed in Mathematica through ExampleData. I shifted it by {-5.3, 8.7} pixels which can be done through ImageTransformation in combination with TranslationTransform.

{dx, dy} = {-5.3, 8.7};

gImg = ImageResize[
   ColorConvert[ExampleData[{"TestImage", "Lena"}], "Grayscale"], 128];
{nx, ny} = ImageDimensions[gImg];
fImg = ImageTransformation[gImg, TranslationTransform[{dx, dy}/{nx, ny}]];

{gData, fData} = ImageData /@ {gImg, fImg};

Keep in mind that shifting the image by fractions of pixel units will resample it and introduce some blur that degrades the image. Therefore, I would not expect that the perfect solution can be found.

What we have in gData and fData are the pixel values of the images which look like this:

Fourier transform

I don't know how fft2 in Matlab exactly works, but usually, taking the Fourier transform of an image does not result where the frequencies are centered. Let me give you an example, by showing how the spectrum looks if you just use FourierTransform

Image[Fourier[gData] // Abs]

To make the upscaling by padding work, we need to center the transform. You can do this, by multiplying each pixel by $(-1)^{x+y}$ which gives a chessboard-like fashion. This is explained in detail in the very easy to read Digital Imageprocessing by Gonzalez and Woods. A pdf of the 2nd edition can be found here and you should follow the discussion on page 154.

In Mathematica, a function to transform an image into this chessboard-like form can be written down as

centerFrequency[data_] := With[{d = Dimensions[data]},
   Table[(-1)^(i + j)*data[[j, i]], {j, Last[d]}, {i,
      First[d]}]
   ];

If we now look at the same spectrum, with this function applied before using the FFT, we get the expected result:

Image[Fourier[centerFrequency@gData] // Abs]

Normalized Cross Power Spectrum

To create an upscaled version of the normalized cross-power spectrum, we follow exactly what is written in the paper

Calculate the centered fourier transform of fData and gData into f and g
Conjugate g
Calculate f*g/Norm[f*g]
Pad the result by a factor k of the original image dimension (I assume a square image in this example)

In code this looks like this

normCrossPow[ref_, tmpl_, k_] := Module[{
  f = Fourier[centerFrequency@ref],
  g = Conjugate[Fourier[centerFrequency@tmpl]], nx
  },
  nx = First[Dimensions[ref]]/2;
  ArrayPad[f*g/Norm[f*g], k*nx]
]

To get the resulting spectrum, we only need to apply InverseFourier, apply centerFrequencies again and look at the real part of the complex array.

spec = Re[
   centerFrequency@InverseFourier[normCrossPow[fData, gData, 3]]
];
Image[spec] // ImageAdjust

Above, I used a factor of k=3.
The white spot in the upper right corner is our maximum we are interested in. Let us find the position of this maximum

max = FirstPosition[spec, Max[Flatten[spec]]]

This gives {36, 22}. Keep in mind, that the y-axis in image coordinates and matrix-coordinates is reversed and that in a matrix the y-values comes first. Additionally, if images have no shift at all, the maximum would be at {1,1} and we need to divide out our factor k again. Note that we divide by 4 (!!) because with a factor of k=1 we would divide by 1, with k=1 by 2, and so on.

({36, 22} - {1, 1})/4.0

This gives {8.75, 5.25} for {dy, dx} but remember that it is actually dy=-8.75 because of the reversed y-axis.

Now, we can use this result and back-transform our fImg

res = ImageTransformation[fImg, TranslationTransform[{5.25, -8.75}/{nx, ny}]]

and look at the combined image of gImg and res

Fits nicely and the only error we see is the black border that was introduced by initially transforming fImg.


Recover an OFDM signal affected by AWGN I have an OFDM signal affected by AWGN, if I know the power of the gaussian noise, can I recover the OFDM signal? Can I perform deconvolution to recover the original signal?
thanks

Yes, when we can recover OFDM in fading scenarios along with gaussian noise, then yes we can also recover it in just gaussian noise. You can think of fading coefficients having a magnitude of 1 and phase zero


How do I interpret this block diagram correctly? I am reading notes ahead of the class and have encountered this particular slide:

While I completely agree with the first block diagram, I am at a loss trying to understand how the second block diagram is equivalent to the first.

Here's how I interpret the first diagram:
x[n] enters H1 and H1 outputs b0 x[n] + b1 x[n-1].
The output of H1 enters H2, and H2 outputs -a1 y[n-1] + b0 x[n] + b1 x[n-1]

Here's how I interpret the second diagram:
I am actually lost! I have no idea how to make sense of it.

I understand how the "simpler" diagrams work, for example the multiplication and addition components. But this is Greek to me, and it seems like I cannot find any list of rules to interpret this.

It is better to "parse" these networks from the output back towards the input, calling their input some general $x$ and performing substitutions and/or compositions.

So, let's call these networks $U$pper and $L$ower.

From the upper diagram:

$$UH_2[n] = x[n] + x[n-1] \cdot -a_1$$

and

$$UH_1[n] = x[n] \cdot b_0+x[n-1] \cdot b_1$$

Now, the output of $U$pper is given by the composition of the two:

$$UH_y[n] = UH_2[UH_1[n]]$$

This is because, the output of one, becomes the input to the other.

Or...

$$UH_y[n] = UH_1[n] + UH_1[n-1] \cdot -a_1$$

...and if you substitute...

$$UH_y[n] = x[n] \cdot b_0 + x[n-1] \cdot b_1 + (x[n-1] \cdot b_0 + x[n-2] \cdot b_1) \cdot -a_1 \Rightarrow \\ x[n] \cdot b_0 + x[n-1] \cdot b_1 + x[n-1] \cdot b_0 \cdot -a_1 + x[n-2] \cdot b_1 \cdot -a_1$$

Notice here that if you are trying to get to the $-1$ of the $x[n-1]$, then that would be the $x[n-2]$.

And this concludes the $U$pper part.

For the $L$ower part, we are not going to go through the whole thing, because of two reasons:

If you notice, the $H2, H1$ are reversed (The $L$ower network calls $H2$ what the $U$pper network calls $H1$).
We have already mapped $H2, H1$.

So, the lower network's response is:

$$LH_y[n] = LH_1[LH_2[n]]$$

Or...

$$LH_y[n] = LH_2[n] \cdot b_0+LH_2[n-1] \cdot b_1$$

And if you substitute:

$$LH_y[n] = (x[n]+ -a_1 \cdot x[n-1]) \cdot b_0 + (x[n-1]+ -a_1 \cdot x[n-2]) \cdot b_1 \Rightarrow \\
   x[n] \cdot b_0 + -a_1 \cdot b_0 \cdot x[n-1] + x[n-1] \cdot b_1 + b_1 \cdot -a_1 \cdot x[n-2]$$

Which, after you re-arrange, looks exactly the same.

Hope this helps.


Covariance between real and imaginary parts of Fourier transform of a stationary time series Since Fourier transform of a random stationary process in time (in the case of existence) is not necessarily real, my question is what is the relation between the covariance of real and imaginary parts of the transformed random variable.

More specifically, according to Wiener-Khinchin theorem for the time function $x(t)$ with some conditions the Fourier transform $X(\omega)$ exists.

The question is can one recover  $$\mathbb{E}\left[\Re\{X(\omega)\},\Re\{X(\omega')\}\right]-\mathbb{E}\left[\Re\{X(\omega)\}\right]\mathbb{E}\left[\Re\{X(\omega')\}\right]$$
and $$\mathbb{E}\left[\Im\{X(\omega)\},\Im\{X(\omega')\}\right]-\mathbb{E}\left[\Im\{X(\omega)\}\right]\mathbb{E}\left[\Im\{X(\omega')\}\right]$$

from $$\mathbb{E}\left[\bar{x}(t),x(s)\right]-\mathbb{E}\left[\bar{x}(t)\right]\mathbb{E}\left[x(s)\right]\;?$$

I would be happy to be guided to a reference as well.

This is the most general case but the simpler case where the input signal is real is the main concern for me. I would be happy to be guided to a reference as well.

Note that in general the Fourier transform of a stationary process $x(t)$ does not exist. The Wiener-Khinchin theorem only states that under certain conditions the power spectral density of $x(t)$ exists, and it can be computed as the Fourier transform of the autocorrelation function of $x(t)$.

Having said that, if for some reason one assumes that the Fourier transform of $x(t)$ exists, then you can do the math and see if you get a useful result. So, let's see. I assume that $x(t)$ is real-valued. This is not at all necessary, but just simplifies things a bit.

Since $\text{Re}\{X(\omega)\}=\frac12[X(\omega)+X^*(\omega)]$ and $\text{Im}\{X(\omega)\}=\frac{1}{2j}[X(\omega)-X^*(\omega)]$, we can compute the desired expectations from $E\{X(\omega)X(\omega')\}$ and $E\{X^*(\omega)X(\omega')\}$. For the first of these expectations we have

$$E\{X(\omega)X(\omega')\}=E\left\{\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt
\int_{-\infty}^{\infty}x(t')e^{-j\omega' t'}dt'  \right\}=\\
=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}E\{x(t)x(t')\}e^{-j\omega t}e^{-j\omega't'}dtdt'$$

With the substution $\tau=t'-t$ and with the autocorrelation function $R_x(\tau)=E\{x(t)x(t+\tau)\}$ (assuming wide-sense stationarity of $x(t)$) we obtain

$$E\{X(\omega)X(\omega')\}=\int_{-\infty}^{\infty}R_x(\tau)e^{-j\omega'\tau}d\tau\int_{-\infty}^{\infty}e^{-j(\omega+\omega')t}dt$$

Interpreting the second integral as a distribution, and with the power spectrum $S_x(\omega)=\mathcal{F}\{R_x(\tau)\}$ we finally get

$$E\{X(\omega)X(\omega')\}=2\pi\delta(\omega+\omega')S_x(\omega')=
2\pi\delta(\omega+\omega')S_x(-\omega)=2\pi\delta(\omega+\omega')S_x(\omega)\tag{1}$$

because $S_x(\omega)$ is an even function (we assumed $x(t)$ to be real-valued). The other expectation $E\{X^*(\omega)X(\omega')\}$ can be derived in a completely analogous manner. The result is

$$E\{X^*(\omega)X(\omega')\}=2\pi\delta(\omega-\omega')S_x(\omega)\tag{2}$$

From (1) and (2) we can obtain the desired expectations

$$E\{\text{Re}\{X(\omega)\}\text{Re}\{X(\omega')\}\}=
\pi S_x(\omega)[\delta(\omega-\omega')+\delta(\omega+\omega')]\\
E\{\text{Im}\{X(\omega)\}\text{Im}\{X(\omega')\}\}=
\pi S_x(\omega)[\delta(\omega-\omega')-\delta(\omega+\omega')]$$

It remains to compute the quantities $E\{\text{Re}\{X(\omega)\}$ and $E\{\text{Im}\{X(\omega)\}$. We can easily derive them from the expectations $E\{X(\omega)\}$ and $E\{X^*(\omega)\}$:

$$E\{X(\omega)\}=E\left\{ \int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt \right\}=
\int_{-\infty}^{\infty}E\{x(t)\}e^{-j\omega t}dt=
\mu_x\int_{-\infty}^{\infty}e^{-j\omega t}dt=2\pi\mu_x\delta(\omega)$$

with $\mu_x=E\{x(t)\}$, where we have again assumed wide-sense stationarity of $x(t)$. Obviously, we get the same result for $E\{X^*(\omega)\}$:

$$E\{X^*(\omega)\}=2\pi\mu_x\delta(\omega)$$

This results in

$$E\{\text{Re}\{X(\omega)\}\}=2\pi\mu_x\delta(\omega)\quad\text{and}\quad
E\{\text{Im}\{X(\omega)\}\}=0$$

So we have

$$E\{\text{Re}\{X(\omega)\}\}\cdot E\{\text{Re}\{X(\omega')\}\}=4\pi^2\mu_x^2\delta(\omega)\delta(\omega')\\
E\{\text{Im}\{X(\omega)\}\}\cdot E\{\text{Im}\{X(\omega')\}\}=0\tag{3}$$

Combining (3) with (1) and (2) you obtain the desired result. From the result you can see that the covariances of the real parts of $X(\omega)$ and $X(\omega')$ vanish everywhere except for $\omega=\pm\omega'$. The same is true for the covariances of the imaginary parts. The difference is an additional term for $\omega=\omega'=0$ for the covariance of the real parts.


Why is Gaussian noise called so? Can you please explain: why is a specific type of noise called "Gaussian noise"? Why is it relevant to call it Gaussian? Please, explain in layman's terms.

Noise is random, but like most random phenomena, it follows a certain pattern. Different patterns are given different names.

Consider rolling a die. This is clearly random. Roll the die 1000 times, keeping track of each result. Then, calculate the histogram of the result; you'll find that you got each of 1, 2, 3, 4, 5 and 6 approximately the same number of times. This pattern is called "uniform", and throwing a die can be modeled by a "uniform random variable".

The same experiment can be repeated with thermal noise. Heat up a resistor, amplify the resulting voltage and measure its instantaneous power multiple times. Then calculate the histogram. This time, you won't find a uniform histogram; it will be shaped like a bell curve, with values near zero more common than values far away from zero. This kind of histogram is called Gaussian, after K. F. Gauss.

Gaussian random phenomena are very common in nature. It turns out that whenever the random thing you observe is the aggregate of many, independent random events, the overall random variable is Gaussian (this is technically called the central limit theorem). In the case of thermal noise, you're measuring the aggregate effect of millions or billions of randomly oscillating electrons, excited by the heat.

There's an easier way to create Gaussian randomness at home (or simulated in a computer): take many dice, say 100, throw them all many times, and keep track of the total sum of each throw. If you find the histogram again, you'll see it follows a bell curve. The reason is intuitively easy to grasp: with 100 dice, you're very unlikely to get a total of 100 (all dice would have to land in 1), but it's very easy to get a number around 350, because many different combinations add up to such a number.

To summarize, there are many different kinds of noise that can affect a signal or an image, each with different statistical properties. Gaussian noise is a particularly important kind of noise because it is very prevalent. It is characterized by a histogram (more precisely, a probability density function) that follows the bell curve (or Gaussian function). As you study it more, you'll find that it also has several other important statistical properties.

See also:

why Gaussian noise is usually used?
what's gaussian noise look like?

and many other related questions on this website.


Sense of zeropadding in a time domain I have the task related to Radon transform which contains a subtask which uses resampling by means of DFT.

Let's consider the non-periodical discretized signal (Fig.1) (for example the string of pixels) having 515 pixels length. In my implementation for resampling contains following steps:

Cyclic left shift (Fig.2).
Add zeros to the center in order to the length of signal became 2^n (in our case 1024-515 = 509 zeros we must have add) (Fig. 3).
Get DFT from this signal (Fig. 4).
Cyclic right shift. (for shifting low frequencies to the center) (Fig.5)

Fig.1

Fig.2

Fig.3

Fig.4

Fig.5

The main question:

Why we must perform cyclic shift of the signal and add zeros exactly in the center? (I assumed what this made the signal periodic)
Zeropadding makes interpolation DFT spectrum, is it correct? (I asked and someone says what it is not quite so)
Maybe someone can explain in simple way what happens with signal after zeropadding.

I have made some experiments in a Matlab and found out that any other sequence of actions can not give required result.

Now let's consider two cases:

a) (THIS CORRECT VARIANT) We has the non-periodical discretized signal (for example the string of pixels) which will be cyclic shifted to left and filled zeros in the center after that will be obtained DFT from this and to shift it back.

b) We has the non-periodical discretized signal (for example the set strings of pixels) which will be filled zeros from left and right  after that will be obtained DFT from this.

What the difference these DFT spectrums?

I have read some books but not found the answer of this zeropadding's case. It seems this can be found only by own experience.

Answer in book:

A. C. Kak and Malcolm Slaney, Principles of Computerized Tomographic Imaging, Society of Industrial and Applied Mathematics, 2001 on page 25

Shifting the data points (fftshift) and zero-padding the exact center of the FFT aperture has the property that all the even (symmetric) components with reference to the center of the original data set end up in the real part of the complex FFT result, and all the odd components end up in the imaginary part. e.g. the evenness to oddness ratio is preserved, which allows phase (with reference to the center of the window) to be easily interpolated.

Being able to interpolate phase is important in the case of a zero-padded FFT because the zero-padding interpolates the spectral magnitude as well.  Since zero-padding results in an interpolated FFT result, that means that any original non-interpolated FFT result points will need to be interpolated from the zero-padded result.  Without an FFT shift of both the non-padded and zero-padded data, these interpolated results will be different from the non-interpolated FFT results (in phase).

This technique is a simple (homework or quiz level) result of the FT property that a shift in one domain is a complex frequency modulation in the other domain.


Log energy calculation I have some code here that I don't quite understand.
It calculates the energy for a given frame the following way:

En=log(1+energy(32767*frame))

The energy function is defined as

energy = sum(s.*s)/max(size(s))

Where 's' is te given frame.

What I don't understand is why the frame needs to be multiplied by 32767, which is the max value for unsigned int 16. I couldn't find the theory behind this. Does anyone of you know why this multiplication is taking place? Does it have to do anything with the log calculation?max value for unsigned int 16. I couldn't find the theory behind this. Does anyone of you know why this multiplication is taking place? Does it have to do anything with the log calculation?

First of all, the amplitude of a signal is expressed in arbitrary units, so there is no "right" answer. It's unlikely you'll be able to compute an absolute acoustic energy value - this depends on the gains of the recording chain, and the characteristics of the transducers.

Some audio analysis/feature extraction software consider that the full scale of the input audio signal is [-32768, 32767] (due to the predominance of the 16-bit format), some other implementation expect their input to be in the [-1.0, 1.0] range. Matlab's wavread function can return both formats when reading data, with the 'double' or 'native' argument.

I suspect the multiplication by 32767 in your code is simply there to keep the output consistent with a reference implementation of the algorithm which expected data to be in the [-32768, 32767] range, or maybe to keep the results consistent with a previous version of the software that used an input signal with the 'native' instead of 'double' scale.


instantaneous frequency for a sum In the literature, the instantaneous frequency for a signal $s(t)=e^{-j \phi(t)}$ is given by: $$f(t)=\tfrac{1}{2\pi}\frac{d\phi(t)}{dt}$$

This my problem: I have a signal $$u(t)=e^{-j \alpha_1 \phi(t)}+e^{-j \alpha_2 \phi(t)}+e^{-j \alpha_3 \phi(t)}$$ with $\alpha_1,\alpha_2,\alpha_3 \in \mathbb{Z} > 0$. How to find the instantaneous frequency for the signal $u(t)$?

Eventhough Laurent has given a broader sense of the answer, let me put here the communications theory sense ot it.

The concept of instantaneous freqency emerges when you consider Frequency Modulation or Phase Modulation systems, where the message is embedded into the change of the frequency or phase of a carrier signal.

This carrier is typically a single cosine wave or a complex exponential equivalent which posses a single well defined angle, whose time derivative gives you the instantaneous frequeny; hence the message. That's why they are also called angle modulation or exponential modulation systems.

The understanding is that there must exist a single carrier to represent the modulated message. Your example includes a general sum of three complex exponentials (unless you speciy those $\alpha$ pars) whose simplification into a single carrier seems not possible on the surface. Then there won't be a unique single instantaneous frequency definition for $u(t)$.

But you can consider three different instantaneous frequencies per complex exponential for your signal $u(t)$. Indeed this would be the case when multiple FM modulated signals (at different frequencies) were existing simultaneously as they are transmitted at different channels.


Guitar pitch detection with autocorrelation This link provides code for an autocorrelation-based pitch detection algorithm but says:

  Cons: Not as accurate, doesn't work for inharmonic things like musical instruments, this implementation has trouble with finding the true peak

Why is that? What implications will it have if I want to create a guitar pitch detector (let's assume that I am only working with monophonic audio for now). Is it the noise that creates those implications? I thought that musical instruments are harmonic things, not inharmonic.

Is the YIN algorithm an improvement in this aspect? If not, what does it improve?

Cons: Not as accurate

This is just compared to the other methods.  I was measuring frequency very accurately to look for clock drift, etc: 1000.000004 Hz for 1000 Hz, for instance.  For guitar pitch detection it will be fine.

  doesn't work for inharmonic things like musical instruments

I should have said "it can't find an accurate fundamental if there is inharmonicity".  It will be shifted up because the plucked string partials are slightly sharp.

In your case, you don't want to find the fundamental, though, you want to find the perceived pitch, which is what autocorrelation is good at.  (though I don't know how accurate, and if it needs weighting, etc.)

http://www.utdallas.edu/~assmann/hcs6367/yost09.pdf

https://ccrma.stanford.edu/~dattorro/Humans.pdf

  this implementation has trouble with finding the true peak

because autocorrelation will produce peaks at the correct frequency but also at all its sub-harmonics, and this function just naively chooses the highest peak.  You need to plot the correlation curve and find a way to choose the correct results for your waveforms.


Derive DTFT of $x[2n]$ If the DTFT of discrete sequence $x[n]$ is $X(e^{j\omega})$, what is the DTFT of $g[n] = x[2n]$?

I see the textbook answer is

\begin{align*}
  G(e^{j\omega}) &= \frac{1}{2} \left( X(e^{j\omega/2}) + X(e^{j(\omega-2\pi)/2}) \right)
\end{align*}

My start of the problem:

\begin{align*}
  X(e^{j\omega}) &= \sum\limits_{n=-\infty}^\infty x[n] e^{-j\omega n} \\
  G(e^{j\omega}) &= \sum\limits_{n=-\infty}^\infty g[n] e^{-j\omega n} \\
  G(e^{j\omega}) &= \sum\limits_{n=-\infty}^\infty x[2n] e^{-j\omega n} \\
\end{align*}

How do I derive the given textbook answer?

UPDATE:

I try the obvious $m=2n$ substitution, however I don't see how you can change the variable on the summation so that it counts by 2?

\begin{align*}
  G(e^{j\omega}) &= \sum\limits_{n=-\infty}^\infty x[m] e^{-j\omega m/2} \\
\end{align*}

HINT:

If you have a sequence

$$\hat{x}[n]=\begin{cases}x[n],&n\text{ even}\\0,&n\text{ odd}\end{cases}$$

then the DTFT of $x[2n]$ can be written as

$$\sum_{n=-\infty}^{\infty}x[2n]e^{-jn\omega}=\sum_{n=-\infty}^{\infty}\hat{x}[n]e^{-jn\omega/2}=\hat{X}(e^{j\omega/2})\tag{1}$$

You can obtain $\hat{x}[n]$ as

$$\hat{x}[n]=\frac12\left(x[n]+(-1)^nx[n]\right)\tag{2}$$

Taking the DTFT of $(2)$ to express $\hat{X}(e^{j\omega})$ in terms of $X(e^{j\omega})$, and using $(1)$ will give the desired result.


Determining whether or not a system is BIBO stable I'm trying to determine whether or not a system with impulse response

$$
h(t) = \sum_{n=-\infty}^{\infty} \delta(t-2n)
$$

is BIBO stable. I haven't touched this material for a very long time -- could anyone lend a helping hand? I recall needing to show that

$$
\int_{-\infty}^{\infty} |h(t)| \, dt < \infty,
$$

but $h(t)$ is defined in a very weird way, and calculating the integral is proving to be difficult.

EDIT: Does the following argument seem reasonable? Since were summing an infinite number of $\delta$ functions, each will have an argument of $0$ for some unique value of $t$. Thus, there are infinitely many values of $t$ for which $\delta(t-2n)=1$, so the integral diverges.

In your case

$$\int_{-\infty}^{\infty}|h(t)|dt=\int_{-\infty}^{\infty}h(t)dt=
\int_{-\infty}^{\infty}\sum_{n=-\infty}^{\infty}\delta(t-2n)dt=
\sum_{n=-\infty}^{\infty}\int_{-\infty}^{\infty}\delta(t-2n)dt$$

And since

$$\int_{-\infty}^{\infty}\delta(t-2n)dt=1$$

(as correctly pointed out by Dilip Sarwate in his comment to Phonon's answer), the above sum is

$$\sum_{n=-\infty}^{\infty}1=1+1+1+\ldots=\infty$$

I believe that this is basically the argument you gave in your question.


What would produce this "hamburger" corruption on images? I have an image that has gone through an unknown number of jpeg and webp conversions (at least one but more likely 2+). Currently this image is a JPEG and contains corruptions in single color or mild gradient backgrounds. What could cause this?

JPEG projects $8\times 8$ blocks of images onto $64$ 2D cosine patterns:

The one in column $1$ and row $5$, once quantized, may look like your hamburger. Luminance and chroma components may get different subsampling patterns. I suspect that the low varying background is nearly horizontal, and due to the different processing steps, it ends up with a mid-frequency non-zero component. My main issue is its $10\times 10$ apparent size. Apparently, WebP also does some prediction on blocks.

You might be interested in compression forencics, that try to recover the history of compression. For JPEG, JPEG CHEst (JPEG Compression History Estimation for Color Images) does the job, based on paper JPEG Compression History Estimation
for Color Images.


Lock-In Amplifier Sampling Rate For a digital lock-in amplifier (synchronous demodulator) I have a sinusoidal reference signal in the range of 10-15kHz, the amplitude and phase I want to observe changes very slowly <10Hz.

My question is do I gain any benefit by oversampling the input signal, or is twice the reference frequency enough, i.e ~30kHz in my case?

I ask because the A/D converted signal is mixed with the reference oscillator and low-passed filtered to almost DC.

Oversampling here would be extremely easy, so do it!

The benefits are of a noise nature (which is probably why you're doing the whole lock-in amplifier dance in the first place); you get a part of your oversampling factor as SNR gain.


Polyphase Analysis Filter Bank I have a question regarding the polyphase DFT filter bank implementation in this page. (https://cnx.org/contents/Peqc-TK2@16/Uniformally-Modulated-DFT-Filterbank)

In figure 2, to analyze the kth filter bank branch, the signal $x(n)$ was multiplied with $e^\dfrac{j2\pi kn}{M}$. But to downconvert a signal to baseband and low pass filter it, the signal should be multiplied with $e^ -\dfrac{j2\pi kn}{M}$. I made this change and continued derivation, I ended with getting IDFT block instead of getting DFT block as shown in Figure 7.

I referred multiple books and papers for this. Everywhere I found a DFT block in the block diagram. But, I couldn't find where I am doing wrong. Please help me in knowing the correct block diagram.

With every polyphase filter bank I have worked with, the first block in the analysis phase is an IFFT, and the block in the synthesis phase is a DFT. These operations essentially cancel one other, so it should be fairly intuitive. Think of the DFTs acting as sinusoids modulators in this case rather than an operation to convert to the frequency domain. By having complementary DFTs, the modulation is essentially performed in analysis and then removed in synthesis. From looking at the figure in the text you linked, it appears that they’re may be a misprint, at least in comparison to my industry experience implementing this types of filter banks.

The following is a great, compact and concise resource on polyphase filter banks, multi-rate identities, and DFT based filter banks:
https://apps.dtic.mil/dtic/tr/fulltext/u2/a457390.pdf


Power Spectrum Estimation of three sinusoids in white noise Let's assume we have a random process consisting of three sinusoids in white noise:
$$x[n] = 3 \cdot \sin(ω_1 \cdot n + ϕ_1) + 5 \cdot \cos(ω_2 \cdot n + ϕ_2) + 2 \cdot \sin(ω_3 \cdot n + ϕ_3) + v[n]$$
where $ϕ_1$, $ϕ_2$, $ϕ_3$ are uncorrelated uniformly distributed random variables in $[0,2\pi]$

I want to estimate the power spectral density of $x[n]$. For that purpose I apply two different methods:

Periodiagram of $x[n]$
Welch method (Averaged Modified Periodogram) using Hamming window

I generated the following code and ran in Matlab:

N = 1024;
w1 = 0.3 * pi;
w2 = 0.6 * pi;
w3 = 0.7 * pi;
phi1 = unifrnd(0, 2*pi);
phi2 = unifrnd(0, 2*pi);
phi3 = unifrnd(0, 2*pi);
u = randn(1, N);
n = 0:N-1;
x = 3 * sin(w1*n + phi1) + 5*cos(w2*n + phi2) + 2*sin(w3*n + phi3) + u;

fig = figure;
plot(n,x);
title('signal x[n]');
xlabel('samples');
ylabel ('Magnitude');

fig=figure;
[p,f]=periodogram(x);
plot(f,p);

fig=figure
[pw,f]= pwelch(x, N/8, N/16);
plot(f,pw);

Which is the difference in these two methods about the resolution and the variance of the power spectrum ?

The standard Periodogram uses all the data and computes its power spectrum spectrum estimatation at once; therefore it provides more spectral resolution but less estimation reliability (larger estimation variance).

The Welch's modification to the standard Periodogram is the concept of dividing the signal into shorter blocks and averaging the computed per block Periodograms, and therefore it has less spectral resolution (due to shorter blocks) but provides smaller estimation variance, due to averaging.


Understand diagram of a process I've found this diagram in a slide, but i have no backround to properly understand it. Could you explain to me, in simple terms, how this process behaves and what's the purpose of the blocks $z^{-1}$ ?

That's a simple process, where $y(k)=u(k-3)$, as the $z^{-1}$ blocks just perform a delay by one sample. $x_1(k)$, $x_2(k)$ and $x_3(k)$ are just arbitrary names. In other words: output will be exactly like input, only delayed by three samples.


Explanation relevant to the filter response ($\mathrm{sinc}$ interpolation) - Equation $(2)$ in my question At present, I am learning the theory of operation of resampling for bandlimited periodic discrete signals using $\mathrm{sinc}$ interpolation. I am developing a design flow and having difficulty in understanding a particular expression/equation, hence I am stuck! Here I go:

Suppose we have a samples $x(nT_s)$ of a continuous time signal $x(t)$, where $t$ is time in seconds and $n$ is sample points that ranges over integers and ofcourse $T_s$ is the sampling period. $x(t)$ is bandlimited to $\pm F_s/2$, where $F_s = 1/T_s$.

Now, we know that the original signal $x(t)$ can be reconstructed or synthesized by the samples of $x(nT_s)$ using the equation:

$$
x'(t) = x(nT_s)+H_s(t-nT_s) = x(t), \quad \text{for}\quad n =\ldots,-3,-2,-1,0,1,2,3,\ldots\tag1
$$

and

$$
H_s(t) = \mathrm{sinc}(\pi F_s t)
$$

To resample $x(t)$ for the new rate i.e. $F'_s = 1/T'_s$, we have to evaluate equation $(1)$ at integer multiples of $T'_s$ i.e. $nT'_s$.

We know that, when the new sampling rate $F'_s$ is less than the original    sampling rate or frequency $F_s$, the lowpass cutoff must be placed below half the new lower sampling rate.

Doubt - Based on the highlighted statement made above, the ideal low pass is given by:

$$
H_s(t) = \min\left\{1,F'_s/F_s\right\}\mathrm{sinc}\left(\pi\min\left\{F'_s,F_s\right\}t\right) \tag2
$$

where the scale factor maintains unity gain in the passband.

What does the equation $(2)$ represent? I understand that, if $F'_s$ is lower than $F_s$, then the filters cut off frequency should be half the $F'_s$, however, I am having trouble understanding this equation $(2)$ and how to benefit from it in the design process. Any explanation regarding this would be appreciated.

Source/Reference:

$\scriptstyle{\textrm{Julius O.Smith & Phil Gosett "A Flexible Sampling Rate Conversion Method" IEEE, 1984}}$.

Eq$(1)$ doesn't look right. It should be:

$$x'(t) = \sum_{n=-\infty}^{+\infty} x(nT_s) \mathrm{sinc}(t-nT_s)$$

For your doubt, consider subsampling (lowering $F_s$) as equivalent to reconstructing $x(t)$ and then sampling it at a lower frequency.

Eq$(2)$ comes from the combination of $2$ lowpass filters:

The ideal reconstruction filter (from eq.$1$)
a new anti-alias filter previos to sampling at the lower $F_s$

These $2$ filters combined should result in eq$(2)$ (several lowpass filters in cascade is equivalent to one filter with the minimum bandwidth, and the product of the gains).


Detecting Individual Photos on a Scanned Page I'm working with a collection of old photographs. They were scanned on brown paper in various configurations. Some scans have a single photo, others have three or four. Some of the photos are mounted on white cards, others are not:

I developed an algorithm to find the individual photographs on each scan, which is described here. The gist is to:

Find the median color (independently for each RGB channel, 0-255)
Blur the image.
Binarize according to whether each pixel's RGB values have an RMSE of <20 from the median.

This leaves a mostly-black image with a few large white rectangles on it, like so:

To find the rectangle coordinates, I use the following procedure:

Pick a random white pixel, (x, y) (statistically, this is likely to be in a photo)
Call this a 1×1 rectangle.
Extend the rectangle out in all directions, so long as you keep adding new white pixels.
If this rectangle is larger than 100×100, record it as a photo. (The originals are much larger than the photos I've uploaded for this question.)
Color the rectangle black.
If <90% of the image is black, go back to step 1.

This has worked on most (90+%) of the scans that I've tested it on, but it is quite ad-hoc, particularly the RMSE 20 threshold in the first portion and the entire procedure to find the rectangle coordinates.

Is there a more principled way to find these photo rectangles? How would someone with experience in image processing solve this problem?

Most of the image processing/computer vision algorithms are ad hoc. Your algorithm is good. However, I think that a few points can be improved:

Instead of finding median color, you could do use the geometry of your images, which is more robust. For example, you know that the background part (The one that touches the boundaries is black, and the one inside it has the color you are looking for.
In order to find the rectangle coordinates, you should use a blob analysis procedure. In Matlab it is done by using regionprops command. It will give you properties like

area of blob (To remove the small ones)
solidity (Can be used to separate the rectangles from non-rectangles)
coordinates of pixels

Also, you could try to use morphological pre-processing like opening, to remove small blobs.


Stability of open-loop transfer function from its Nyquist plot I am facing a confusion on understanding system open-loop transfer function stability from its Nyquist plot.

According to the formula, for open loop transfer function stability: $$Z=N+P=0$$
where $N$ is the number of encirclements of $(0,0)$ by the Nyquist plot in clockwise direction.

Now how can I verify this formula for the open-loop transfer function: $$H(s)=-\frac{1}{s^3(s+1)}$$
The Nyquist plot is attached in the image.

Here $N=1$. Now how can I find $Z$ and $P$ and verify the formula?

As mentioned by Tendero in a comment, you should be interested in the stability of the closed-loop transfer function

$$Q(s)=\frac{H(s)}{1+H(s)}=\frac{s^3(s+1)}{s^4+s^3-1}\tag{1}$$

where I've assumed unity feedback. Since the stability of $(1)$ is determined by the zeros of the denominator $1+H(s)$, we're interested in the encirclement of the point $-1+j0$ (and not in the encirclement of the origin). In your example, the number of clockwise encirclements of the point $-1+j0$ happens to equal the number of encirclements of the origin, i.e., we have $N=1$. Since the given $H(s)$ has no poles in the right half-plane, the number of zeros in the right half-plane equals the number of clockwise encirclements, i.e., $Z=N=1$, which means that $Q(s)$ has one pole in the right half-plane, and, consequently, the closed-loop transfer function is unstable.

This is also easily verified by computing the poles of $Q(s)$, which are the roots of the polynomial $s^4+s^3-1$:

roots([1,1,0,0,-1])
ans =

  -1.38028 + 0.00000i
  -0.21945 + 0.91447i
  -0.21945 - 0.91447i
   0.81917 + 0.00000i

This shows that one of the four roots is in the right half-plane, as predicted by the Nyquist plot.