What information does fourier transform carry? As one starts learning signal processing, then comes inevitably the topic of Fourier Transforms. Unfortunately I have difficulties not in computing but in interpreting the results of the Fourier Transforms, in particular the one being Continuous-Time Fourier Transform, CTFT, of the signal $x(t)$ which is: $$X(j\omega) = \int_{-\infty}^{\infty}{x(t)e^{-j\omega t}dt}$$

Now I wonder what kind of information does this $X(j\omega)$ give about the signal $x(t)$?  An example is highly appreciated, if possible.

There are a variety of Fourier Transforms (and Series) such as Continuous Time FT, Discrete Time FT, Discrete FT all of which are generally attributed to the fundamental assertion made by J.B.Fourier about at the beginning of 19th century, which claims (without proof) that "if a continuous-time signal (function) x(t) is periodic with T, then it is possible to represent that signal x(t) as an infinite sum of harmonically related trigonometric functions (sines and cosines) as $ x(t)= \sum{ a_k \sin { 2\pi kt \over T} + b_k \cos { 2\pi kt \over T}}$  in which the weights $a_k$ and $b_k$ (the Fourier Coefficients) represent the amount of that particular harmonic in the signal x(t) being analysed"

In fact the above argumentation is strictly for Continuous-Time Fourier Series. But the core idea is generalised into Fourier-Transforms of aperiodic and periodic (with the help of impulse $\delta (t)$ functions) signals. There are conditions on which signals can have such a representation.

In essence, computing a Fourier Transform means finding those coefficients $a_k$ and $b_k$ for which the method is suggested by the analysis equation of the Fourier transform, while the equality in the first paragraph is noted as the synthesis equation.

Eventhough it is quite solid to understand the meaning of those sinusoids inside a periodic signal, when it comes to non-periodic signals, for which we use the Fourier Transforms, the exact meaning of what a single sine wave represents inside such a signal is a little vague and instead we emphasize the transient character of that signal under concern and the necessity of existance of a continuum of infinetely many sine waves.


Does converting an FIR filter from causal to non-causal remove phase delay? Let's say I have a causal 3-tap FIR filter:
$$y[n] = ax[n] + bx[n-1] + cx[n-2]$$
Assuming I have access to all the data prior to filtering and generate a similar non-causal filter:
$$y[n] = ax[n+1] + bx[n] + cx[n-1]$$
Would this new filter have the same frequency response but with different phase delay?

Correct. The magnitude response would be the same and the phase response would be summed with $\omega$ (omega).


what is $1/f$ property of the visual world? I was watching a talk by Edward H. Adelson on the perception of materials and surfaces. He mentions a $1/f$ property of the visual world. He says:

  low frequencies are prevalent than high frequencies.

Also,

  the simplest model of the world is that it is random noise with a $1/f$ spectrum.

I didn't understand what he meant. Anyone willing to explain the concept he is talking about?

"Essentially, all models are wrong, but some are useful" (attributed to G. Box)

One of the simplest morphological model for 2D images decomposes objects into structure and texture. Both can be used to detect or understand objects. Even though image formation models are still complicated, objects often appear as smooth intensity gradient, filled with decoration, bounded by edges (some being false).
Smooth intensity gradients are often the most energetic, thus low frequencies can be prominent. Edges, though energetic, are sparser and local, thus possibly less present in a spectrum.

So, a simple "low frequencies are prevalent than high frequencies" model is a decreasing spectrum. Now what would be the decreasing simplest model in terms of elementary functions? Possibly $1/f$. And yes, $\exp^{- |f`|}$  could seem better, especially with the singularity at $0$. More involved models involve $1/f^p$ or $\exp^{- p |f|}$ shapes. Indeed, such models are dependent of the transformations used to reshape the data, and Fourier is a very simple reshaping tool. More local tools like wavelets use Generalized Gaussian Distributions: $\exp^{-\alpha x^\beta}$ in the time/space domain.

Simpler can also be seen from a computational point of view, see Frequency Synthesis of Landscapes (and clouds):

  Frequency synthesis is based upon the observation that many "natural"
  forms and signals have a 1/fp frequency spectra, that is, their
  spectra falls off as the inverse of some power of the frequency where
  the power is related to the fractal dimension.

  This leads naturally to a method of generating such fractals:

  Generate a random (white noise) signal.
  Transform this into the frequency domain.
  Scale the resulting spectra by the desired $1 / f^p$ function.
  Inverse transform.

Such simple models can reproduce realistic shapes. From the frequency domain:

one gets:


Efficiently calculating autocorrelation using FFTs I'm trying to calculate an autocorrelation on a platform where the only accelerated primitive I have available is the (I)FFT. I'm having a problem though.

I prototyped it in MATLAB. I am, however, slightly confused. I assumed that it works simply as follows (this is from memory so apologies if I've got it slightly wrong).

 autocorr = ifft( complex( abs( fft( inputData ) ), 0 ) )

However I get a different result than I get from using the xcorr function.  Now I'm fully expecting not to get the left hand side of the auto correlation (as it's a reflection of the right hand side and thus not needed anyway). However, the problem is my right hand side appears to be, itself, reflected around the halfway point. Which effectively means that I get about half the amount of data that I'm expecting.

So I'm sure I must be doing something very simple wrong, but I just can't figure out what.

pichenettes is right, of course. The FFT implements a circular convolution while the xcorr() is based on a linear convolution. In addition you need to square the absolute value in the frequency domain as well. Here is a code snippet that handles all the zero padding, shifting & truncating.

%% Cross correlation through a FFT
n = 1024;
x = randn(n,1);
% cross correlation reference
xref = xcorr(x,x);

%FFT method based on zero padding
fx = fft([x; zeros(n,1)]); % zero pad and FFT
x2 = ifft(fx.*conj(fx)); % abs()^2 and IFFT
% circulate to get the peak in the middle and drop one
% excess zero to get to 2*n-1 samples
x2 = [x2(n+2:end); x2(1:n)];
% calculate the error
d = x2-xref; % difference, this is actually zero
fprintf('Max error = %6.2f\n',max(abs(d)));


phase, magnitude, audio file and information I've been doing some thinking lately. I know that for image files a larger portion of information is contained in the phase of the signal. How does that go for the audio? Does the same hold or is the amplitude the one with the more information content? Thx all.

As the comments said. "monoaural phase" is not very important and monoaural phase changes are often inaudible. "binaural" phase is very important for localisation and imaging, particularly in stereo files.

"monoaural" means the phase of a single channel or that the same phase shift or modification is applied to all channels the same way.

"binaural" means specifically the difference in the signal at the left and the right ear.


"Revolving" QPSK Constellation I am trying to code a QPSK transmitter in VHDL but first I want to create one out of Simulink blocks and translate it into VHDL code. I am using an example project for the PlutoSDR as a guide.

This is the block diagram:

I ran into the following error for which i used Simulink's "fix" option and got this project to compile:

When i run this project this is the output i am seeing on Constellation Diagram 1 & 2:

What's happening is the points on the diagrams are continuously moving in a "circular" way. I've uploaded two GIFs of the constellation diagrams 1&2 from my Simulink project.

Constellation Diagram 1:

Constellation Diagram 2:

The output of Constellation Diagram, however, does not move and looks like this:

This is the output of the example on the Mathworks website (guide) seen as "Constellation Diagram" in the block diagram above.

My questions: Why did i have to put "Single Task Rate Transmission" boxes into my project as suggested by Simulink? Why do Constellation Diagram1 and Constellation Diagram2 "dance" around like this? What am i missing?

Thank You!

You are seeing the "Digital IF” carrier assuming you are modelling what you want to implement digitally in VHDL. The sine and cosine blocks together with the multipliers and combiner perform a frequency translation from DC to the carrier mentioned.  Thus the constellation at this point is really not observable as it will be "spinning" at this carrier rate. (The constellation as observed in the upper figure is at DC and to note it appears to be sampled at 2 samples per symbol, which is typical, so it is really every other sample (the ones on the diagonals) that represents the QPSK constellation at the correct symbol sampling locations).

To observe the constellation you would need to frequency shift the signal back to DC while being careful not to introduce phase and frequency offsets, or shifts in sample timing (a receiver does this with carrier and timing recovery loops so that it can run completely on an asynchronous clock from the transmitter).

(To note, an optional approach is "Zero-IF" where you implement the frequency translation in the analog domain using an IQ mixer (which is the same functional block with analog components)


Hilbert transform of sinusoid -- apparent contradiction The arrow "$\xrightarrow{\mathcal H}$" denotes a Hilbert transform:

  $$\cos(\omega t)\xrightarrow{\mathcal H} \sin(\omega t)$$

happiness.

But

  $$\cos(\omega t) = \cos(-\omega t) \xrightarrow{\mathcal H} \sin(-\omega t) = -\sin(\omega t)$$

sadness.

This is a contradiction. What have I done wrong?

The error lies in the assumption that if $g(t)$ is the Hilbert transform of $f(t)$, then the Hilbert transform of $f(-t)$ must be $g(-t)$. This is not the case.

Let $f^-(t)=f(-t)$. Then we have

$$g(t)=\mathcal{H}\{f\}(t)=\frac{1}{\pi}\text{p.v.}\int_{-\infty}^{\infty}\frac{f(\tau)}{t-\tau}d\tau\tag{1}$$

and

$$\begin{align}\mathcal{H}\{f^-\}(t)&=\frac{1}{\pi}\text{p.v.}\int_{-\infty}^{\infty}\frac{f(-\tau)}{t-\tau}d\tau\\&=\frac{1}{\pi}\text{p.v.}\int_{-\infty}^{\infty}\frac{f(\tau)}{t+\tau}d\tau\\&=-\mathcal{H}\{f\}(-t)\\&=-g(-t)\end{align}\tag{2}$$

So if $f(t)=\cos(\omega t)$ and $g(t)=\mathcal{H}\{f\}(t)=\sin(\omega t)$, then the Hilbert transform of $f^-(t)=\cos(-\omega t)$ equals $-g(-t)=-\sin(-\omega t)=\sin(\omega t)$. Happiness again.

EDIT:

In response to your comment about substituting for $\omega$ instead of reversing time, note that

$$\mathcal{H}\{\cos(\omega_0 t)\}=\sin(\omega_0 t)$$

is only valid for $\omega_0>0$. This can most easily be seen in the frequency domain, where the Hilbert transform corresponds to a multiplication with $-j\;\text{sign}(\omega)$. The Fourier transform of $f(t)=\cos(\omega_0t)$ is given by

$$F(\omega)=\pi\left[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)\right]$$

If $\omega_0>0$, multiplication by $-j\;\text{sign}(\omega)$ gives the Fourier transform of $\sin(\omega_0t)$. However, if $\omega_0<0$, we obtain the Fourier transform of $-\sin(\omega_0t)$, because then $\delta(\omega+\omega_0)$ appears at positive frequencies instead of $\delta(\omega-\omega_0)$.


Low pass filter algorithm origin I have been doing some research, because I used this algorithm that's very spread in the web:

y += a * (x - y)

or

y = (a * y) + (x - (a * x)

Where x is input, and y is the filter output, and a is actually alpha, which there's some math behind it to define it's value.

This algortihm is seen on wikipedia, this article and many more.

But none actually shows why this algorithm is a low pass filter and from which formula it actually came from(wikipedia has no source for it). I'm doing a small academic paper, so I need to know exactly from where it came.

I tried to make my derivation match up exactly with embeddedrelated's. The filter you cited uses an approximation, which further obfuscates intuitively understanding it.

embeddedrelated does cite the frequency transfer function that the IIR filter is based off of:
$$
H(s)=\frac{1}{1+s\tau}
$$

The transfer function for a simple RC low pass filter is:
$$
H(s)_{RC}=\frac{1}{1+sRC}
$$
$$
\tau=RC
$$

Schematic:

We can derive a difference equation from the above schematic.

$$
I=C\frac{dV_{out}}{dt}
$$
$$
I = \frac{V_{in} - V_{out}}{R}
$$
$$
C\frac{dV_{out}}{dt}=\frac{V_{in} - V_{out}}{R}
$$
$$
C\frac{V_{out}[n]-V_{out}[n-1]}{ts}=\frac{V_{in}[n] - V_{out}[n]}{R}
$$
$$
V_{in}[n]+\tau f_s V_{out}[n-1]=\tau f_s V_{out}[n] + V_{out}[n]
$$
I have tested this form:
$$
\frac{V_{in}[n]+\tau f_s V_{out}[n-1]}{\tau f_s + 1}=V_{out}[n]
$$
$$
\tau f_s >> 1
$$
I have not tested this form:
$$
\frac{V_{in}[n]}{\tau f_s} + V_{out}[n-1]\approx V_{out}[n]
$$


Output of LTI (in time and frequency $\omega$ domain) : when input goes through LPF I would like to raise a mathematical question :
Let's say we are been given :
$$x(t) = \begin{cases} \cos(\pi t) & |t| \leq 0.5 \\
                                       0 & \textrm{otherwise} \end{cases}$$
If $x(t)$ goes through an LPF : $H(\omega) = \begin{cases} \pi - 0.5|\omega| & |\omega| \leq 1.5\pi\\
0 & \textrm{otherwise} \end{cases}$
Compute and sketch :

Fourier coefficients of $x(t)$, as well as $X(\omega)$
$y(t)$ and $Y(\omega)$ where $y$ is the output of the LTI. Is y(t) periodic? If so compute it's period and its trigonmetric fourier anlysis coefficients
I did the math and I've got this:
$a_n = \frac{2}{\pi}\frac{\cos(\frac{n\pi}{2})}{1-n^2}$ which can not be defined for $n=-1$ or $ n=1$
update:* We have to compute $a_1$ and $a_{-1}$ in seperate: If we do the computations we gain : $a_1= a_{-1}=0$
I then did :
$$c_m = \frac{1}{2}(a_m + jbm) \rightarrow c_m = \frac{a_m}{2}$$
$$x(t) = \sum_{m=-\infty}^{\infty}c_me^{jm \omega_o t} \rightarrow
X(\omega) =2\pi \sum_{m=-\infty}^{\infty}c_m\delta(\omega - m\omega_o)$$

What I tried:
Now we know for an LTI : $y(t) = x(t) * h(t) \rightarrow Y(\omega) = X(\omega)H(\omega)$
but since $H(\omega)$ exists only in $\delta=[-1.5\pi, 1.5\pi]$ and $\omega_o = \pi$ then we only need $X(\omega)$ in $\delta$ .
But , $X(\omega)$ is $0$  for $|m|=1 \rightarrow |\omega| =\pi$ and therefore $X(\omega) \neq 0 \rightarrow m=0$ in $\delta$ so :
$X(\omega) = 2\pi a_0 \delta(\omega)$ in $\delta$
we also know $\int_{-a}^{a}f(x)\delta(x)dx=f(0)$  so:
\begin{align}
y(t) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}Y(\omega)e^{j\omega t}d\omega \rightarrow \int_{-1.5\pi}^{1.5\pi}\frac{1}{2\pi}Y(\omega)e^{j\omega t}d\omega \\
&\rightarrow \int_{-1.5\pi}^{1.5\pi}\frac{1}{2\pi}(\pi^2a_0 -0.5a_0|\omega|)e^{j\omega t}\delta(\omega)d\omega = \frac{1}{2\pi}\pi^2 a_0 = \frac{1}{2\pi}\pi^2 \frac{2}{\pi} = 1
\end{align}
I assume that it should have been periodic but it's not. So where did my solution go wrong?

Since this appears to be homework, I can give you a few generic pointers

LTI systems maintain periodicity: if the input is periodic so is the output and vice versa.
If a signal is finite in length it can't be periodic
If a signal has finite bandwidth, it must be infinite in time.
You can model the input as $x(t)$ as cosine multiplied with a rectangular function.
Multiplication in the time domain is convolution in the frequency domain.


Scale FFT frequency range for a bars graph I'm trying to show frequency spectrum for audio data as a bars graph using FFT.

If I have data for a sine wave at 440 Hz and with sampling rate 44100 Hz and analyze its first 1024 samples I get an array with 1024 FFT values.
I know that the first half of those values (512) represents the frequency range equal to half the sampling rate (22050 Hz).

So to draw a bars graph I just take the first 15-20 values from the FFT array, calculate the frequency range for each element and plot it.

In the case of FFT size 1024 I frequencies each ~43 Hz wide.
But if I change FFT size to 2048 each frequency range is ~21.5 Hz. So the peak of the spectrum of the sine wave will be outside my 15-20 bars.

So I wanted to ask if it's possible to scale/normalize the X axis for a frequency range and a constant graph bars count?

The example limiting data would be:

frequency range 0-500 Hz
16 bars on the graph

I guess it's some kind of bins/magnitudes grouping. But should I sum them or get max? And what could be done with magnitudes when grouped bin number isn't integer? (I hope my explanation is not too confusing, it's hard to choose the correct terminology)

Thanks.

Here is my test code in Python:
(I used a static int array imported from another file that represents a sine wave)

import scipy.io.wavfile as wavfile
import scipy
import scipy.fftpack
import numpy as np
from scipy.signal.windows import *
from matplotlib import pyplot as plt

from imp.sine_440_hex import *

def run():
  # fp = 'e:/sine_440_hz.wav'
  # FS, data = wavfile.read(fp)

  # data - array of int for a sine wave at 440Hz and amplitude 0.5
  DS = len(data)
  FS = 44100

  FN = 1024
  # FN = 2048
  MAXS = 32768

  MN = 15   # number of first frequencies in the FFT result

  ndata = np.array(data)
  ndata = ndata / MAXS

  w = np.hamming(FN)
  y = ndata[:FN] * w

  freqs = scipy.fftpack.fftfreq(FN, 1/FS)
  mags = abs(scipy.fft(y))

  dbfs = 20 * np.log10(mags * 2 / np.sum(w))
  dbs = dbfs + 120

  # -------- plot --------
  fr_res = freqs[:MN]
  db_res = dbs[:MN]

  ax1 = plt.subplot(211)
  plt.plot(fr_res, db_res)
  plt.grid(True)
  plt.xlabel('Frequency [Hz]')
  plt.ylabel('Amplitude [dB]')

  ax2 = plt.subplot(212)
  plt.bar(fr_res, db_res, fr_res[1]-10)
  plt.grid(True)
  ax2.set_ylim([80, 120])
  ax2.set_xticks(fr_res)
  ax2.set_xticklabels(fr_res.astype(int))

  plt.show()

# -------
run()

And example graphs:

For FN = 1024

For FN = 2048

To achieve that, use frequencies rather than bins. To convert from bin number to the frequency it represents use:

$$f = \frac{k}{N_{FFT}} \cdot Fs$$

Where:

$k$ is a bin number whose actual frequency you want to recover
$Fs$ is the sampling frequency
$N_{FFT}$ is the number of points of the FFT

You can work out your $N_{FFT}$ to ensure that your scanned frequency range will be divided into 16 segments. For example, to cover the frequency range $f_l=300 f_h=3000$ Hz with $N_{bars}=16$ bars and $Fs=44100 Hz$, you would have to use:

$$N_{FFT} = \frac{Fs \cdot N_{bars}}{|f_h-f_l|}$$

Where basically, you work out the spacing between the two frequencies $f_l, f_h$ for 16 segments and then "translate" it to the $Fs$ (but the ratio has been simplified there and that's why $Fs \cdot N_{bars}$)

This works out to be:

$$N_{FFT} = \frac{44100 \cdot 16}{|3000-300|} \approx 261.3333 $$

Which will have to be rounded to $262$ points.

So, now you know that your FFT has 262 points with a spacing of $\frac{Fs}{N_{FFT}}=\frac{44100}{262}\approx 168.3206$ Hz.

Your first bin is at $\frac{300}{168.3206} \approx 1.7823$ which will have to be rounded to $2$ and 3000 Hz should be at bin $2 + 16 = 18$.

Let's see, $\frac{18}{262}\cdot 44100 \approx 3029.7709 $ Hz.

Hope this helps.


Testing filter code w/ Octave I want to test some code generated by this site. I selected Bessel LP 1st sample rate $600\textrm{ Hz}$ corner $8\textrm{ Hz}$ long $10\textrm{ bit}$. If I adjust the code for octave to be:

clear all;
close all;

N = 1024

v0=256;
v1=256;

m = [];

for mi = 1:128
    ip=round(rand(N,1)*512);
    for i = 1:N
        v0 = v1;
        tmp = ((((ip(i) * 2699550)/32)
              + ((v0 * 3856860)/2))
              + 1048576) / 2097152;
        v1 = tmp;
        op(i) = v0 + v1;
    endfor
    op = fft(op,N);
    m = [m; abs(op(2:N/2))];
endfor

x=[2:N/2];
y = mean(m);

semilogx(x,20*log(y)-170);

I get this plot:

which looks vaguely correct but the slope is $\approx 13\textrm{ dB/octave}$ and not the expected $6\textrm{ dB/octave}$.

I have limited experience with DSP so I need a method to verify the actual code. Can someone recommend a procedure for verifying code like this (that works)?

I think the solution is to use $\log_{10}(y)$ instead of $\log(y)$.  log(y) in octave is the natural $\log$, but you want a log base $10$, which is log10(y).

Confirming this:

$20\log_{10}(0.5) = -6 \textrm{ dB/octave}$

$20\ln(0.5) = \frac{-13.86}{20} \textrm{ nepers/octave}$

I am using $\textrm{ln}$ for natural log to avoid any confusion.


What is a complex helical sequence? Matlab says

  x = hilbert(xr) returns a complex helical sequence, sometimes called the analytic signal, from a real data sequence.

How did the name complex helical come about? How is it related to analytic signal?

The description complex helical sequence indicates that it is a sequence of complex numbers that, when plotted in three dimensions (real part, imaginary part, position in the sequence), resembles a screw, or helix.

See this illustration (by RobHar, on Wikimedia Commons):

From a mathematical point of view, all signals can be decomposed into complex sinusoids (sine waves), which are terms of the form $A(x) = e^{(\omega x +\phi)i}$. This is an often useful perspective. As $x$ is varied, $A(x)$ describes a circle in the complex plane; plotted in three dimensions, it is the helix from above.

Real-life signals normally consist of real numbers only, i.e. the imaginary part is always zero. This means that the complex sinusoids of which it is composed come in mirror pairs, where one runs clockwise and one anticlockwise; otherwise (frequency, amplitude, phase), they are identical. In such a pair, the imaginary parts cancel, and the result is a real-only signal.

An analytic signal is made from a real signal by removing one half from each of those pairs so that all remaining sinusoids move in the same direction. It is no longer real-valued, i.e. the imaginary part is in general non-zero, but the real part is identical to the original signal.

(Technically, the sense of rotation is considered the sign of the frequency $\omega$, and an analytic signal only contains sinusoids with non-negative frequency.)


Contradictory results from power spectral densities and autocorrelation in shot noise I have a white noise generator circuit which involves a Zener diode meant to undergo breakdown across a resistor $R1$, and am trying to compare power spectral densities and autocorrelation of the output I get. I recorded the data on an oscilloscope with microsecond, milisecond, and half second timescales.

First, I tried looking at what sample rates that I chose from was going to best characterise the noise. I figured a lower sampling time would be better, sort of arbitrarily, but the noise histograms I got for lower sampling rates tended to look more Poissonian in nature (well, most of all Gaussian looking):

I guess I should take it that longer sampling rates tend to characterize the shot noise better?

As I say, Zener diode breakdown noise is shot noise, so I was ideally trying to see which combination of Zener diode voltage and $R1$ I could use to generate the "whitest" looking noise. That is, that could generate a flat power spectral density and most random-looking autocorrelation (all lags except for $1$ near $0$).

However, I have some contradicting graphs when I plot them. For instance, for a Zener diode voltage of $12 V$ and $R1 = 100 k \Omega$ my PSD looks like this (using plt.psd(data,Fs=1e6))

which is somewhat flat? However, the autocorrelation certainly disagrees with this being perfectly random (the axes are the canonical ones - lags and ACF):

which, if I am to assume the PSD shows flatness, contradicts it. Meanwhile, for a milisecond time between samples with the same Zener voltage and $R1$ resistance..

which again seems to contradict itself, but now in the other direction -- the near 0 autocorrelation supports the idea that noise is random, yet the PSD is not flat.

Meanwhile, with a half-second time between samples..

my data seems to be the most self-consistent here, if I can say that the PSD is more-or-less flat. Is this essentially telling me that the proper time between samples should've been one half-second? If not, what's going on here?

I'm very shaky to the idea that the faster sampling rate failed to capture the noise distribution as well as a slower one.

Your Zener diode will have an intrinsic/paracetic capacitance (you might also have connected a separate capacitor in your circuit as well, across zener' output or nearby) and assocaited to it will be a RC time constant with the resistance across it. For the noise to show it's random nature accross this circuit's tapping point, we would have to wait for a time that is  certain order of the time greater than the time constant. That is why a slower sampling rate is able to better capture that.

That is why you observe that when the sampling rate was 1MHz, the autocorrelation was high all accross and when the sampling interval is of the order of Milliseconds, you have very much acheived a good correlation value because you are nearing the typical time constant of an RC circuit, for the half second sampling the noise fluctuations have indeed taken effect across the tappping point as now the time is much greater than the time constant and you observe true randomness, that is short noise that is white.


Sampling low pass filtered white noise If we filter out ideal white noise using an ideal LPF of cutoff frequency 10 KHz and then sample it at 30 KHz , is the resulting discrete signal statistically independent? I would like to know the statistical behaviour of the output signal.

I was attending an on-line test. My answer was told wrong.Below is the snapshot from the on-line test.

I would also like to know what would be the result if sampling frequency is below Nyquist frequency?

The output of your filter is what is sometimes called band-limited white noise. In
your particular case, the autocorrelation function of the output noise is a sinc
function whose zeroes are every $100$ microseconds, that is, samples taken at
the rate of $10^4$ samples/second are uncorrelated.  Your samples at $3\times 10^4$
samples per second are closer together and thus are correlated.  In fact,
the correlation coefficient between successive samples is
$\displaystyle \frac{\sin(\pi/3)}{\pi/3} = \frac{3\sqrt{3}}{2\pi} \approx 0.827$

Note added after question was edited: Samples taken $0.03$ milliseconds apart
are at a frequency of $33.33\ldots$ kHz, not $30$ kHz as you say in the
part of the question that you typed in yourself. Regardless, the answer (B) is
incorrect and the reasoning given in support of answer (B)
is bogus. What is being sampled is
not white noise but filtered white noise, and Answer A is correct.

For (B) to be the correct answer,
the sampling rate must be a sub-multiple of $10$ kilosamples per second,
that is the samples must be spaced $100$ microseconds apart
or an integer multiple of $100$ microseconds
apart.  Sampling every $0.3$ milliseconds, that is, every $300$ microseconds
meets this criterion; sampling every $0.03$ milliseconds does not.
Now sampling at $10$ kilosamples/second  makes the samples uncorrelated
(answer (C)) and to get from this to the stronger result that
the samples are statistically independent (answer (B)), we need
further assumptions about the noise. The standard assumption is
that the noise is Gaussian.

  In summary, (B) would be the correct answer
  if the sampling were done every $0.3$ milliseconds (since the problem
  statement already includes the assertion that the noise is Gaussian)
  and (C) would be the correct answer if it did not say that the noise
  is Gaussian since we could not the make the specialization from
  uncorrelated to independent.
  But when the sampling is done every $0.03$ milliseconds, then (A)
  is the correct answer.


Finite length truncated exponential sequence $\mathcal Z$-transform, zeros over circle explanation I'm looking at an example on how to obtain the $\mathcal Z$-transform from a finite length truncated exponential sequence, namely:
$$x[n] =
\begin{cases}
a^N &\text{for} & 0 \leq n \leq N-1\\[2ex]
 0 & \text{otherwise}
\end{cases}
$$

The $\mathcal Z$-transform from that sequence is:
$$
\frac {1}{z^{N-1}} \cdot \frac {z^{N} - a^N}{z - a}
$$

It can be seen that there are $N-1$ poles at the origin and $N-1$ zeros at $z = a$ (because the pole in $z=a$ is cancelled by a zero). My doubt is on the location of the zeros over the $z$-plane.

I know that the $N-1$ zeros at $z = a$ can be expressed as:

$z = a\cdot1\ \text{and}\ 1 = e^{j2\pi k}$ so:
$$
z_k = ae^{j2\pi k}\quad\text{with}\quad k = 1,2,\ldots,N-1
$$
Yet in the example I'm looking at, they say that the zeros are located at
$$z_k = ae^{j\frac{2\pi k}{N}}\quad\text{with}\quad k = 1,2, \ldots, N-1$$
So they are spaced at $\frac{2\pi}{N}$ instead of being all at the same place and I don't understand where that $N$ came from.

The $z$-transform is

$$\frac{1}{z^{N-1}}\cdot\frac{z^N-a^{N}}{z-a}$$
The zeros will therefore be given by
$$\begin{align}
z^N&=a^{N}\cdot1\\
&=\left(a^{N}\right)e^{j2k\pi}\\
&=\left(ae^{\frac{j2k\pi}{N}}\right)^{N},\,k=1,2,\cdots,N-1
\end{align}$$
Hence,
$$\boxed{z_k=ae^\frac{j2k\pi}{N},\,k=1,2,\cdots,N-1}$$


DFT Polyphase Filter bank What is the purpose of DFT & IDFT blocks which are placed after the Analysis filter bank & before the Synthesis filter bank.
The reason i am asking this question is because I tried to design Polyphase DFT filter bank & got the same output with & without DFT-IDFT block.

I guess you are referring to an implementation like the one shown here. If you leave out the DFT/IDFT blocks, then you don't have the modulation step which is necessary for actually implementing a filter bank, and you will not be able to do appropriate sub-band processing between the analysis and synthesis stages.

When you test your analysis and synthesis stages without sub-band processing in between you will obviously get the same result if you leave out the DFT/IDFT blocks. Analysis and synthesis together are then a pure delay and you might as well replace everything by a pure delay, but this defeats the purpose of the filter bank which is to split the input signal into sub-bands for further processing, and then to reconstruct the original signal from these separate sub-bands.


Question related to Discrete signal convolution I have two finite signals $x[n]$ and $h[n]$, I want to convolve these signals

using the method defined in Example 2.2 of Signals and Systems by Oppenheim. I am getting this answer

$$ y[-3]=4, y[-2]=11, y[-1]=12, y[0]=8, y[1]=0, y[2]=-2, y[3]=-2, y[4]=-1 $$

and when I tried to solve this question using the method defined in example 2.1 of Signals and Systems by Oppenheim, I am getting different answer. Shouldn't  it be the same?
$$ y[-5]=4, y[-4]=11, y[-3]=12, y[-2]=8, y[-1]=0, y[0]=-2, y[1]=-2, y[2]=-1 $$
its 2 unit shifted right side

Your solution based on Example 2.2 of A.Oppenheim's Signals & Systems is wrong.

In that example, the output
$$y[n] = \sum_{k=-\infty}^{\infty} {x[k]h[n-k]}$$

is computed by the graphical method (a.k.a. flip-and-drag method).

Based on the nonzero range of input $x[n]$, this becomes:
$$y[n] = \sum_{k=-2}^{2} { x[k]h[n-k] }$$

Also based on the nonzero range of $h[n]$, it can be seen that $y[n]$ has nonzero samples in the range from $n=-5$ to $n=2$.

In the flip & drag method, it's best to draw the shifted functions $h[n-k]$, to guide the computation process.

The key is to recognise the nonzero range for the flipped & dragged function $h[n-k]$ ; time reversed $h[-k]$ shifted by "$n$", drawn on an axis of $k$ and thus considered as a function of $k$.

You can see that the overlap between the functions $x[k]$ and $h[n-k]$ begins at the shift of $n = - 5$, from the left side of $x[k]$ (the sample $x[-2]$) and continues until the shift for $n = 2$ at the right side of $x[n]$ (sample $x[2]$). For any other "$n$" value, there is no overlap between $x[k]$ and shifted $h[n-k]$ and output $y[n]$ becomes zero.


Autocorrelation Computation with MATLAB I have done the computation of the autocorrelation of a sequence (vector), here called "frame", and I have plotted it. Here the code:

The resulting graph is:

I do not understand the position of the peak in this graph: I know that the peak of the autocorrelation should be in the origin, why is it there?

The x axis is simply an index number, and this is showing your the correlation for leading and trailing shifts. The shift at index 80 is lag = 0, the origin you are interested in. Type help(xcorr) in Matlab to get further detail in how to return the lag. Following that, if you do the following you will get the result you are looking for:

r, lag = xcorr(frame, frame);
figure()
plot(lag, r)

I did it as shown above to highlight how to return lag from xcorr specifically while minimizing change to your code, but xcorr(frame) will return the same result as xcorr(frame,frame) so can be further simplified for the case of autocorrelation.


Understanding a lowpass - comb filter implementation I try to understand the implementation of the low-pass comb filter of the Freeverb reverberation algorithm:

https://ccrma.stanford.edu/~jos/pasp/Lowpass_Feedback_Comb_Filter.html

The original implementation is as follows:

inline float comb::process(float input)
{
    float output;
    output = buffer[bufidx];
    filterstore = (output*damp2) + (filterstore*damp1);
    buffer[bufidx] = input + (filterstore*feedback);
    if(++bufidx>=bufsize) bufidx = 0;
    return output;
}

I drew the closed loop of this algorithm:

In the link, it says:

  Inspection of comb.h in the Freeverb source shows that Freeverb's
  ``comb'' filter is more specifically a lowpass-feedback-comb filter
  (LBCF4.11--§2.6.5). It is constructed using a delay line whose output
  is lowpass-filtered and summed with the delay-line's input. The
  particular lowpass used in Freeverb is a unity-gain one-pole lowpass
  having the transfer function

  $\displaystyle H(z) = \frac{1-d}{1-d\,z^{-1}}. $

  When $ d=0$ , the LBCF reduces to the feedback comb filter (FBCF) of
  §2.6.2 in which the feedback was not filtered. The overall LBCF
  transfer function is then

  $ \displaystyle \hbox{LBCF}_{N}^{\,f,\,d} \;= \; \frac{z^{-N}}{1 -f\frac{1-d}{1-d\,z^{-1}}\,z^{-N}}. $

Apparently, this transfer function is implemented here. But can someone tell me, how to derive this implementation out of this transfer function?

Your diagram looks correct. Let's call the transfer function in the feedback loop $G(z)$. Consider the signal $w[n]$ at the input to the delay line. Its $\mathcal{Z}$-transform satisfies

$$W(z)=X(z)+G(z)Y(z)\tag{1}$$

where $X(z)$ and $Y(z)$ are the $\mathcal{Z}$-transforms of the input and output sequences, respectively. The output is just a delayed version of $w[n]$, i.e.,

$$Y(z)=W(z)z^{-M}\tag{2}$$

From $(1)$ and $(2)$, the transfer function of the complete system can be derived as

$$H(z)=\frac{z^{-M}}{1-G(z)z^{-M}}\tag{3}$$

Now it remains to find $G(z)$. It is quite straightforward to show that

$$G(z)=\frac{d_2f}{1-d_1z^{-1}}\tag{4}$$

Combining $(3)$ and $(4)$ gives the desired result with $d_1=d$ and assuming that $d_2=1-d_1$.


What is the difference between convolution and cross-correlation? I've found on multiple sites that convolution and cross-correlation are similar (including the tag wiki for convolution), but I didn't find anywhere how they differ.

What is the difference between the two? Can you say that autocorrelation is also a kind of a convolution?

The only difference between cross-correlation and convolution is a time reversal on one of the inputs. Discrete convolution and cross-correlation are defined as follows (for real signals; I neglected the conjugates needed when the signals are complex):

$$
x[n] * h[n] = \sum_{k=0}^{\infty}h[k] x[n-k]
$$

$$
corr(x[n],h[n]) = \sum_{k=0}^{\infty}h[k] x[n+k]
$$

This implies that you can use fast convolution algorithms like overlap-save to implement cross-correlation efficiently; just time reverse one of the input signals first. Autocorrelation is identical to the above, except $h[n] = x[n]$, so you can view it as related to convolution in the same way.

Edit: Since someone else just asked a duplicate question, I've been inspired to add one more piece of information: if you implement correlation in the frequency domain using a fast convolution algorithm like overlap-save, you can avoid the hassle of time-reversing one of the signals first by instead conjugating one of the signals in the frequency domain. It can be shown that conjugation in the frequency domain is equivalent to reversal in the time domain.


Implementing an offline guitar monophonic pitch detector I am implementing an offline guitar monophonic pitch detector. I was wondering about the trade-off of the standard methods: Autocorrelation, HPS, Cepstrum, picking the peak from an STFT, the YIN algorithm.

I am working with audio that could will be segmented according to onset times, but let's assume this is done, so let's assume that I am given 5 segments which have been correctly segmented (start at the onset time, finish a frames afterwards).

My first question is: How long should a be? I guess that for algorithms like Autocorrelation, the bigger the segment the better?

Also, when working with guitar in particular, is there any issue I should take care of and is one of those methods preferable? I would assume HPS is good as the guitar produces a lot of harmonics.

I tried implementing Autocorrelation but I get a lot of octave errors. Would the YIN algorithm correct many of those? I also implemented a simple peak picking from an STFT but I often get strongest harmonics instead of the perceived pitch.

Each of the "standard methods" you listed has a failure mode (or modes) and an estimation error (which may vary with the guitar and the note). Characterize the failure modes or each method, and evaluate the estimation errors of the more successful methods, then match them to your samples, and perhaps pick the least bad (or other suitable statistical fit.)

With "the guitar in particular", one thing to watch for is that some notes on some guitars can produce a slightly (but measurably) inharmonic overtone series.  If so, the waveform will never be perfectly periodic.  This can cause successive periods of the fundamental mode waveform to never exactly match (which can throw off both autocorrelation and AMDF lag estimation methods, etc.).  Thus, one may need a psychoacoustic model to determine what might likely be the perceived pitch to a human (assuming that's what you care about).  Also, the ratios between overtone magnitudes will usually change (evolve) over time, not just decay linearly as per an ADSR envelope model, which might also affect periodicity and your choice of method.


How to get the FFT of a sine wave I am trying to sample a sine wave and plot it's frequency components, but I am having problems implementing it.The result of taking 65536 samples of one cycle of a sine wave with max amplitude 1 and a frequency 100 can be seen below.Where the Y-axis this the magnitude of the complex Fourier sum, and the x-axis is the sample number.How can I see what frequency the sine wave has from this plot?, I would like to have a plot of freq vs Magnitude and not sample Number vs magnitude.

I was expecting one big spike of amplitude 1 or 2 spikes of amplitude 0.5 each, but I seem to be getting a large amplitude, I don't know what i might be doing wrong.

Matlab Code:

Amp = 1;
freq = 100;
dt = 2 * pi /65536;
index = 1;
for t = 0:dt:2*pi
   sine(index) = Amp * sin(2*pi*freq*t);
   sampleNumber(index) = index;
   index = index + 1;
end

transform = fft(sine);
magTransform = abs(transform);
plot(sampleNumber,magTransform);

It looks like you're getting your frequency in Hertz confused with in radians/sec since you have the factor of $2\pi$ in both your sampling period dt and your signal. I rewrote a bit of your code to clarify what I think you really want.

Amp = 1;
freqHz = 12000;
fsHz = 65536;
dt = 1/fsHz;
t = 0:dt:1-dt;
sine = Amp * sin(2*pi*freqHz*t);

N = 65536;
transform = fft(sine,N)/N;
magTransform = abs(transform);

faxis = linspace(-fsHz/2,fsHz/2,N);
plot(faxis/1000,fftshift(magTransform));
axis([-40 40 0 0.6])
xlabel('Frequency (KHz)')

If your sampling frequency is 65536 samples/second, and you want for example a tone at 12 KHz, you can create it as shown. So here your sample period is 1/65536 seconds.

Your expectation of getting two spikes with an amplitude of 0.5 each was correct -- just your generated tone was not.

As for scaling the x-axis to be in Hertz, just create a vector with the same number of points as your FFT result and with a linear increment from $-fs/2$ to $+fs/2$. Note also the fftshift I used in the plot. That's because the output of Matlab's FFT function goes linearly from 0 to fs. I find it easier to visualize having DC centered but either way is fine. Without the fftshift the faxis vector would go from 0 to fs.


About the fourier transform of Periodic Signal I cannot understand some part of the period signal's Fourier transform.

Here this my note's methods,

  For periodic signal with period $T_0$, define as $s_{T_0}(t)$ as

  $$
s_{T_0} (t) =
\begin{cases}
s(t) for -T_0/2 <t<T_0/2\\
0 & \text{otherwise.}
\end{cases}
$$

  $$s(t)=\sum_{n=-\infty}^{\infty}s_{T_0}(t-n T_0) = s_{T_0}(t) * \sum_{n=-\infty}^{\infty} \delta (t-n T_0)$$
  So the fourier transform will be
  $$S(f)=S_{T_0}(f) \cdot f_0 \sum_{n=-\infty}^{\infty}\delta(f-n f_0) = f_0 \sum_{n=-\infty}^{\infty}S_{T_0}(n f_0) \delta(f-n f_0)$$

and the other methods is

  $$x(t)= s(t) \cdot \sum_{n=-\infty}^{\infty} \delta (t-n T_0)$$
  $$X(f)=S(f) * f_0 \sum_{n=-\infty}^{\infty}\delta(f- n f_0) = f_0 \sum_{n=-\infty}^{\infty}S(f-n f_0)$$

I cannot understand why I have "$\cdot f_0$" after Fourier Transform and why "$\delta(f-n f_0)$" can transform to "$\delta(t-nT_0)$"? And the last questions is, for the periods signal, why the note can use fourier transform directly with out considering the fourier serials?

It is not obvious what exactly OP Samuel wants to know, and so I am posting a generic answer.

Let $x(t)$ denote a finite-energy continuous-time signal by which we mean that
$\displaystyle \int_{-\infty}^\infty |x(t)|^2\,\mathrm dt = \mathbb E_x < \infty$.
It is important to understand that a finite-energy signal cannot be a periodic signal, that
is, there is no nonzero real number $T$ for which the following equation
can hold for all time instants $t$:
$$x(t+T) = x(t)~~\text{for all}~ t, -\infty < t < \infty.$$
Of course, the above might hold for some, possibly even
many, values of $t$,
but it cannot hold for all $t$.  For example, if
$x(t) =\begin{cases}\sin(200\pi t), &-1 \leq t \leq +1,\\
0, &\text{otherwise,}
\end{cases}$ is a sinusoidal pulse of duration $2$ seconds
containing $200$ cycles
of a $100$ Hz sine wave, then
$x(t+0.01) = x(t)$ for all $t$ except for $t \in (-1.01,1)$ and
$t \in(0.99,1)$.
On the other hand, we do not insist that $x(t)$ be of finite duration,
or that $x(t)$ be causal, etc. In these circumstances, $x(t)$ has
a Fourier transform $X(f)$ given by
$$X(f)
= \int_{-\infty}^\infty x(t)\exp(-j2\pi ft)\,\mathrm dt,
~~ -\infty < f < \infty \tag{1}$$
where $j = \sqrt{-1}$,
and $x(t)$ can be recovered from its Fourier transform using the inverse
Fourier transform which gives
$$x(t)
= \int_{-\infty}^\infty X(f)\exp(j2\pi ft)\,\mathrm df,
~~-\infty < t < \infty.\tag{2}$$

All this, of course, is old hat to most readers, so now let us
construct a periodic signal from $x(t)$.  Let $T_0$ denote
some positive number, and define

  $$s(t) = \sum_{n=-\infty}^\infty x(t - nT_0).\tag{3}$$

Then, we claim that $s(t)$ is a signal of period $T_0$.
Note that for every choice of real number $t$, we have that
$$\begin{align}
s(t+T_0) &= \sum_{n=-\infty}^\infty x(t+T_0-nT_0)\\
&= \sum_{n=-\infty}^\infty x(t-(n-1)T_0) &\text{and now replace}~ n-1~\text{by}~m\\
&= \sum_{m=-\infty}^\infty x(t-mT_0)\\
&= s(t).
\end{align}$$
But, what does one period of $s(t)$ look like? Well, for
$-T_0/2 \leq t \leq T_0/2$, $s(t)$ is not the same as $x(t)$
except in the special case when the support of $x(t)$
is a subset of $[-T_0/2,T_0/2]$. Instead, one period (what OP
Samuel calls $s_{T_0}(t)$) can be found by slicing up
the graph of $x(t)$ into pieces of duration $T_0$ centered at
integer multiples of $T_0$, moving all the pieces to be centered
at $0$, and then adding them up.  This sum, one period of the
periodic signal that we have created, is now replicated at
intervals of $T_0$ seconds all along the time axis. Thus,
we have another representation of $s(t)$, namely,

  $$s(t) = \sum_{n=-\infty}^\infty s_{T_0}(t - nT_0).
\tag{4}$$

Let us at this point introduce Dirac deltas or impulses $\delta(t)$
which have the property that for any $g(t)$ that is continuous
at $t = a$,
$$\int_{-\infty}^\infty g(t)\delta(t-a)\,\mathrm dt
= \int_{-\infty}^\infty g(t)\delta(a-t)\,\mathrm dt = g(a).$$
Now, the convolution of $g(t)$ and $h(t)$ is a signal $(g\star h)$
whose value $(g\star h)(\tau)$ at particular time instant $\tau$ is given by
$$(g\star h)(\tau) = \int_{-\infty}^\infty g(t)h(\tau-t)\,\mathrm dt$$
and so if $h_n(t) = \delta(t-nT_0)$, we have that
$$\begin{align}
(x\star h_n)(\tau) &= \int_{-\infty}^\infty x(t)h_n(\tau-t)\,\mathrm dt\\
&= \int_{-\infty}^\infty x(t)\delta(\tau-t-nT_0)\,\mathrm dt\\
&= \int_{-\infty}^\infty x(t)\delta((\tau-nT_0) - t)\,\mathrm dt\\
&= x(\tau - nT_0).
\end{align}$$
This gives us yet another representation of the periodic signal $s(t)$,
namely,

  $$s(t)= x(t) \star \sum_{n=-\infty}^\infty \delta(t-nT_0)
= \sum_{n=-\infty}^\infty x(t-nT_0). \tag{5}$$

It is worth reminding ourselves at this point that $s(t)$ is
a continuous-time finite amplitude signal that, despite
what it looks like in the middle sum in $(5)$ does not
contain any impulses.

Returning to Fourier transforms (bet you thought I forgot about
them, didn't you?), we come up against the uncomfortable fact
that $s(t)$ does not have a Fourier transform in the usual
sense of the word because it is not a finite-energy signal.
What it does have is a Fourier series representation,
as pointed out by Fourier himself almost two hundred years ago.
We can incorporate Fourier series into Fourier transforms
if we use impulses in the frequency domain. So, what is
the Fourier series for $s(t)$? We write
$$s(t) = \sum_{k=-\infty}^\infty c_k \exp(j2\pi kt/T_0), -\infty < t < \infty$$
where
$$\begin{align}
c_k &= \frac{1}{T_0}\int_{-T_0/2}^{T_0/2}s(t)\exp(-j2\pi kt/T_0)\,\mathrm dt\\
&= \frac{1}{T_0}\int_{-T_0/2}^{T_0/2}\sum_{n=-\infty}^\infty x(t-nT_0)
\exp(-j2\pi kt/T_0)\,\mathrm dt  &\text{using}~ (3)\\
&= \frac{1}{T_0}\sum_{n=-\infty}^\infty \int_{-T_0/2}^{T_0/2}x(t-nT_0)
\exp(-j2\pi kt/T_0)\,\mathrm dt &\text{now change variables}~ \tau = t-nT_0\\
&= \frac{1}{T_0}\sum_{n=-\infty}^\infty \int_{nT_0-T_0/2}^{nT_0+T_0/2}
x(\tau) \exp(-j2\pi k(\tau+nT_0)/T_0)\,\mathrm d\tau \\
&= \frac{1}{T_0}\sum_{n=-\infty}^\infty \int_{nT_0-T_0/2}^{nT_0+T_0/2}
x(\tau) \exp(-j2\pi k\tau/T_0)\,\mathrm d\tau
&\text{now combine the integrals into one}\\
&=\frac{1}{T_0}\int_{-\infty}^\infty x(\tau)\exp(-j2\pi k\tau/T_0)\,\mathrm d\tau\\
&= \frac{1}{T_0}X\left(\frac{k}{T_0}\right) &\text{using the definition}~ (2)
\end{align}$$

So we have the Fourier series representation
$$s(t)
= \frac{1}{T_0}\sum_{k=-\infty}^{\infty}
X\left(\frac{k}{T_0}\right)\exp(j2\pi kt/T_0).\tag{6}$$
What should $S(f)$ be so that when it is plugged into the
inverse Fourier transform $(2)$ will give us
$$s(t)
= \int_{-\infty}^\infty S(f)\exp(j2\pi ft)\,\mathrm df
= \frac{1}{T_0}\sum_{k=-\infty}^{\infty}
X\left(\frac{k}{T_0}\right)\exp(j2\pi kt/T_0)??$$
The answer is readily worked out to be

  $$\begin{align}
S(f) &= \frac{1}{T_0}\sum_{k=-\infty}^{\infty}
X(f) \delta\left(f-\frac{k}{T_0}\right)
= \frac{1}{T_0}\sum_{k=-\infty}^{\infty}
X\left(\frac{k}{T_0}\right) \delta\left(f-\frac{k}{T_0}\right)\\
&= f_0 \sum_{k=-\infty}^{\infty} X(kf_0)\delta(f-kf_0)
\end{align}$$

I will not say much about OP Samuel's "other method" since
it is comparing apples and oranges. The signal
$$x(t)\cdot \sum_{n=-\infty}^\infty \delta(t-nT_0)
= \sum_{n=-\infty}^\infty x(nT_0)\delta(t-nT_0)\tag{7}$$ is a
impulse train in the time domain (one way of
representing a sampled-data or discrete-time signal with sampling interval
$T_0$), and is very different
from the signal
$$x(t) \star \sum_{n=-\infty}^\infty \delta(t-nT_0)
= \sum_{n=-\infty}^\infty x(t-nT_0)$$ in $(5)$ which is a continuous-time
signal containing no impulses.
There is no reason for the two to have similar or identical Fourier transforms,
and the OP's Procrustean efforts to equate both is leading him
into a great deal
of confusion.
In fact, duality suggests that, just as the Fourier transform of
a periodic signal is a set of equally-spaced impulses (of different
amplitudes) in the frequency domain, the Fourier transform of
a set of equally-spaced impulses (of different amplitudes) in
the time domain is a periodic function in the frequency domain.
For details, see  this answer
to a different question in which the Fourier transform of $(7)$ has been
described in more depth than the hand-waving in the previous sentence.


Removing 'high-frequency' components from an array? I am new to the stream of digital signal processing, and hence my doubt may seem extremely stupid, but:
In tools such as matlab, signals are stored in a single array. If i convolute it with an LPF array, i get another array/signal, that has certain elements attenuated, can i find out which elements these are and completely remove them from the new array.

Does my question even make sense? :/

You mention convolution which implies that your arrays represent time domain signals. As this is digital signal processing, the time is discretized and every element of the array corresponds to a sample. Two consecutive samples have a fixed time distance, say $T$.

If you convolve the discrete time signal with the impulse response of a low-pass filter (LPF) the resulting ouput signal will be longer than the input signal, because in general, linear time-invariant (LTI) systems have a memory (length of impulse response). You cannot reduce the array length without distorting the output signal (unless you use block-wise processing with overlap-add/-save, but that's another story). So, no matter what type of filter the impulse response of your LTI system represents, you cannot just remove elements from the output array.

To find out which frequency components are attenuated by a system you should analyze its transfer function (sometimes also refered to as frequency response). The transfer function is obtained by transforming the impulse reponse into the frequency domain. A tool for doing that with DSP is the discrete Fourier transform (DFT) which almost always comes in form of its FFT implementation (Fast Fourier transform). In case that the system in question is a black box and its impulse response is unknown the frequency domain representation of input and output signal can be compared to find out what frequencies are attenuated or amplified.


How does an Image filter really sweep an Image? Suppose, we have the following image:

Now, we want to apply a 5x5 filter to the image.

How would the filter sweep the image?

Like this?

Or, like this?

Or, anything else?

I assume by image filter sweeping, you mean the image is convolved with a kernel.  In that case, the method is more similar to what you show in your second image.

Conceptually, the NxM kernel is multiplied by an NxM block of pixels in the image, and the results are summed to a single value-- this is the first value in your 2D result array.

Then the next NxM block of pixels is taken by offsetting the window one pixel to the right, the kernel is multiplied by this block, the results summed, and this is the second value (top row, 2nd from left) in your result array.

This continues until you reach the end of the row, then you offset down one pixel and start again at the beginning of the row.  This result is then the first result in the second row of your result array.

And so forth.

Of course, it doesn't have to be done in order, and many optimizations or parallel programming techniques may perform better if you don't step through sequentially overlapped like this, but conceptually this is what's happening-- take a block, convolve, get one pixel output, put it into the result for that offset.

There is also the question of how to handle edges-- if you want the convolution result to be the same size as the input image, you need to assume values for pixels outside the border of the actual image.  Zero is a common assumption.

Considering edges, and in reference to the question in the comments, where you start the kernel depends on what exactly you want out.  First, different kernels may have different "origins"-- typically the result of the kernel calculation is referenced to the center of the kernel, but that's not a requirement (and can't happen if the kernel is an even number of units across).

Then the tradeoff needs to be resolved between the size of the output and the validity of the data.  If you want an output image of the same size as the input image, then starting (3,3) of the kernel at (1,1) in the image, assuming 1 based indexing, will give you that result.  The challenge is what value to use for parts of the kernel that don't overlap the image.  For example, if you are using a kernel to average a neighborhood of pixels, and you assume 0 for values outside the image, then the edge-most pixels of your result image will darken because you are averaging in those zeros with the existing data.  That may be reasonably harmless.

Other kernels, however, may give bad data at the border if you assume values outside the image.  Assuming zero outside the image when you're taking a picture of sky could cause "phantom" edges to appear with an edge detection kernel, for example, and this could cause problems later in your processing pipeline.  In this case you might prefer to start with (3,3) of your kernel at (3,3) of the image, and then your result image will be 4 pixels narrower and 4 pixels shorter.


Is sound affected by wind force? Suppose wind direction is from point A to point B and some sound source like glass shattering is at point B, will sound waves going to point A be affected by that wind coming in opposite direction, or will the sound waves cut right through?

will sound waves going to point A be affected by that wind

Yes. There are two major effects of wind on sound.

Noise: Wind is moving air and moving air is sound. It's mostly DC air movement so it's not audible per se, but it when the airflow hits objects it often creates wide band turbulences which can be quite loud. Particularly problematic is wind hitting the receivers i.e. a microphone or your ears. It can generate a lot of pick up noise. Wind screens are one way to deal with this.
Waveform shaping: Sound waves are moving with resect the fluid media (i.e. air). If the media moves as well, than the overall movement of the sound wave will be the sum of both motions and this deforms the wave shape. If the wind is relative constant you still get roughly spherical waves, but they are not concentric anymore, i.e. the center is moving with the wind. So your glass shatter would be significantly louder downwind then upwind.


What is the autocorrelation of a Dirac pulse? What is the autocorrelation of $x(t) = \delta(t)$?

Can you explain to me how to calculate it?

Well, by definition of the $\delta$ distribution, you have:

$\int_{-\infty}^{\infty} f(t) \delta(t-T)\, \textrm{d}t = f(T)$

The autocorrelation of a function $g(t)$ can be computed via:

$\int_{-\infty}^{\infty} g^{*}(t)g(t + \tau)\, \textrm{d}t$, with $g^*$ as the complex conjugate of $g$. Since $\delta(t)$ is real-valued, this is conjugation can be skipped. So you are left with:

$\int_{-\infty}^{\infty} \delta^{*}(t)\delta(t + \tau)\, \textrm{d}t = \int_{-\infty}^{\infty} \delta(t)\delta(t + \tau)\, \textrm{d}t = \delta(-\tau)$.

The first = sign comes from the autocorrelation of the real-valued $\delta$, the second from the definition of the $\delta$-distribution.

So, the autocorrelation function of the $\delta$-distribution is the distribution itself. A eigenfunction of the autocorrelation function, so to say ;)

Think about it, this does make sense: the only perfect match is achieved with no time shift, ie at $\tau = 0$. All other shifts would end up with one of the arguments of the $\delta$ being different from 0, hence with the $\delta$-function being 0 there.

BTW: $\delta(-\tau) = \delta(\tau)$, since the function/distribution is symmetric.


Help in sketching a CT signal with discontinuities This is the Original signal:

And I'm supposed to sketch
$[x(t) + x(-t)] .u(t)$
I know that unit step signal will make the signal starts from 0 to $+\infty$
My problem is at key point (0) and key point (1)

You are right, the solution manual is wrong and here is a piecewise linear definition of the signal $[x(t)+x(-t)]u(t)$ based on the acceptance of $\text{u}(0) = 0.5$ (i.e. the value at a discontinuity is the half way between right and left ends)
$$[x(t)+x(-t)]u(t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  t < 0}}\\&[1.5] &\scriptstyle{\text{, t = 0}}\\&3 &\scriptstyle{\text{, 0<t<1}}\\&[2] &\scriptstyle{\text{, t = 1}}\\-&2t+3 &\scriptstyle{\text{,  1<t<2}}\\&[-0.5] &\scriptstyle{\text{, t=2}}\\&0 &\scriptstyle{\text{,  t > 2}}\\
\end{align}
\end{cases}$$

Where the values within $[.]$ are points of discontinuities.

Let me put the steps to reach this conclusion in more details:
Looking at the provided graph of $x(t)$ we can deduce a piecewise definition of the signal. However we must be causious about the possible discontinuities, for the sake of mathematical (continuous) function theory, those isolated points of discontinuities can be omitted entirely or as an engineering convention the value of the signal can be defined to be the half of the sum of its values at the right and left limits of the function at the discontinuity point. I'll accept this second definition, as DSP is concerned more with engineering than with pure mathematics.

Therefore, we have the definition of $x(t)$ as follows:
$$x(t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  t < -2}}\\
&[-0.5] &\scriptstyle{\text{, t = -2}}\\
&t+1 &\scriptstyle{\text{, -2<t<-1}}\\
&[0.5] &\scriptstyle{\text{, t = -1}}\\
&1 &\scriptstyle{\text{,  -1<t<0}}\\
&[1.5] &\scriptstyle{\text{, t=0}}\\
&2 &\scriptstyle{\text{,  0<t<1}}\\
&[1.5] &\scriptstyle{\text{,  t=1}}\\
&-t+2 &\scriptstyle{\text{,  1<t<2}}\\
&0 &\scriptstyle{\text{,  t>2}}\\
\end{align}
\end{cases}$$

From which we can also deduce $x(-t)$ as follows:
$$x(-t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  -t < -2}}\\
&[-0.5] &\scriptstyle{\text{, -t = -2}}\\
&-t+1 &\scriptstyle{\text{, -2<-t<-1}}\\
&[0.5] &\scriptstyle{\text{, -t = -1}}\\
&1 &\scriptstyle{\text{,  -1<-t<0}}\\
&[1.5] &\scriptstyle{\text{, -t=0}}\\
&2 &\scriptstyle{\text{,  0<-t<1}}\\
&[1.5] &\scriptstyle{\text{,  -t=1}}\\
&t+2 &\scriptstyle{\text{,  1<-t<2}}\\
&0 &\scriptstyle{\text{,  -t>2}}\\
\end{align}
\end{cases}$$

Observe carefully how every occurence of $t$ is replaced with $-t$ both in the defining formula of the function and also in the definition of its intervals. We shall reformat intervals of $x(-t)$ with positive $t$ as follows:

$$x(-t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  t > 2}}\\
&[-0.5] &\scriptstyle{\text{, t = 2}}\\
&-t+1 &\scriptstyle{\text{, 1<t<2}}\\
&[0.5] &\scriptstyle{\text{, t = 1}}\\
&1 &\scriptstyle{\text{,  0<t<1}}\\
&[1.5] &\scriptstyle{\text{, t=0}}\\
&2 &\scriptstyle{\text{,  -1<t<0}}\\
&[1.5] &\scriptstyle{\text{,  t=-1}}\\
&t+2 &\scriptstyle{\text{,  -2<t<-1}}\\
&0 &\scriptstyle{\text{,  t<-2}}\\
\end{align}
\end{cases}$$

Now adjust the order of intervals of $x(-t)$ so that they are in left to right order:

$$x(-t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  t<-2}}\\
&t+2 &\scriptstyle{\text{,  -2<t<-1}}\\
&[1.5] &\scriptstyle{\text{,  t=-1}}\\
&2 &\scriptstyle{\text{,  -1<t<0}}\\
&[1.5] &\scriptstyle{\text{, t=0}}\\
&1 &\scriptstyle{\text{,  0<t<1}}\\
&[0.5] &\scriptstyle{\text{, t = 1}}\\
&-t+1 &\scriptstyle{\text{, 1<t<2}}\\
&[-0.5] &\scriptstyle{\text{, t = 2}}\\
&0 &\scriptstyle{\text{,  t > 2}}\\
\end{align}
\end{cases}$$

Now we can add these two signals for each respective interval. In fact we are here a bit lucky that the signals do have matching intervals, which would not be the case if the original signal $x(t)$ involved nonsymetric pieces in the first place. The addition produces the following signal:

$$x(t)+x(-t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{,  t<-2}}\\
&[-0.5] &\scriptstyle{\text{,  t=-2}}\\
&(t+1)+(t+2)=2t+3 &\scriptstyle{\text{,  -2<t<-1}}\\
&[0.5]+[1.5]=[2] &\scriptstyle{\text{,  t=-1}}\\
&1+2=3 &\scriptstyle{\text{,  -1<t<0}}\\
&[1.5]+[1.5]=[3] &\scriptstyle{\text{, t=0}}\\
&2+1=3 &\scriptstyle{\text{,  0<t<1}}\\
&[1.5]+[0.5]=[2] &\scriptstyle{\text{, t = 1}}\\
&(-t+1)+(-t+1)=-2t+3 &\scriptstyle{\text{, 1<t<2}}\\
&[-0.5] &\scriptstyle{\text{, t = 2}}\\
&0 &\scriptstyle{\text{,  t > 2}}\\
\end{align}
\end{cases}$$

Finally we shall multiply this with the funciton $\text{u}(t)$ whose piecewise definition is:
$$\text{u}(t) =
\begin{cases}
\begin{align}
&0 &\scriptstyle{\text{, t<0}}\\
&[0.5] &\scriptstyle{\text{, t=0}}\\
&1 &\scriptstyle{\text{, t>0}}\\
\end{align}
\end{cases}$$

And the resulting function $[x(t)+x(-t)]\text{u}(t)$ was already provided in the first part above...


Continuously changing pitch of stable-pitch audio sample I have an 44100-rate 4 second sample of a legato violin playing at 440 Hz.  I would like to modulate this sample so that it starts at 440 Hz, then transforms its frequency continuously and exponentially to 660 Hz, while maintaining the original timbre.
I have tried pitch shifting directly from 440 to 660 and this works fine, but trying to simply cut the sample into extremely short bits and increase the pitch shift for each consecutive bit resulted in a lot of artifacts.
Do you have any suggestions?
Bonus points if it would be easy to implement in python!

You're trying to do the time warp dance.  If you were copying the samples point by point, your playback time (destination) and your look up time (source) would be the same.  If you wanted to speed it up by a factor of two(one octave), your lookup time (source location in the array) would be moving twice as fast as your sample time (location in your output buffer).  For variable speed time, you have to figure out where in the lookup table to grab the next value (function range).  Where you are going to write it to is the next spot in the buffer, as in uniform time (function domain).  When you calculate the location to get your signal value, it may not land on an integer.  You could get some rough sound if you just use the nearest one.
So, you got two choices for better quality.

Increase the density of your lookup table

Use some method (called interpolation) to figure out what the value should be

A. Linear interpolation means draw a line between the two nearest samples and use that height whereever your fraction takes you
B. Polynomial interpolation means use the nearby samples to draw a nice smooth curve, and use that height
C. Sinc interpolation means use the sinc function as a weighted average of many nearby points to get the value
Your musician friends will tell you it doesn't sound right.  Instruments are complicated and their timbre varies by frequency and the middle name of the nearest bartender.
That's really the precursor to DSP.  Trying to fix the sound so it sounds right is a challenge you can have hours of fun pursuing, as others are.  There is lots of material if you look for it.  I am not really an expert on that at all.
In a synthesizer wavetable the source location just goes round and round and round...  how much faster each step is related to the frequency.
Keep a steady tone by not changing the step size.
Linear growth in time.  Add the same amount to each step.
Exponential growth in time.  Add a multiple of the previous step size.
For instance,
ThisStep = 1.00001 * PreviousStep
SourceLocation += ThisStep
The step size grows exponentially and so does its accumulation.
Speeding up by a factor of $2^{1/12}$ is one semitone.
There's math problem and solution in there somewhere, I'm sure you'll spot it and solve it to get this done.  Or maybe two.


Is sampling a Fourier transformed signal and fourier transforming a sampled signal the same? I'm having a hard time understanding an assignment that states:

  Draw the complex spectrum of the sampled signal $f(t)$ (periodic and continuous). Do this, by first calculating the Fourier transformation and sample it afterwards by multiplying with the impulse train.

The way I understand it, I need to calculate $$\begin{equation*} \mathcal{F}\left[f(t) \cdot \sum_{k=-\infty}^{\infty}\delta(t-kT)\right]
\end{equation*}$$
however the second sentence suggests this is the same as
$$\mathcal{F}[f(t)] \cdot \sum_{k=-\infty}^{\infty}\delta(t-kT)
$$
However, Wikipedia says this is not the case and instead suggests
$$\begin{equation*}\mathcal{F}\left[f(t) \cdot \sum_{k=-\infty}^{\infty}\delta(t-kT)\right] = \mathcal{F}[f(t)] *\left[\frac{1}{T} \cdot \sum_{k=-\infty}^{\infty}\delta\left(t-\frac{k}{T}\right)\right]
\end{equation*}$$
are maybe both formulas correct? and how do I draw a complex spectrum? is it 3D?

If I number your equations (1), (2) and (3) starting from the first, then your Eq (1) is correct as it represents a conceptual model of the ideal sampling process: multiplying the continuous signal with an impulse train. Consequently, Eq (3) is also correct since multiplication in time domain is convolution in frequency domain (recall that the Fourier transform of an impulse train in time domain is also an impulse train in frequency domain with a different scaling factor). Convolving the two is what causes the spectrum to repeat at multiples of the sample rate.

Your Eq (2) however is not entirely correct. Fourier transforming the signal first and sampling it later implies sampling in frequency domain. Interestingly, this should create 'aliases' in time domain, i.e., repetition of the same waveform at multiples of 'sample rate'.

Why it is not completely incorrect here is that your signal is given to be continuous and periodic. A periodic signal has a sampled spectrum due to the reason described above; technically, it is said to have a Fourier series. So when you compute $F(f)$, either you should have in time domain those complex sinusoids with Fourier series coefficients, or you should have in frequency domain a frequency domain impulses with those coefficients. This might have caused confusion to whoever made the assignment.

Update: A complex spectrum can be drawn in 3D with frequency, I and Q dimensions, but the most common practice (and that reveals the most information) is to draw two separate graphs: one for magnitude response against frequency and the other phase response against the frequency.


Non-separable, higher-order generalized cosine windows in 2D for images Higher-order generalized cosine windows define a useful parametric family of windows functions, embedding Hamming, Hann, Blackman, Nutall, etc.:
$$w(n) = \sum_{k = 0}^{K} a_k\; \cos\left( \frac{2 \pi k n}{N} \right)$$
Are there non-separable 2D similar forms, i.e. not a tensor product of the 1D form, usable on the standard rectangular lattice (excluding designs such as in Multidimensional windows over arbitrary lattices and their application to FIR filter design)?

The form:

$$w(x, y) = \sum_{j = -K}^{K}\sum_{k = -K}^{K} a_{|j|,|k|}\; \cos(\pi jx + \pi ky),$$

where:

$$a_{j,k} = a_{k,j},\\
\sum_{j = -K}^{K}\sum_{k = -K}^{K} a_{|j|,|k|} = 1,\\
x \in -1..1,\\
y \in -1..1$$

yields window functions that are 2D similar. Window functions of that form are not generally separable, but some of them are, for example:

$$w(x, y) = \frac{\cos(\pi x - \pi y)}{8} + \frac{\cos(\pi x + \pi y)}{8} + \frac{\cos(\pi x)}{4} + \frac{\cos(\pi y)}{4} + \frac{1}{4}\\
= \left(\frac{\cos\left(\pi x\right)}{2} + \frac{1}{2}\right)\left(\frac{\cos\left(\pi y\right)}{2} + \frac{1}{2}\right)$$

For this window function, if the balance between the zero-frequency term and the higher terms is altered (as in going from Hann to Hamming), separability breaks.


Spectrum in 1/3 octave bands from FFT I’m looking to use an FFT to generate a frequency spectrum in 1/3 octave bands. After reading many posts on this site as well as others, I believe the appropriate approach is set out below. I’ve dry-run the calculations in Excel and the results look reasonable, but reasonable is not necessarily correct.
To ask a specific question, is the method I’m using valid? In particular, I question if the method applied in Step 2 is correct.
Background.
I’m sampling the electric signal on the output of a pre-amp and my goal is a frequency spectrum in dB. My reference (denominator in the dB calc) is arbitrary and chosen to set an appropriate scale for the spectrum output display. I’m not attempting to take measurements from the output and don’t need to calibrate it against an objective reference. My project will be implemented in C++.
Step 1.
Start with a 1024 sample FFT at 48,000 Hz (513 bins at 46.9 Hz each). Result is the DFT, Fn for n = 0 – 512. Units are in V.
Note, normally, I’d normalize the FFT output by multiplying each bin by 2/N (N is number of samples, i.e. 1024). However, because we eventually plot the spectra in dB relative to an arbitrary amplitude, any normalizing I do here will be washed away when I compute the results in dB so I don’t bother with it in my project.
Step 2.
Convert the DFT into 1/3 octave bands.
I obtained the following formulas from here:
https://www.ap.com/technical-library/deriving-fractional-octave-spectra-from-the-fft-with-apx/
The amplitude in V for each 1/3 octave band is:

Where Lb is the amplitude for each 1/3 octave band in V, for b = 1 to 32, and gn,b is the gain multiplier for FFT bin n and 1/3 octave band b:

Where fn is the frequency of bin n and fb is the center frequency of band b. k is the octave bandwidth designator, 1 for full octave and 3 for 1/3 octave.
Result is Lb for b = 1 to 32 in V.
Step 3.
Disregard lower bands.
Below about 250 Hz, for an FFT of length 1024, there are not enough FFT bins within each band to make the results meaningful. Therefore, we disregard the lower bands and keep only the 1/3 octave bands starting with b13 = 251 Hz.
Step 4.
Convert the amplitude data in each band to dB.
For this computation, I’m computing the dB value relative to the approximate maximum amplitude likely to appear in any frequency band. This will be Lref.
Final value in dB for each 1/3 octave frequency band above 250 Hz is:

Below is a sample of individual gn,b terms (FFT bins down the left, 1/3 octave bands across the top):

Assistance, input, corrections and comments appreciated.
Thanks!

It seems like you’ve by and large got it, but I did notice a couple things. You’ll probably want to include support for windowing and the corresponding correction factor as shown by APs article on FFT noise.  Rectangular windows aren’t great for RMS type measurements (at least in my opinion) for all the side lobe type explanations that are common with windows.  Second, I’m pretty sure you need to sum from 1 to N, not N/2.  You may have realized that because the input signal is real, that there is redundant data that doesn’t need to be calculated twice. That is true, with the exception of the 0 and N/2 (if N is even) bins.  Granted this probably wouldn’t make the output too wrong, but it would be off a bit.
Edit:  you know what, AP may have designed those coefficients so that windowing either isn’t necessary, or would do more harm than good.  You may still want to try it, since it’s easy enough to support, but it may not help.  Also, I don’t see how this would be more computationally efficient than doing the calculations directly in the time domain.


How can I calculate the Number of bins per octave? I am trying to calculate the Constant Q Transform for Monophonic audio using Python.
The Python has the library to calculate the CQT , the default value of num_bins_per_octave is 12.

How would one decide the number bins?
Newbie to the Topic. Please enlighten me.

I do not know how to determine the ideal number of bins. However, I can shed some light on why in Librosa's CQT function, the default values are num_bins_per_octave=12 and n_bins=84. (This is my first answer so any feedback would be helpful. And because this does not directly answer the question, should I perhaps post this as a comment?)
The Librosa CQT function you mention can be found in the Librosa documentation here.
First some music theory
Western music consists of 12 identifiable types of pitches. Each of these 12 pitches is considered a new note. The 12 pitches are as follows:
C
C# / Db
D
D# / Eb
E
F
F# / Gb
G
G# / Ab
A
A# / Bb
B
The notes with slashes can be called by either name, depending upon which key we’re in. (It doesn't matter what a key is right now.)
These 12 pitches repeat in the same order throughout the complete span of human hearing. So after B, there would be another C, then another C#, another D, and so on and so forth, giving us multiple C pitches, C# pitches, D pitches, etc.
A specific group of these 12 pitches is defined by something called an octave. An octave is the distance between one note (like C#) and the next note bearing its same name (the next C# that’s either higher or lower).
The reason why an octave is divided into 12 pitches rather than, let's say 11 or 13 pitches, can be found here.
Librosa's default num_bins_per_octave is 12 specifically because there are 12 pitches in an octave.
The piano
A great visual of multiple consecutive octaves can be seen in this 76-key piano here:  Notice that the set of 12 unique pitches repeats over and over again on this keyboard as you move left to right (or vice versa), counting both the black and white keys.
A normal piano has 88-keys, 7 octaves plus four extra notes. Looking at the 7 octaves, each octave has 12 pitches giving us a total of 7*12=84 pitches on this part of the piano. This is where the 84 from in the n_bins=84 parameter of the Librosa CQT method comes from.
There are multiple reasons why the standard piano has 88 keys instead of 84 that you can read about here.
Hope this helps.


How large FFTs can pull signals out of the noise floor? I am trying to detect unknown RF tones around -140 dBm and my scan BW is 5 MHz, going through the Noise power calculations the signal is below the thermal noise based on the scan BW. I read that using large FFTs can help to pull signals out of the noise floor. My question is how large FFTs can accomplish this?

The receiver used is an RFIC (LNA, Filter and ADC), and I was thinking FFTs of the order of 2¹⁶ or higher running on a hardware accelerator.

The unknown RF signals are ERP emissions from electronic devices which have RF circuitry.

If you look at the formula of a single DFT bin

$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j2\pi k\frac nN}\text,$$

you'll notice that his is essentially a correlation of $x$ with the complex sinusoid $e^{-j2\pi k\frac nN}$.

That means the DFT can just be understood as a filter bank of matched filters for single tones that fall in the DFT "raster".

Hence, you simply get FFT length-based processing gain: The length of the sum.

But: you probably don't have perfect knowledge of the exact frequency of the signal you're trying to detect¹! So, you can't put things into that perfect DFT raster.

Now, the larger you choose the FFT length $N$, the finer that raster will get, but also, the longer your observation has to be, and the more compute power you'll need.

At some point, the DFT stops being the best possible tone detector, and superresolution techniques become relevant. In this case (weak tone, you're sure that you've only got exactly one tone in your signal), the ESPRIT algorithm with a long observation period leading to the autocovariance matrix estimate that it takes as input, would probably work very nicely.

¹ There's inevitably frequency error in your receiver, and in your transmitter. Papers that start with We assume perfect synchronization typically skip the hard part of making a system work...


Does orthogonal beamforming mean $\mathbf f_a \cdot \mathbf f_b =0$? If i said that these two beamformings $\mathbf f_a $ and $ \mathbf  f_b $ are orthogonal to each other,does it mean $\mathbf  f_a \cdot \mathbf f_b \approx 0$ ?

Because i have see two paper about orthogonal beamforming,but they both doesn't mention this relation.i just want to make sure about that

Orthogonality might be considered in literature as background knowledge. The definition is somehow a converse version. If $\cdot $ denotes an inner product (or a dot or scalar product), and whenever
 $\mathbf  f_a \cdot \mathbf f_b = 0$, then $f_a$ and  $f_b$ are said to be orthogonal.

Note that, if $\mathbf  f_a$ is zero, it is orthogonal to everything.


Best way to find an object in a picture I'm trying to find an object in a picture, my solution is to take the picture and a photo of the object and find the maximum of the mutual covariance, this is my project on GitHub: ImageFinder

It works but if the object is rotated or transformed it doesn't work well. Which other method should I use?

I thought that a good way can be to look for the object through the colors and then try to rotate the image or to resize it while does not match, but this idea need an improvement.

There is no "best method", it will vary according to application. There are many trains of thought regarding this subject.

I've done something in that area using Cepstral Analysis. Take a look at this:

https://8b14d84720d60d831e70ec440417feb774595227-www.googledrive.com/host/0B4-fFSwiy1mGdS00SHB2Nlo3bkU/imreg_translation/imreg_translation.html

Note: I am the author of the above.

Basically, you do the following:

Sum the target image and the object
Calculate the Cepstrum of said sum.
Find the peaks (there will be 2, you have to resolve the ambiguity using something like SAD)
The location of the peak will tell you the translation you must apply to the image.

I based it off this article: http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=6128658&tag=1

I'm actually in the middle of creating a barcode (1D and 2D) open-source program using this kind of processing. The link to my bitbucket acount is here: https://bitbucket.org/DSevero/

Although the repo is private for now, I should be releasing a first version around January

:)


Help in proper notations and mathematical formulation For a given series, $S = \{-1,0,-2,1,etc\}$. The number of elements in $S$ is $N = 100$. Each $s_i$ belongs to a alphabet from a finite alphabet set $\mathcal{A} =  (a_1=0, a_2=1, a_3=2, a_4=3, a_5=-1, a_6=-2, a_7=-3, a_8=-4)$. The cardinality of $\mathcal{A}$ is $m = 8$. Each $a_i$ is associated with a probability $p_i$ and $\sum_{i=1}^m = 1$.

Part 1: I need help in the correct notations used to mathematically express the above example. This is what I have tried:

Let, $S =\{s_1,....,s_N\}$ be a random sequence / series represented as a succession of $N$ symbols from a finite alphabet set $\mathcal{A} = \{a_1, \ldots, a_m\}$, $|\mathcal{A}|= m$ where $a_i$ are real or complex numbers and  $s_n \in \mathcal{A}$, $n=1,2,\ldots,N$.

Is the above formulation and notation correct? What is the proper way to express this for a general case and not just for $m=8$

UPDATE : Part 1 has been answered.

Part 2: For the given series which is representing the data (information), I need to estimate the probability of occurrence of each symbol using MLE formulation. I think it should be like this but unsure. LEt $c_i$ be the letter counts since the data sequence $S$ appears to follow multinomial distribution. Multinomial distribution is based on counts.

The likelihood is
\begin{eqnarray*}
l(p_1,\ldots,p_{m};s_1,\ldots,s_{N}) &:=& {P}(c_1=a_1,\ldots, c_{m}=a_{m}; p_1,\ldots,p_{m})
\end{eqnarray*}

How to find $\hat{p}_i$ using MLE?

You observe a sequence $(s_1,\ldots, s_N)$ where random variables $S_k$ are i.i.d. and are drawn from a set (alphabet) $\mathcal{A} = \{a_1,\ldots, a_m\}$ with letter probabilities $P(a_i)=p_i$ for $1\leq i \leq m$ and $\sum_{i=1}^m p_i = 1$.

The likelihood function can be written as:

\begin{eqnarray*}
l( p_1,\ldots, p_m; s_1, \ldots, s_N ) &=& P( S_1 = s_1, \ldots, S_N = s_N ) \\
&=& \prod_{k=1}^N P(S_k = s_k) \\
&=& \prod_{k=1}^N \sum_{i=1}^m p_i \mathbf{1}\{s_k = a_i\}
\end{eqnarray*}
where $\mathbf{1}$ is the indicator function which evaluates to one when the condition is true, and zero otherwise.

In part 2 of your question, we need to change the random variable to letter counts instead of the sequence itself. Let $X_i$ be the number of times we observe letter $a_i$ for $1\leq i \leq m$ in our sequence of length $N$. Then $\{X_i\}_{i=1}^m$ has a multinomial distribution:

$$
P(X_1 = x_1,\ldots,X_m=x_m) = \frac{N!}{\prod_{i=1}^m x_i!} \prod_{i=1}^m p_i^{x_i}.
$$

(Convince yourself this is in fact the right probability mass function for letter counts; the argument involves some combinatorics.)

The log-likelihood function can be written as:

\begin{eqnarray*}
\log l(p_1,\ldots,p_m; x_1,\ldots, x_m) &=& \log P(X_1 = x_1, \ldots X_m = x_m) \\
&=& \log(N!) - \sum \log(x_i!) + \sum x_i \log p_i.
\end{eqnarray*}

The MLE of $p_i$ can be found by maximizing the log-likelihood. However, note that we have a constrained maximization problem here because we require $\sum p_i = 1$.

To convert this to an unconstrained maximization problem, let's introduce a Lagrange multiplier $\eta$. So the MLE of $p_i$ can be expressed as:

$$
\hat{p}_i = \arg \max _{p_i} \log(N!) - \sum \log(x_i!) + \sum x_i \log p_i + \eta \left( 1-\sum p_i \right).
$$

Setting the gradient of the objective function to zero gives:

$$
\frac{x_i}{p_i} - \eta = 0 \implies p_i = \frac{x_i}{\eta}.
$$

However, note that in order to meet the constraint, we must set $\eta = N$ (because $\sum x_i = N$).

This finally gives us the MLE:

$$
\hat{p}_i = \frac{x_i}{N}
$$

which can also be expressed as
$$
\hat{p}_i = \frac{ \sum_{k=1}^N \mathbf{1}\{s_k = a_i\} } {N}.
$$

This is intuitively pleasing! To get the ML estimate of $p_i$ count the number of times you see the letter $a_i$ in your original sequence and divide that by $N$.


Why every musical instrument sounds different than any other while listening? This is soft question.I am not expert in music. We know that English language has 26 letters (which are again grouped into vowels and consonants) , while language of music is same all over the world and it has only 7 letters (all are vowels). There are number of musical instruments. one can blindly recognize which instrument is being played just by listening to it.

Each musical instrument sounds (i am not finding proper word here in case of "sounds" ) different than any other musical instrument.Every instrument seems to be unique while listening.My basic question is why is it so?

Also,is it possible to replace one musical instrument by combining two or more other instruments playing together or like that ? (For example is it possible to obtain  same Guitar sound (i.e.sound/tones our ear listen while Guitar is being played) by  combination of other musical instruments say violin and harmonium etc played together ?)

Actually there are classes of musical instruments whose members will have similar sounds such as the string-class of violin, viola, cello or the class of woodwind instruments with flute,picolo,pan flute etc. or sets of drums or the class of percussions including a grand piano, stage piano, console piano etc...

So its clear that each class has a very distinctive sound as an average person can easily classify them while within a given class it is more difficult for an average person to identify every member (but a person with a musical ear will be able to distinguish almost every type of instrument from each other, thats one of the primary differences between a musically able person from a not so one)

Coming to what makes all those instruments (both the classes and members)  sound different from each other is their physical structure and vibration mechanisms. These are specific to almost every class and to each member based on their material, geometry and driving mechanisms so that the resulting acoustical waves will have specific harmonic characteristics for each instrument (indeed this is true for any object of the nature) which will lead to distinguishable timbre / tone differences by humans.

Coming to the ability of simulating the acoustical charactheristics of an instrument by another one. There are already electronic analog/digital synthesizers which can mimic many of the instruments to some degree. There are also more sophisticated acoustical modeling programs which offer the complete re-creation of the acoustical properties of any imaginable or existing object/instrument. I'm not experienced in their results so you can look further.


Radar signal processing - Difference between range gate/range bin In learning about radar signal processing, I see two different terms that sometimes seem to refer to the same thing, but sometimes different things:

Range gate
Range bin

Usually, I interpret both to simply mean "time sample", but sometimes I see range gating used to describe just taking a few samples around the estimated target range during a tracking phase.

Is there a more nuanced definition that I should be aware of?

If we wanted to add a nuance, the "gate" is the actual switching on of the receiver for the duration of observing the reflected signal (which can be as narrow as a single sample, which even a single sample for analog to digital conversions is an integrated observation over a sample time), while the range bin is which "gate" you are in. (Imagine a receiver that can have several gates covering the full possible range of the target). This originated from pulsed radar systems which is what my description above would apply to but can be shown to be equally applicable to other radar methods such as FMCW where an ambiguity in position can occur based on the repetition rate of the waveform.

For tracking, I believe this is where the term "Early-Late Gate" originated: In this use the gate is subdivided into an early half and late half (or can be described as separate gates, one slightly earlier and the other slightly later in vicinity of known target location), and the central portion of the gate position is adjusted earlier or later based on the difference between the early and late gate.  Here too we are "gating" the received signal to limit our observation to a narrow time interval (or even one sample).


How to calculate STFT of a function for a rectangular window How to calculate the STFT (by hand) of $$u(n)\cos(0.2\pi n)$$ for a rectangular window of a length 20, positioned at $n = 5$.

I know that to use STFT I need to divide longer signal to a shorter parts and than calculate Fourier Transform on each part. Also doing it from the definition is very long and I assume that there is a more efficient way to do it on paper.

I don't know how to start with the question so any materials on this topic are appreciated.

I came up with this idea:

First picture is just a rectangular window of length 20 and position 5.
Second picture is how I see the STFT of it. Now to calculate STFT should I provide some coefficients like mainlobe width and highest sidelobe? How I can calculate them from my data?

Are my pictures good?

I don't understand what you mean by "positioned at $n=5$" with length $20$. Indeed:

is "positioned" the beginning  or the center?
is length an half-length around a central position, or a full length?
finally, at a given location, this is not an STFT anymore, but a mere windowed Fourier transform.

If I interpret your question in its most obvious sense, the window starts from $n=5$ to $n=20+5-1=24$, on an interval where $u[n]=1$, so you'd just have to compute a simple DFT of a cosine. The result can be obtained via Euler/De Moivre formulae, with two finite sums of geometric series.

If this is not the case, please provide more information, close forms are at hand


What is the difference between PSD and squared magnitude of frequency spectrum? The power spectrum of a signal can be calculated by taking the magnitude squared of its Fourier transform. Being an audio person, the signal of interest for me would be a time series.

How does this representation differ from a PSD (power spectral density), and importantly, in what practical situations should one use a PSD instead of the power spectrum described above?

The power spectral density describes the density of power in a stationary random process $X(t)$ per unit of frequency. By the Wiener-Khinchin theorem, it can be calculated as follows for a wide-sense stationary random process:

$$
S_{xx}(f) = \int_{-\infty}^{\infty} r_{xx}(\tau) e^{-j2\pi f \tau} d\tau
$$

where $r_{xx}(\tau)$ is the autocorrelation function of the process $X(t)$:

$$
r_{xx}(\tau) = \mathbb{E}\left(X(t)X(t - \tau)\right)
$$

This is only valid for a wide-sense stationary process because its autocorrelation function is only a function of the time lag $\tau$ and not the absolute time $t$; stated differently, this means that its second-order statistics don't change as a function of time.

With that said, if you have a sufficiently-detailed and accurate statistical model for your signal, then you can calculate its power spectral density using the relationship above. As an example, this can be used to calculate the power spectral density of communications signals, given the statistics of the information symbols carried by the signal and any pulse shaping employed during transmission.

In most practical situations, this level of information is not available, however, and one must resort to estimating a given signal's power spectral density. One very straightforward approach is to take the squared magnitude of its Fourier transform (or, perhaps, the squared magnitude of several short-time Fourier transforms and average them) as the estimate of the PSD. However, assuming that the signal you're observing contains some stochastic component (which is often the case), this is again just an estimate of what the true underlying PSD is based upon a single realization (i.e. a single observation) of the random process. Whether the power spectrum that you calculate bears any meaningful resemblance to the actual PSD of the process is situation-dependent.

As this previous post notes, there are many methods for PSD estimation; which is most suitable depends upon the character of the random process, any a priori information that you might have, and what features of the signal you're most interested in.


How to remove heart beats interference from pectoralis major electromyogram? I am using electromyography to detect activity on pectoralis major, however I found an interference from heart beats that is affecting my posterior analysis. The sampling rate is 2000 Hz, and the time of sampling is 150 s. The main signal is the muscle activity as the arm moves (four major big masses in the figure, between $-2*10^{-4}$ and $1*10^{-4}$), and the interfercence can be seen in the figure below as pulses at constant rate on values between $-5*10^{-5}$ and $1*10^{-5}$:

For comparison, I provide a figure of another muscle's (Deltoid) signal , without the heart interference:

As I don't know much about filters, I would like to ask for advice on how to remove the heart interference from this signal. Specifically, some directions on how to do it on matlab or R.

I can not use methods to record simultaneously heart and pectoralis major signals, since all individuals have already been sampled.

I found some papers addressing this issue, but I'm not able to implement the suggested solutions. For example:

Removing ECG noise from surface EMG signals using adaptive filtering

A sample file with this data is available here.

One possible solution is to use a notch filter to remove frequencies commonly associated with heartbeats somewhere between 40 to 100 bpm or .66-1.66Hz

Or if you dont need any low frequencies just use a "highpass" filter where the cutoff is around 2Hz

These solutions assume that anything near the frequency of a heartbeat is unimportant

edit

using the sample data you provided I applied a simple HPF with a cutoff at 2Hz (technically it should be a little past 2, but you will see even as it is it reduced the low frequencies by several orders of magnitude). This may have removed a significant amount of the signal power but it didn't perfectly replicate the cardiac noise. I suppose you have to think the electrical pulse associated with contracting the separate chambers of the heart are much more complex than the resultant heartbeat. In any even here is the data and code, just so you can see the process

x = importdata('pectoralis Major sample.txt');
samples_per_second = 2000;
Ts = 1 / samples_per_second;

%filter design
cuttoff_hz = 2;
cutoff_freq = 2*pi*cuttoff_hz;
tau = 1 / (2*pi*cutoff_freq);
a = Ts / tau;

length_x = length(x);
nfft = 2^nextpow2(length_x);
X=fft(x,nfft)/length_x;

%single pole hpf
x_hpf = filter([1-a a-1],[1 a-1], x);
X_HPF = fft(x_hpf,nfft)/length_x;

figure(1);
subplot(2,2,1);plot(2*abs(X(1:nfft/2+1)));title('original fft')
subplot(2,2,2);plot(2*abs(X(1:nfft/2+1)));title('original fft')
xlim([0 5])
subplot(2,2,3);plot(2*abs(X_HPF(1:nfft/2+1)));title('filtered fft 0-5Hz')
subplot(2,2,4);plot(2*abs(X_HPF(1:nfft/2+1)));title('filtered fft 0-5Hz')
xlim([0 5])

figure(2);
subplot(2,1,1);plot(x);title('original samples')
subplot(2,1,2);plot(x_hpf);title('filtered samples')


Convolving Room Impulse Response with a Wav File (python) I have written the following code which is supposed to put echo over an available sound file. Unfortunately the output is a very noisy result which I don't exactly understand. Can anybody help me with regard to this? Is there any skipped step?

#convolving a room impulse response function with a sound sample both of stereo type
from scipy.io import wavfile
inp = wavfile.read(sound_path + sound_file_name)
IR = wavfile.read(IR_path + IR_file_name)
if inp[0] != IR[0]:
    print "Size mismatch"
    sys.exit(-1)
else:
    rate = inp[0]
print sound_file_name
out_0 = fftconvolve(inp[1][:, 1], IR[1][:, 1])
out_1 = fftconvolve(inp[1][:, 1], IR[1][:, 1])
in_counter += 1
out = np.vstack((out_0, out_1)).T
wavfile.write(sound_path + sound_file_name + '_echoed.wav', rate, out)

The way you read the files is a little odd.  The method used here is clearer with respect to the sampling rate and the signal.
You will note that the above example calls pcm2float after reading.  This may be the cause of your problem.  The wave reader returns an array of ints.  Processing those (quite large) values may well cause problems.  Try converting your wave data to floats before doing the convolution.


Determining the causality of a signal with it's pole-zero plot I have the following question:

  Pole-zero plot of x(t) and y(t) are given below:

  The signal $g(t)$ and $h(t)$ are defined as $g(t)=x(t)e^{-3t}$ and
  $h(t)=y(t)*e^{-t}u(t)$. If $g(t)$ and $h(t)$ are both absolutely
  integrable, determine whether the signals $g(t)$, $h(t)$  are
  left-sided/right-sided.

My try:

I take the laplace transform of both the signals and get $G(s)=X(s+3)$ and $H(s)=Y(s)\cdot \frac{1}{s+1}$

Also because both $g(t)$ and $h(t)$ are absolutely integrable, their transforms must be stable so both must have $jw$-axis in their respective ROC.

As we can see $X(s+3)$ shifts the pole-zero plot to the left by $3$ units so we have all the poles in the left $s$-plane and ROC would be $Re\{s\}>-1$ hence $g(t)$ is right sided.

Similarly $H(s)$ has all the poles in the left $s$-plane and ROC is again $Re\{s\}>-1$ hence $h(t)$ is right sided.

Is this reasoning correct?

Looking at the pole-zero plots of the continuous-time signals $x(t)$ and $y(t)$, and the new signals $g(t)= x(t)e^{-3t}$ and $h(t) = y(t) \star e^{-t}u(t)$, the pole locaitons of $g(t)$ and $h(t)$ are found to be:

$$g(t) = x(t)e^{-3t} \implies G(s) = X(s+3) \implies \text{ Re(poles) } = \{-1,-1 \} $$

$$h(t) = y(t) \star e^{-t}u(t) \implies H(s) = Y(S) \frac{1}{s+1} \implies \text{ Re(poles) } = \{-2,-2,-1 \} $$

From these pole locations we see the following:

1-$g(t)$ has two possible ROCs: $Re\{s\} <-1$, and $Re\{s\} > -1$, and only the second one includes the $j\omega$ axis and hence can be stable.

2-$h(t)$ has three possible ROCs: $Re\{s\}<-2$, $-2 < Re\{s\} <-1$, and $Re\{s\} > -1$, and only the last one includes the $j\omega$ axis and is stable.

So for $h(t)$ and $g(t)$ to be absolutely integrable (stable), their ROC's must include $j\omega$ axis and this means their ROCs are to the right of the largest poles which implies that the signals are causal right sided signals.