Difference between discrete time fourier transform and discrete fourier transform I have read many articles about DTFT and DFT but am not able to discern the difference between the two except for a few visible things like DTFT goes till infinity while DFT is only till N-1. Can anyone please explain the difference and when to use what? Wiki says

  The DFT differs from the discrete-time Fourier transform (DTFT) in
  that its input and output sequences are both finite; it is therefore
  said to be the Fourier analysis of finite-domain (or periodic)
  discrete-time functions.

Is it the only difference?

Edit:
This article nicely explains the difference

alright, i'm gonna answer this with an argument that "opponents" to my rigid nazi-like position regarding the DFT have.
first of all, my rigid, nazi-like position:  the DFT and Discrete Fourier Series is one-and-the-same.  the DFT maps one infinite and periodic sequence, $x[n]$ with period $N$ in the "time" domain to another infinite and periodic sequence, $X[k]$, again with period $N$, in the "frequency" domain.  and the iDFT maps it back.  and they're "bijective" or "invertible" or "one-to-one".
DFT:
$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] e^{-j 2 \pi nk/N} $$
iDFT:
$$ x[n] = \frac{1}{N} \sum\limits_{k=0}^{N-1} X[k] e^{j 2 \pi nk/N} $$
that is most fundamentally what the DFT is.  it is inherently a periodic or circular thing.
but the periodicity deniers like to say this about the DFT.  it is true, it just doesn't change any of the above.
so, suppose you had a finite-length sequence $x[n]$ of length $N$ and, instead of periodically extending it (which is what the DFT inherently does), you append this finite-length sequence with zeros infinitely on both left and right.  so
$$ \hat{x}[n] \triangleq \begin{cases}
x[n] \qquad & \text{for } 0 \le n \le N-1 \\
\\
0 & \text{otherwise}
\end{cases} $$
now, this non-repeating infinite sequence does have a DTFT:
DTFT:
$$ \hat{X}\left(e^{j\omega}\right) = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] e^{-j \omega n} $$
$\hat{X}\left(e^{j\omega}\right)$ is the Z-transform of $\hat{x}[n]$ evaluated on the unit circle $z=e^{j\omega}$ for infinitely many real values of $\omega$.
now, if you were to sample that DTFT $\hat{X}\left(e^{j\omega}\right)$ at $N$ equally spaced points on the unit circle, with one point at $z=e^{j\omega}=1$, you would get
$$ \begin{align}
\hat{X}\left(e^{j\omega}\right)\Bigg|_{\omega = 2 \pi\frac{k}{N}} & = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] e^{-j \omega n} \Bigg|_{\omega = 2 \pi\frac{k}{N}} \\
& = \sum\limits_{n=-\infty}^{+\infty} \hat{x}[n] e^{-j 2 \pi k n/N}  \\
& = \sum\limits_{n=0}^{N-1} \hat{x}[n] e^{-j 2 \pi k n/N}  \\
& = \sum\limits_{n=0}^{N-1} x[n] e^{-j 2 \pi k n/N}  \\
& = X[k]  \\
 \end{align} $$
that is precisely how the DFT and DTFT are related.  sampling the DTFT at uniform intervals in the "frequency" domain causes, in the "time" domain, the original sequence $\hat{x}[n]$ to be repeated and shifted by all multiples of $N$ and overlap-added.  that's what uniform sampling in one domain causes in the other domain.  but, since $\hat{x}[n]$ is hypothesized to be $0$ outside of the interval $0 \le n \le N-1$, that overlap-adding does nothing.  it just periodically extends the non-zero part of $\hat{x}[n]$, our original finite-length sequence, $x[n]$.


Looking for algorithm to extract line network from binary images within LabView I would like to create a LabView VI that receives an input image and the RGB-color codes of the powerlines depicted in the image. The output should (eventually) consist of id'd polygonial connections and their coordinates.

Right now i've finished creating a VI that converts the source image to a binary image with the respective line color types extracted into a binary image (see below).

I've been trying to apply some of edge/line finding VIs supplied by LabView Vision, but the results have not been very convincing.

I have considered implementing some sort of vectorization algorithm like 'Straight Line Extraction Using Iterative Total Least Squares Methods' or even some sort of line-walking-and-mapping algorithm, but before I do so I'd like to know if there already is some kind of fitting algorithm out there.

Attached you can find examples for the source and binary images.

Thanks in advance for your input!

The go-to algorithm for extracting edges is to transform the image the an appropriate Hough space. Luckily, the Hough transformation for extracting straight lines is one of the best documented.

It looks like these folks may have a Hough transformation implementation for LabVIEW you could use. If not, it's a pretty well documented algorithm and you could try to roll your own or try to integrate an implementation from another language, eg: MatLab, OpenCV.

With the appropriate Hough transformation, you would be able to extract line segments from your image (converting to a binary image is a great start). These line segments may not connect perfectly, depending on how straight your line is and how clear your image is, and you may need to process the vectors a bit more to get the results you desire.


Ramp function as derivative in frequency domain? It is said that to get Laplacian of Gaussian in frequency domain, we may multiply the Fourier transform of Gaussian with two differentiating ramp function (1 ramp gives 1 order of derivative).

The description from the material that I was following:

And the file can be found here

So how could we possibly get derivatives by multiplying ramp functions in the frequency domain? How does the math work here?

You may be familiar with the classical properties of the Fourier transform. One of them is the one regarding differentiation in the time domain:

$$x(t)\xrightarrow{\mathscr{F}}X(j\omega) \implies \frac{d}{dt}x(t)\xrightarrow{\mathscr{F}}j\omega X(j\omega)$$

And there appears the ramp function, $j\omega$.

The demonstration of why this is true is straightforward. If we have a signal $x(t)$ with Fourier transform $X(j\omega)$:

$$x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(j\omega)e^{j\omega t} \ \mathrm{d}\omega$$

Then by differentiating both sides:

$$\frac{d}{dt}x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}j\omega X(j\omega)e^{j\omega t} \ \mathrm{d}\omega$$

And that's the (pretty simple) proof.


how to create a harmonic mask from fundamental frequency? Using an algorithm that gives the fundamental frequency across time, example https://www.mathworks.com/help/audio/ref/pitch.html (or some other one?)

how would I then use that, to create a harmonic mask of nHarmonics based on that fundamental frequency? So the mask will contain the f0 + nHarmonics harmonics.

The point of this mask will be to then use it as a filter for the audio mixture.

I want to do this using MATLAB.

edit: I forgot to mention that the f0 is not a constant value, it is a ranged f0 track from fmin to fmax, for example dominant f0 from 100Hz to 800Hz. So the f0 track contains different values, not just one value. I can understand a comb filter working for one value, but this is a bit more complex, as the f0 track has multiple different values, so now each one of those would have to have nHarmonics multiples above the f0 for me to create this harmonic mask.

See this post which describes creating a notch filter using a delay and add:

Moving Average for Notch filtering

The filter used is generically given as

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

This results in passing DC and placing zeros at the normalized radian frequency locations given by:

$$\omega_n = \frac{2\pi n}{N} + \frac{\pi}{N}$$

for all integers n and $\omega_n$ in the range of 0 to $2\pi$ to correspond to the primary frequency range of DC to the sampling rate.

This comes directly out of factoring the polynomial given by H(z), resulting in zeros on the complex unit circle at the angular frequencies given above. Because of the extra $\pi/N$ rotation on the unit circle, this filter works well for filtering odd harmonics based on prudent choice of N and therefore the sampling rate as an integer multiple of the notch location.

If instead a delay and subtraction is used:

$$H(z) = 1-z^{-N}$$

The result would be a zero (notch) at DC and and placing zeros at the normalized radian frequency locations given by:

$$\omega_n = \frac{2\pi n}{N}$$

Which comes from factoring $H(z) = 1-z^{-N}$ which results in the Nth roots of unity (zeros evenly spaced around the unit circle). This approach conveniently allows for placing a notch at every harmonic of a selected frequency.  To introduce the more advanced approaches given below, it is useful at this point to visually plot the zeros in the z domain for this case, using the example $H(z) = 1- z^{-8}$, along with the associated frequency response:

Both of these approaches are referred to as "Comb Filters" and notice in this case how the notch occurs at a fundamental frequency and every integer multiple of that frequency. The requirement here is the sampling rate be an integer multiple of the notch frequency. The filter has linear phase which is important if phase distortion is of concern, but a potential drawback is the gradual loss in frequency components at other than the distinct harmonics.

An improved approach for the amplitude response would be to extend the concept of the second order notch filter which I have derived at this link given below.

Transfer function of second order notch filter

This will result in much flatter amplitude response over the pass frequencies, but the tight notches come at the cost of increased phase distortion in proximity of the notches which may be an issue depending on the application.

To extend this to multiple notches at integer multiples note that the background of the 2nd order notch filter is the simpler DC notch filter given by

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

Which has a zero at $z=1$ and a pole placed inside the unit circle at $z = a$ in close proximity to that zero; the closer the placement the tighter the bandwidth but also the more precision (bit size) is required to achieve the desired performance and maintain stability. Notice below the comparison to a single zero as done in the Comb Filter implementations (with $N=1$) to an example case of the DC notch filter ($a=0.95$ for the example shown):

At this point you may see where I am going with this. We can combine the concept of the DC Nulling filter with the comb filter to achieve tight notches at integer harmonic locations by placing poles in close proximity to the zeros given in the Comb Filter example as depicted in the graphic below.

Following similar polynomial factoring relationships into the complex plane, this is realized with the following transfer function below:

$$H(z) = K\frac{1-z^{-N}}{1-a^Nz^{-N}}$$

The amplitude normalization factor K if desired is determined by solving for the value in the passpand $H(z=\omega_n)$ where $\omega_n = \frac{2\pi n}{N} + \frac{\pi}{N}$ and inverting this, which ends up being:

$$K = \frac{1+a^8}{2}$$

The result for N= 8 and a = 0.95 is shown in the plot below.

The possible challenges with this approach is the significant passband phase distortion as N increases, and the amount of precision required in implementation.

If the phase distortion was an issue and the notches need not be as tight as the previous IIR implementation can provide, then a third option is to combine the linear phase notch response given by $H(z) = 1-z^{-N}$ cascaded with a simple 3-tap linear phase compensator similar to what I detailed in this link (the compensator shown is actually a raised minus cosine response; a raised positive cosine response would be the approach for this case, following the same process outlined at that post.)

how to make CIC compensation filter

In particular see the plot in that post for the interpolated CIC compensator which is implmented by simple zero insert (meaning to interpolate by 8, simply insert 7 zeros between each filter coefficient in the 3 tap compensator). This won't achieve the flatness and sharp nulls as done with the IIR approach above but will have significantly improved flatness over a wider passband range than the simple notch filter approach AND with linear phase. The cascade of the simple delay and subtract comb filter with the interpolated 3 tap FIR compensator (given by the two responses shown below) results in a significantly flatter passband response with minimal added complexity.


FFT- based frequency offset estimation algorithm I am trying to understand and implement FFT based recovery algorithm suggested me in the post and it was also discussed here
Let assume rx_signal -signal at the receiver'input with frequency offset. according to the posts i do the following steps in matlab:
% fft

fft_rx = fft(rx_signal)

% [peak p] = max(abs(fft_rx))

I don't understand what my next step is?

Given the links in the question, the OP is interested in more details on implementing a squaring algorithm for carrier recovery. A squaring algorithm is a popular analog option for the recovery of BPSK and QPSK modulations based on its simplicity and is completely applicable to digital implementations as well, but it is not the only algorithm nor most efficient. Specific to squaring algorithm for carrier recovery, the following will help understand how this approach works and its limitations.
If the modulation is BPSK, then the symbols each have a carrier phase of either 0 or 180 degrees relative to the unmodulated carrier (which the OP is trying to recovery). A squaring loop strips the modulation by using the property:
$$\cos(\alpha)cos(\beta) = \frac{1}{2}cos(\alpha - \beta) + \frac{1}{2}cos(\alpha + \beta)$$
From above observe what happens if we square a BPSK signal given as
$$f(t) = \cos(2\pi f t+\phi[m])$$
Where $\phi[m] = 0, \pi$ for each symbol $m$
$$f^2(t) = \frac{1}{2}\cos(0) +\frac{1}{2}cos(4 \pi ft + 2\phi[m])$$
The first term is a DC offset value and is ignored or filtered out. Notice what happens in the second term from which we can now derive a constant carrier signal: The frequency has doubled (so squaring a tone in general will double its frequency, and frequency doublers are implemented by doing this) AND importantly the modulation has been completely removed, since doubling a signal's frequency will also double its phase. Since the phase was going between $0$ and $\pi$ (180 degrees) before being doubled, it will go between $0$ and $2\pi = 0$ after being doubled.  This is a common analog recovery scheme due to the simplicity of implementing a frequency doubler (using a mixer, or even an XOR gate), and a frequency divider (with a flip-flop) to divide the recovered tone without the modulation by 2 to recover the coherent BPSK carrier.  Pulse shaping will introduce other spectral artifacts, but they will all be sufficiently lower than the strong recovered carrier such that a PLL can be used to clean-up and provide a clean recovered coherent carrier for demodulation.
The same approach can be used for QPSK, but in this case we have four phase states, so the doubling needs to occur twice (resulting in $f^4(t)$).  Doubling once will reduce the QPSK signal to a BPSK signal at double the carrier frequency, and doubling that will then strip the remaining BPSK modulation for a recovered carrier at four times the original frequency. (We also see that for any four quadrature phase states, if we raise it to the fourth power it will reduce to one phase state, meaning the modulation will be stripped: for example, for all phases $(0, \pi/2, \pi, 3\pi/2)$, any of these phases multiplied by 4 will be modulo $2\pi$ which is a phase equal to 0.
This approach could be done using an FFT such as first multiplying the QPSK signal by four and then use the FFT to find the strongest carrier that is at approximately four times the original carrier frequency, but that is highly inefficient and an approach I would not recommend following (this is basically breaking fred harris' rule "Don't copy the analog!" given we have so much more ability to provide direct optimized solutions with high efficiency with more direct DSP approaches targeting the core problem, where analog implementations are often much more constrained by the physical components themselves). Further it immediately requires a sampling rate 4 times higher to handle the 4x carrier increase (even for any given carrier offset). Further once the strongest FFT bin is found, the actual frequency will inevitably be at an smaller offset between the bins, so further resolution will be required to maintain actual coherence. I suggest alternate carrier recovery approaches at the links below, but if it is necessary to pursue a squaring approach, this would typically be done in a loop such that a PLL is used to clean-up the strongest single tone signal (at double the actual carrier for BPSK or four times the actual carrier for QPSK).
For a suggested carrier recovery approach please also see these other links:
High modulation index PSK - carrier recovery
BPSK and QPSK demodulation computational complexity
Demodulation of 4-QAM


Removing DC offset and ensuring zero-sum for 2D filter I'm working on project where I'm using Independent Components Analysis (ICA) on NxN patches from images to separate components in an image signal. Each column of the resulting MxN^2 unmixing matrix is reshaped back to an NxN filter containing the independent components. For my project, this filter is then convolved with new images for feature extraction purposes.

My question is, how do I remove the DC offset of these NxN filters that result from ICA? In order words how can I ensure that (in MATLAB parlance) sum(sum(filter)) == 0 returns 1? Is it fair to just subtract the filter's sum from each component, or will that affect the filter's response?

If the filter sums to some non-zero value X, is it appropriate to just subtract X/number_of_coefficients from each filter coefficient?

Removing the mean (and not the sum) from a linear filter leads to a filter that removes DC.

For example:

N = 10;
h = [0.5 1 0.5];
conv(ones(N,1), h, 'valid')

leads to 2s.

N=10
h = [0.5 1 0.5]; h_high=h-mean(h);
conv(ones(N,1), h_high, 'valid')

leads to zeros.

Alternatively, you can use any high pass filter.


Auto-correlation function, an inverse problem $x[n]$ is a complex function $n=0,1,2,\cdots,L-1 $

we assume $x[n]$ is periodic in its index:  $x[n+L]=x[n]$

Its auto-correlation function $C[n]$ is uniquely defined as:
$$
C[n]=\sum_{i=0}^{L-1} x[i+n]x^*[i]
$$
$C[n]$ also has the periodic property: $$C[n+L]=C[n]\tag{1}$$

And ''conjugate-symmetry'' property: $C[-n]=C^*[n] \tag{2}$

Now my question is:

For given $C[n]$, which satisfies property (1) and (2):

Can we find the corresponding $x[n]$ ?

If yes, is it unique? and what is the method to find $x[n]$?

If no, what other constraint properties should we add to $C[n]$, in order to make it yes?

Let's look at the case $x[n] \in \mathbb{R}$, where $x[n]$ is real.

Autocorrelation is basically convolution of the signal with it's time inverse. This can be easily expressed in the frequency domain.

$$ \mathscr{F}\Big\{ r_{xx}[n] \Big\} = \mathscr{F}\Big\{ x[n]  \Big\} \cdot  \mathscr{F}\Big\{ x[-n] \Big\} $$

$$R_{xx}(\omega) = X(\omega)\cdot X^*(\omega) = \Big| X(\omega) \Big|^2
$$

So it's easy to see that the Fourier Transform of the auto correlation is simply the magnitude squared of the Fourier Transform of the input signal. That's sometimes referred to as the Power Spectrum.

It's also easy to see that information gets lost in the process. There are $N$ unique values going in but because of the symmetry properties of the auto correlation there are only $\frac{N}{2}$ unique (independent) values coming out. Looking in the frequency domain, we can see that the phase is lost.

  If yes, is it unique? and what is the method to find $x[n]$?

No, it's not unique

  Can we find the corresponding $x[n]$ ?

There is an infinite number of $x[n]$.

Take the Fourier Transform of the autocorrelation
Take the square root
Add an aribtiraty (but odd-symmetric) phase function
Do an inverse Fourier Transform

Any signal derived this way will have the same original auto correlation function.

  if no, what other constraint properties should we add to $C[n]$, in order to make it yes?

You can't make it a yes, since it's not unique. No matter what auto correlation you choose, there will be infinite $x[n]$ that will have it as an autocorrelation.


Signal deconvolution from exponential and linear tail I would like to deconvolve the signal shown in the plot below from its linearly decaying tail and the exponential levelling. That means that after the deconvolutin, the tail should be at a constant value of ~410 mV. Does anyone have an idea how to approach this problem?
I uploaded the data here:
https://nextcloud.mpp.mpg.de/nextcloud/index.php/s/RwXEzEa6bxHTByR

Thank you very much in advance!

You probably can't. If I interpret your description correctly, this is the unit step response of your system. This looks very much like the unit step response of a high pass filter, that means the system is filtering out exactly the frequency that you want (which is DC or 0 Hz).
You may be better off creating a parametric model of your system response and try to do a weighted least square error fit of the data.


Filtering performance on Poisson noise with quadratic data-fidelity We recently performed a work on signal filtering/component separation (sparse signal/trend/noise). The cost function contains:

A quadratic data fidelity term,
Some smoothed $\ell_1$ terms for sparsity and positivity promotion.

This was initially designed for Gaussian noise, and performed satisfactorily. Reviewers asked to test the algorithm on the same signals with Poisson noise, and it performed well also.

In many works concerned with Poisson-Gaussian noise mixtures, authors usually use more involved penalty functions and/or use variance stabilization transforms (such as Anscombe's).

Presenting the work to different signal processing folks, several colleagues (signal and image persons equally) said they were "not surprised" that a simple quadratic term would work well with Poisson noise too.

Since I am less a theoretician than a frequentist, and a bayesian (much) less than a frequentist, I do not understand the reason behind the relatively good behavior with both Gaussian and Poisson perturbations.

May a reader offer practical or theorical hints behind this "non surprise"?

For large intensities / large "bins", i.e. "areas for which events are counted and accumulated", Poisson processes lead to nearly Gaussian distributed individual values -- basically, without trying to derive this, I think that's  the application of the CLT on a lot of realization of a point process.

EDIT: Shameless plug: I really really like the wikipedia page on "Relationships among probability distributions; notice the arrow between the Poisson and the normal distribution:


Low-Pass Filtering on Small Chunks of Audio I've got a script (Python) that chops 44.1 kHz realtime audio into 10 ms chunks, processes each chunk, and stitches the chunks back together with no added artifacts.  I'm trying to create a low-pass filter that operates on each chunk, but every filter I've tried introduces unwanted artifacts.

I've tried creating various types filters, including IIR (Butterworth) and FIR (Hann, Hamming, rectangular windows, etc.), and Remez.  The attempted methods of applying the filters to the input signal include just a straight transfer function, convolution, and a forward-then-backward filter (namely scipy.signal.filtfilt).  I've also tried the moving average approach, but I don't know what to do near the endpoints of the chunk.

Is there such a filter that is designed for this kind of application?

INPUT WAVE (lots of harmonic content)

|<-- 1 CHUNK  -->|
:                :                :                :
:                :                :                :
: -----       ---:-       -----   :   -----       -:
:|     |     |   : |     |     |  :  |     |     | :
:|     |     |   : |     |     |  :  |     |     | :
:|     |     |   : |     |     |  :  |     |     | :
:       -----    :  -----       --:--       -----  :
:                :                :                :
        V                V                 V
        V                V                 V
     FILTER            FILTER            FILTER
        V                V                 V
        V                V                 V
:                :                :                :
:   ---         -:-         ---   :     ---        :
:  /   \       / : \       /   \  :    /   \       :
: /     \     /  :  \     /     \ :   /     \     /:
:/       \   /   :   \   /       \:  /       \   / :
:         ---    :    ---         :--         ---  :
:                :                :                :

OUTPUT WAVE (less harmonic content)

EDIT: Following the answer from @hotpaw2, I wish to try the IIR approach first (I'm a novice at this, but my gut tells me try IIR first).  For discussion purposes, let's define these arrays:

B = transfer function numerator coefficients
A = transfer function denominator coefficients
I0 = previous input array (same size)
I1 = current input array (441 samples)
R0 = previous output, post-filtered (DO WE NEED?)
R1 = current output, desired post-filtered result array

Do I generate the initial conditions from I0 only, or do I need to figure R0 into the equation?  To ask another way, here's the Python code I'm trying:

initCon = scipy.signal.filtic(B, A, I0)

Once I have this set of initial conditions, are they directly applied to the filter, or is there more that I need to do with them?  Again, the Python code:

R1, _ = scipy.signal.lfilter(B, A, I1, z = initCon)

As the code stands now, it executes but there is a lot of off-frequency buzzing.

EDIT 2:  Just thought I'd post the output of the filter with a sine wave input.  It's quite obvious something goes wrong every 441 samples (the beginning of each chunk):

Filtering a data block usually results in more data than fits in the size of the block.  If you throw away this added data, that will produce artifacts across blocks.

For FIR filters, you need to pad each chunk or block with at least the length of the impulse response of your filter before filtering (>= N+M-1).  Then use overlap-add or overlap-save (FFT fast convolution algorithms) to process all the extra results from convolution of the added padding (usually carry them over to the next chunk or chunks).  Otherwise you are either getting circular convolution artifacts in each block, or are throwing away part of the filter response from applying your filter to each block.

For IIR filters, you need to carry over the filters internal state between chunks, and not start over from zero or other default initial state of the IIR.  That state still carries energy from the filter result, and if you throw it away, you will get glitches between chunks.


Is cutting EQ scientifically different than boosting EQ? This is a fairly standard belief in the audio production scene that cutting EQ is somehow "cleaner" than boosting EQ. You hear arguments like there being different phase characteristics or whatever. The more careful arguers refer to psychoacoustics (that it's easier to listen to the effects of cuts, rather than boosts).

However, is there a difference in the filter in terms of cutting and boosting? And which one is "cleaner"?

There is a very compelling reason for preferring cutting over boosting. Assuming you are using the typical bell-response parametric EQ bands, then you will see a significant phase distortion around the resonant frequency, i.e. the frequency you're boosting or attenuating. (There is an exception to this rule: Linear phase EQs that are often used for mastering do not show this phenomenon. Those only exist as digital implementations however and have other problems.)

The most noticable effect of the phase distortion is that it comes with a certain group delay around the resonant frequency, and especially narrow changes can result in audible post-ringing and temporal distortion of transients.

Now if you're boosting what you want to hear more, you introduce this kind of problem exactly in the band that you want to give more prominence. Cutting also introduces the same issues, but you usually cut what you don't want to hear as prominently. So cutting will move the phase artefacts to the sonically less important frequency range while boosting will place it right where you don't want it.

That's the theory at least. In practice, the rule can be relaxed in most cases. The negative effects of resonant phase distortion greatly depend on the derivative of the phase response with respect to the frequency. That means that resonances with very small bandwidth will result in greater issues than wide band resonances. And this is absolute bandwidth, not relative bandwidth. This means for a typical EQ band that has a bandwidth proportional to the center frequency ("constant-Q filter"), the bandwidth automatically gets larger for higher frequencies.

So a bell EQ at a few kHz will have a large bandwidth, even if it has a rather high resonance. Therefore cutting vs boosting doesn't make a huge difference for higher frequencies. It does however get relevant for frequencies in the lower mid and the bass range, especially when you're using high Q-factors.


Covariance matrix explanation I am trying to understand and visualize the concept of a covariance matrix. Suppose I have a matrix:

$A = \begin{pmatrix} 2 & 3 & 4\\ 5 & 5 & 6 \end{pmatrix}$

how do I calculate its covariance matrix, and what is its practical meaning?
(All I was able to understand by now is that on the diagonal of the covariance matrix, variances for particular variables are placed, and on the upper and lower fields are the correlations between those variables.)

You can't really have a covariance of a matrix.  What you can have, is a covariance matrix of a set of vectors.  So, if you think of the rows of your matrix A as two vectors in 3D: [2 3 4] and [5 5 6], then the covariance matrix of this set of two vectors is C = A' * A (A transpose times A).  Note that if you shuffle the rows of A in a different order, C does not change.

Now that you have C, you understood it correctly.  You can also do interesting things with the eigenvalues and eigenvectors of C.

Oops... Totally forgot.  Before you do anything, you have to compute the mean vector (column-wise mean of A) and subtract it from every row.


How to determine SNR required for adding AWGN noise to clean IQ samples? I'm doing an experiment w.r.t PSD output. I have generated some 11a  IQ samples and checking its PSD output in MATLAB using pwelch and comparing with the PSD output obtained using hardware.

Here are the block diagram  and output snapshots

wfm from VSG is transmitted at -10 dBm

I realized that VSG---> loopback cable ---> VSA   adds noise to the samples which results in noise floor whose level is at -90ish dBm/Hz . And when I was using pwelch function , noise floor goes down to -140ish dBm/z

To Make it more similar , I am adding some AWGn noise before I input to the pwlch function using MATLAB

iqSamples = awgn(iqSamples,SNR,'measured')

I was doing some trial and error for SNR value and using SNR =50 dB gave me the following output

Few questions here

1)  In this case ,On what parameters does SNR depends?
2)  I' guessing this value might change with hardware and also loopback cable ? But what if I don't change the hardware , loopback cable used and VSG tx_power_level ?

Thanks for going through my simple question. Still trying to learn communications/DSP.

Note that the MATLAB result is floating point (or may be fixed point but with much higher precision) while the VSA hardware is fixed point with lower precision that also has dynamic range limitations in the analog portion of the hardware as well. To see this directly and more clearly, try these two experiments:

Create a single carrier tone and capture results at different signal levels of that carrier. Ensure the sampling rate is much higher than the tone frequency (such as > 12x higher to capture the 3rd and 5th harmonic).

Observe the output with a properly windowed FFT (I recommend the Kaiser window) and increase the capture duration until the higher harmonics are clearly visible. (increasing the duration will lower the FFT noise floor).

From such testing you can determine the SNR (based on the wideband noise floor) as well as spurious free dynamic range (SFDR) of your test set-up. The SNR is estimated as the difference between the tone and the integrated power that is not at the tone frequency (and outside the main lobe of the window used).  The SFDR is the difference in power between the tone and the highest spur which may or may not be at a harmonic of the signal (In some definitions of SFDR the harmonic components are not considered and that distortion is captured in Harmonic Distortion). In any event find the maximum tone power that maximizes SFDR, and this sets a mark for "Full Scale" from which the other metrics can be computed.

If you proceed with this, please post your results. What would be interesting too is a plot of Pin vs Pout for your tone, from which you could also determine where it starts to compress. The point where Pout is 1 dB lower than the otherwise linear relationship is known as the 1 dB compression point.

To be complete you could also do two tone testing, where you generate two sinusoids closely spaced (but wider than the width of the main lobe of your window in the frequency domain) and measure the two-tone third order intercept point. See Wikipedia for more details on that test:

https://en.wikipedia.org/wiki/Third-order_intercept_point

I would want to understand all these parameters in conjunction with the use of any test equipment such that I understand what is being introduced by the test equipment rather than in the waveform itself.

Finally when adding AWGN to the signal to simulate the SNR for a discrete system, the AWGN is a power spectral density that is flat across the entire sampling bandwidth. The SNR of interest would be the portion of the spectrum that is in-band---thus if you are sampling at a higher rate, the AWGN with the same variance would be spread across a wider bandwidth resulting in lower SNR given the portion of the noise that would occupy the same bandwidth of the signal. Any noise outside of the signal bandwidth can be filtered so is not part of the SNR computation.


Why taking derivative amplifies noise In image processing we may use derivatives to help us detect the edges. While at mean time, this operation would also make the data noisier. But why do we have this effect?

My intuition is that if we denote the image data as

$I_{actual}(x, y) = I_{ideal}(x, y) + I_{noise}(x, y)$.

Then apparently when we do $I_{actual}(x_2, y) - I_{actual}(x_1, y)$, the noises will cancel out with each other? But how come they end up adding up to each other?

Differentiation in one axis will amplify high-frequency components. Assuming noise at Fs/2 (i.e maximum frequency), the gain at this frequency will be 6 dB. Intuitively, if you have a sequence alternating between +1 and -1 i.e (1, -1, 1, - 1), Then the difference will be +2, - 2, +2, -2, a gain of 6 dB

Another reason is that by differentiating, you remove the low-frequency components, so your eyes see mostly the high-frequency components giving you the impression of having more noise.


Generating a time series given a transfer function I'm trying to work my way through a paper I found online, titled "Three Models of Wind-Gust Disturbances for the Analysis of Antenna Pointing Accuracy" by W. Gawronski, 2002.

http://ipnpr.jpl.nasa.gov/progress_report/42-149/149A.pdf

In the Appendix of the paper, Section I, the author outlines an approach to generate a time series.

  (1) Apply white noise to the ﬁlter, H, with the following transfer
  function: $$H_d = \dfrac{0.1584 \cdot z^3 −0.3765 \cdot z^2 + 0.2716 \cdot z −0.0534}{ z^4 −2.9951 \cdot z^3 +3.0893 \cdot z^2 −1.1930 \cdot z + 0.0988}$$ and obtain time series of wind $\Delta v(t)$.

How do you apply white noise to a filter? How do you use the transfer function to create a time series?

I think that procedure should be as follows:

1.Generate the white (Gaussian) noise  signal - randn function in MATLAB.

2.Having a transfer function of your filter divide it by $z^4$ (highest power in denominator) to get a proper representation for IIR filter:

$$H_d(z)=\dfrac{0.1584\cdot z^{-1} −0.3765\cdot z^{-2} + 0.2716\cdot z^{-3} −0.0534\cdot z^{-4}}{1 −2.9951\cdot z^{-1} + 3.0893\cdot z^{-2} −1.1930\cdot z^{-3} + 0.0988\cdot z^{-4}}$$

3.Now you can extract the coefficients according to equations that can be found here: IIR.

$$ H(z) = \dfrac{\sum_{i=0}^{P}b_i z^{-1}}{1+\sum_{j=1}^{Q}a_j z^{-j}}$$

This will give you following vectors:

a=[1 -2.9951 3.0893 -1.1930 0.0988] % denominator
b=[0 0.1584 -0.3765 0.2716 -0.0534] % numerator

4.These are coefficients of your filter that you can apply to the white noise signal. If you are using MATLAB/Octave, then use the filter function. Something like:

y=filter(b,a,x)

where x is your white noise signal generated previously.

The amplitude response of your filter is shown below - because you are going to apply it to white noise (having all frequencies equally probable), then spectrum of your filtered signal should be corresponding to amplitude response of your filter. There is a high gain in lower frequencies which makes somewhat sense (wind gusts are mostly energetic in lower range):

Judging from your article, now you only have to normalise the filtered white noise and scale it accordingly.
Good luck!


Hilbert transformer's phase in .AC simulation I don't understand this: the HT works as it should but in an .AC simulation, the phase shows up like this (wrapped):

Why doesn't it look like a constant -90deg phase, like in a differentiator?

[edit]

The thing that also confuses me is that a regular low-pass FIR has the same phase, so does that mean the FIR also adds constant -90deg? Obviously not, and I can't see the explanation. Can someone please shed some light?

You have a linear phase because you probably implemented a causal filter. This just means that you have the desired phase shift of $-\pi/2$ plus a constant delay. The additional delay, which corresponds to a linear phase term, usually causes no problems. It is necessary to make the filter causal and therefore implementable.

The difference between an FIR lowpass filter and a Hilbert transformer is the phase value at frequency zero: it should be zero for a lowpass filter, and $-\pi/2$ for a Hilbert transformer.


How to adjust receiver gains to avoid saturation and quantization noise to optimise post digital processing? I am using a Hack RF One as SDR and the software GNU Radio to receive and process signals emitted from a drone. Given that the drone is moving and is emitting signals from a wide zone, received signal strengths vary dramatically : if the signal is too large, it over-ranges the ADC input and if the signal is too small, it gets lost in the converter's quantization noise.

At the reception, I have three different gains : RF gain, IF gain and BB gain. RF gain can be set to 0 or 14dB, IF gain from 0 to 40dB by 8dB steps and BB gain from 0 to 62dB by 2dB steps. At the transmission, I keep the transmitting power constant.

For now, my algorithm to manage the different values of gains is the following : Every x seconds, I look at the received digitised signal peak-to-peak value $V_{pp}$ (I take the maximum value V_pp for real or imaginary received signal) and I try to maintain this value into a range such that $V_{pp,lim1} < V_{pp} < V_{pp,lim2}$ , to avoid quantization noise and saturation. To do so, I simply increase or decrease the gain in function of the value $V_{pp}$, starting with IF and BB gains to finish with the amplifying RF gain.

The problem is the following : My algorithm is not optimised. Given that the possible changing steps are not linear and different for every type of gain (14dB for RF, 8dB for IF and 2dB for BB), I sometimes fall into situations where my gains are continuously changing between two situations and do not manage to stabilise even if I keep the distance constant between emitter and receiver.

The question is : how can I manage to have the more constant output amplitude so that my dynamic range is optimised for the ADC?

Constraints :

Gain steps restricted : 14dB for RF (0 or 14dB), 8dB for IF (from 0 to 40dB) and 2dB for BB (from 0 to 62dB)
Order of changing gains : "Increase or decrease the IF and baseband gain controls roughly equally to find the best settings for your situation. Turn on the RF amp if you need help picking up weak signals" (source : https://github.com/mossmann/hackrf/wiki/FAQ )

To avoid quantization or limit cycle oscillations, you can add hysteresis to your AGC decision loop (a form of "debouncing").


Number of samples in one period of subcarrier after IFFT I have a quick question about the number of samples in one period of a sub carrier in OFDM systems. The IFFT process generates N-parallel sub carriers in discrete time domain which are sent simultaneously over the sub carrier frequencies.

What is the number of samples in one period of a sub carrier after IFFT?

Thanks for your reply. Let me give you more details about my query with an illustration as below. The picture shows the sub carriers as analog but lets take them as discrete before DAC and CP addition. The system has 4 sub carriers and hence they have different frequencies. My questions is here at this point; again assuming they are just the output of the IFFT and digital-
Does each sub carrier has 4 samples in every cycle? So lets say the duration of OFDM symbol is T seconds. The subcarrier1 has 1 cycle and I am visualizing that the sub carrier has 4 samples during the T seconds.
Subcarrier2 has 2 cycles (period) and hence it has 2x4=8 samples during the OFDM T seonds.
And for subcarrier3 it has 3x4=12 samples.
For subcarrier4, there are 4x4=16 samples during the OFDM symbol (T seconds length).

And in total, the OFDM symbol has 36 samples in the symbol length T seconds.

Am I correct interpreting this?

If you'd like to think analog, an OFDM signal can be written as a sum of weighted complex sinusoidals:
$$
x(t)=\sum_{k=0}^{N-1}X_k \exp{\left( \mathrm j 2\pi\frac{kf_\mathrm{s}}{N}t \tag{1}\right)},
$$
where $N$ is the number of subcarriers, $f_\mathrm{s}$ is the sampling frequency and $X_k$ denotes the subcarrier values. For a digital implementation, (1) is sampled at $t=n/f_\mathrm{s}$:
$$
x(n/f_\mathrm{s})=\sum_{k=0}^{N-1}X_k \exp{\left( \mathrm j 2\pi\frac{kn}{N} \tag{2}\right)},
$$
which corresponds to the inverse discrete Fourier transform (IDFT), implemented by an IFFT. The IFFT has $N$ input values and $N$ output values. That's ok, because the IDFT is N-periodic and thus any additional sample would be redundant. Consequently, every subcarrier in (1) is represented by $N$ samples, regardless of its number of periods during the OFDM symbol length.

If you try to imagine the digital OFDM time domain signal as a superposition of $N$ sampled sinusoidals, the figure you posted has two problems, in my opinion:

The first subcarrier in the image obviously has frequency $f_\mathrm{s}/N$ and the fourth subcarrier has frequency $f_\mathrm{s}$. That's wrong. From (1) you can see that the first subcarrier should have frequency 0 and the fourth should have frequency $(N-1)f_\mathrm{s}/N$. You can reorder the subcarriers (e.g. from negative to positive) but the DC subcarrier must be included and there must not exist a subcarrier with frequency $f_\mathrm{s}$. (One could argue that subcarriers with frequency 0 and $f_\mathrm{s}$ are actually the same but I think it's a poor choice for an illustration)
The image only shows the real (or imaginary) part of the sinusoidals. That could give you the impression that some aliasing is going on when sampling. I will explain in a second what I mean.

Please consider the following figure that I have created instead

It contains the subcarriers from $k=0$ to $k=3$ (lines) and their samples (dots). It also shows the real part (blue) and imaginary part (red). Now consider the red line in the first and the third subfigure: their samples are exactly the same. However, the blue (imaginary) part is different. That is, $N=4$ samples are sufficient here to represent all subcarriers.

Here is the Matlab code I used to create the above figure:

N = 4; % num subcarriers
os = 128; % oversampling
t = linspace(0, 1-1/(N*os), N*os); % time
n = 0:N-1; % discrete time

figure;
cnt = 1; % helper
for k = 0:N-1;
    subcarrier = exp(1i * 2*pi*k*t);

    subplot(N,1,cnt);
    plot( t, real(subcarrier), 'b-');
    hold on;
    plot( t(n*os+1), real(subcarrier(n*os+1)), 'b.', 'MarkerSize', 20);

    plot( t, imag(subcarrier), 'r--');
    hold on;
    plot( t(n*os+1), imag(subcarrier(n*os+1)), 'r.', 'MarkerSize', 20);

    if(k ~= N-1)
        set(gca, 'XTick', []);
    else
        xlabel('t/T');
    end
    ylim([-1 1]);
    title( sprintf('k=%d', k));

    cnt = cnt + 1;
end


Factored form vs partial fraction form? I have already understood partial fraction and here is link for my relevant DSP SE question

Finding inverse z transform for two sided ROC?

But now i want to know, is there any difference between partial fraction form and factored form in signal processing context?

For example I have a z transform $$Y(z)=\frac{(z^2−z)}{(z^2+1.3z+0.3)}$$

a)What will be its partial fraction form?

b)What will be its factored form?

The partial fraction form helps in calculating the z-transform inverse since we can get the inverse of each term in partial fraction by inspection and hence get the inverse of the whole transfer function. The factored form can directly give poles and zeros by equating the numerator and denominator to zero.


Head Phone Transfer Function MATLAB Is it possible to model the acoustic transfer function in MATLAB for a common set of headphones? i.e. I want to model the summation of a noise sample and the audio from the output of the headphone e.g.

signal = exteriorNoise + Music

Where signal would be the information the user would hear? I would imagine headphones would implement some attenuation and low pass filtering, but I cannot find any documentation. I am aware every headphone / in-ear would be different, but I merely wish to make a more realistic model that basic signal addition.

Many Thanks.

It's actually a band-pass characteristic, since very low frequencies will be inaudible on headphones. Firstly, searching for "headphone frequency response" brings up a lot of results, especially on an image search. Then, the site headphone.com seems to offer frequency characteristic plots for many of the headphones they sell.

Unfortunately, you often only get graphs. However, you can use tools like g3data (plot digitizing tools) to assist you in recovering data from pictures or even scans of plots. I used such tools a lot to trace plots from scientific papers for which data is usually not available.


Is there a strategy for discrete control of a system with dynamics near sample rate? I'm trying to control a system where the controller sample rate is physically fixed and the plant has significant dynamics on the same order as the sample rate. I understand that one would prefer to have the sample rate considerably faster than the plant dynamics, but the physics of this system are such that this is inherently impossible. I can get some control of the system by hand-tuning a PID, but it seems considerably sub-optimal.

Is there a strategy for developing a controller like this?

Edit: I should add that the plant is this situation is a relatively complex LTI plant with reverberations caused by pure time delays. The delay cycles are a bit slower than the Nyquist frequency. The plant is stable.

Depending on your plant you might be able to do something. If the plant is unstable, but dominantly linear, you should be able to get a very accurate model using exact discretization (zero-order hold). Either you will be able to stabilize your plant using PID, or PID will be of too low order, in which case you need to switch to a control law matching you plant order. The H-infinity synthesis framework (robust control) will do that for you. Here you can also include uncertainty weights to limit control law gain, which might help esp. if your system is not linear. Depending on the uncertainty you might not be able to find a feasible control law for stabilization, however...

If your plant is stable, or stabilizable, you can further improve performance using feed-forward control. Here the sampling time should not be of much concern, and you can use something like zero phase error tracking.


Segmentation of microscopic images This is a classical image segmentation problem: I need to make a visual measurement of particles.
I've been using watershed, kmeans and mean shift but it doesn't work well on this image. Could you please recommend me some other methods/open source code?

P/S: Any manual/semi automatic segmentation is also ok with me. I just need the correct segmentation in a reasonable not time consuming way.

Since it is very hard to guess what result you need, proposing any automatic segmentation method is very hard.

Let me instead propose two different semi-automatic segmentation with markers.

Segmentation using Minimum Spanning Forests with markers

This method is based on a version of watershed [Cousty et. al.: Watershed Cuts: Minimum Spanning Forests and the Drop of Water Principle, 2009].

Basically, you calculate a version of watershed (there are many versions), called watershed cuts, which gives you a (pretty severe) over-segmentation.

After that, you mark your one object (or all of your objects, 1 marker per object + 1 for the background if there is any) you want to find, and it calculates the most likely segmentation based on pre-calculated watershed.

I don't want to go into details, since there is pseudo-code and example images provided (one is even a cell image) in the article, but if you want me please leave a comment and I can add some more information.
Segmentation using $\alpha$-trees (or $\alpha$-connected components) with markers

Here, the "starting over-segmentation" is actually simply the segmentation to flat-zones (equal-intensity areas in the image). I am going to explain a little bit of theory here.

An $\alpha$-connected component of a pixel $p$, or $\alpha$-$CC(p)$ is obtained as follows:

first you choose a connectivity (4- or 8-conn.; 4-conn simply means the neighbors of a pixel are up, down, left, and right)
if the intensity difference between two neighboring pixels is $\leq \alpha$, those pixels are connected
$\alpha$-$CC(p)$ contains all the pixels connected to pixel $p$ as well as all the pixels that can be reached from pixel $p$ through a connected path

Now, for various reasons (I can expand the answer again if somebody asks), if you follow the $\alpha$-$CC$ around one pixel, its area usually grows slowly as you increase the $\alpha$, then at some point you have a sudden and rapid increase in the area, and then it continues to grow slowly. Usually, if you look at the $\alpha$-$CC$ just in the moment before the rapid size increase, it will give you your object.

Again, you can use markers to mark which $\alpha$-$CC$'s you want to "follow". The advantage here is that you can calculate the $\alpha$-hierarchy once and then do all the calculations on the same hierarchy.

I know it sounds a bit confusing with all the new terms unless you know something about $\alpha$-trees, but here's some references to help:

Najman, Soille: On Morphological Hierarchical Representations for Image Processing and Spatial Data Clustering, 2012 - a good reference for $\alpha$-hierarchies in general
Serna, Marcotegui: Attribute Controlled Reconstruction and Adaptive Mathematical Morphology, 2013 - these guys have done something very similar to what I am proposing. One of the proposed and explored usages in their article is segmentation (they're just using more attributes, not just area).
Merciol, Lefevre: Fast Image and Video Segmentation Based on Alpha-tree Multiscale Representation, 2012 - this is more or less exactly the usage I described.


What is the name of this very simple spectral subtraction technique? Let $S = X + N$ be the sum of two audio signals $X$ and $N$ which are both stationnary (let's think X is a constant volume 440 Hz sinusoid and N is constant volume noise).

If the sum S has a -20 dB volume and N has a volume of -30 dB, what is the volume of X?   (could be RMS or peak volume, it doesn't matter here).

The answer is

20*log10(10^(-20/20)-10^(-30/20)) ~ -23.3

i.e. X has a peak volume of -23.3 dB.

(of course this is not true if X and N have phase-cancellation, but except this case, it works).

Question: what is the name of this computation? i.e. find the volume of X given the volume of N and X + N?

I'm using an application of this for STFT denoising (with a noise template which is N):

if a FFT bin has -20 dB amplitude (signal S), and the noise has amplitude -30 dB amplitude for the same bin (signal N), what would be the amplitude of the denoised signal X ? Answer: -23.3 dB, thus this FFT bin should be lowered by 3.3 dB.

It works quite well for my noise reduction application (again here X is nearly a constant sinusoid and N constant noise), but I haven't found a name for this simple technique. What would be the name?

Estimating the power of a sum of signals depends on coherence or incoherence assumptions. Details are given for instance on incoherent signal summing or coherent signal summing. Those can provide an estimate of the noise level. They can be called power sums or voltage sums. Other details in Voltage and power add differently. I'll be digging for other names.

The expectation of the energy of a sum $v$ of variables $v_1$ and $v_2$ is:

$$ E((v_1+v_2)^2) = E(v_1^2)+E(v_2^2)+2E(v_1 v_2)\,.$$

Three specific cases. If $v_1=v_2$, you get $+10\log 4 = +6 $ dB. If sources are uncorrelated, $E(v_1 v_2)=0$, so you get  $+10\log 2 = +3 $ dB. If sources are oppositive, they cancel each other, and the total energy becomes $-\infty$ dB.

Then, subtracting the noise level to the observed signal is an instance of (scalar) hard thresholding, called spectral subtraction in the Fourier domain. As eluded to in your Spectrogram with square or non-square magnitude of STFT: power vs. magnitude question, scalar spectral subtraction on the observed STFT ${S}(\tau,\omega)$ can be performed with a power law, as:

$$|\overline{X}(\tau,\omega)|^\alpha=\max\left(|{S}(\tau,\omega)|^\alpha - \lambda|\overline{N}(\tau,\omega)|^\alpha,0\right)\,.$$

The recovery of the noiseless phase is a full different world.


Spectrum of windowed version of original continuous signal Suppose we have the complex signal $x(t)= \exp(j\omega_0 t)$. Using the properties of Fourier transform we can prove its CTFT is Dirac $\delta$ function.

If any one ask me about the spectrum of $x(t)$ "Does $x(t)$ has continuous spectrum or discrete spectrum", my answer will be "The spectrum of $x(t)$ is discrete".

Now, if I apply a rectangular window to the complex exponential $x(t)$ in the time domain and then take CTFT I will end up with $\mathrm{sinc}$ function. Now the spectrum of the windowed complex exponential is continuous and not discrete.Is this interpretation true?

If you don't understand the difference between the Continuous Time Fourier Transform (CTFT), the Discrete Time Fourier Transform (DTFT) and the Discrete Fourier Transform (DFT), now would be a good time to read about them. The very short version is that the DTFT yields the (continuous-valued) spectrum of a sequence (i.e., a sampled signal). The DFT computation results in a sampled version of the DTFT. To apply the DFT requires a finite number of samples (i.e., a time-domain window) whereas such restrictions are not placed on the DTFT in general. On the other hand, the CTFT deals with continuous time signals. There is a lot more to all of this, and I recommend you read more.

It is true that a complex-valued signal $x(t) = \exp\left(j \omega_0 t\right)$ and the Dirac delta function form a CTFT pair. I can agree with you that the spectrum of this signal is discrete (nonzero for a finite number of frequencies, in this case a single frequency).

After applying a rectangular window in time to the complex-valued signal $x(t)$, the CTFT is a frequency-translated sinc function centered at $\omega_0$ with a lobe width inversely proportional to the size of the time domain window. This shows us that the result of truncating a time domain signal with infinite support in time and a discrete spectrum in frequency can lead to a new time domain signal with finite support and a continuous spectrum. So, the answer is yes, the interpretation given in the question is true.


Two basic questions related to complex Fourier series Doing excersies from Richard Hammings book about digital filters I've got two questions about them:

1) Fourier expansion of $g(x) = \sin^5(x)$. Provided answer is: $\sin^5(x)=5\sin(x)-20 \sin(3 x)+16\sin(5 x)$. Which is obviously contradicts well-known solution. Is it mistake in a book or author meant something another?

2)The second excersice which confusing me is expansion in compex notation of $g(x) = \exp^ {\alpha x}$. I got the following form of coefficients: $c_{k} = {\frac {1} {(2\pi)}} (\exp^{\alpha\pi}  - \exp^{-\alpha\pi}){\frac {(-1)^k} {(\alpha-ik)}}$. These coefficients looks a little bit frightening to me and I don't know how to check them... Any ideas how to validate solution?

Thanks in advance!

First, or second, with complex notation and the binomial theorem,

$$\sin^5 x=\left(\frac{e^{ix}-e^{-ix}}{2i}\right)^5 = \frac{e^{i5x}-5e^{i3x}+10e^{ix}-10e^{-ix}+5e^{-3ix}-e^{-i5x}}{32i} = \frac{\sin(5x)-5\sin(3x)+10\sin x}{16}$$

On the other hand,

\begin{align}
\sin(5x) &= Im(e^{i5x}) = Im([\cos x+i\sin x]^5) \\
&= 5\cos^4x\sin x - 10\cos^2x\sin^3x + \sin^5x\\
&= 5(1-2\sin^2 x+\sin^4 x)\sin x-10(1-\sin^2x)\sin^3x+ \sin^5x\\
&= 5\sin x - 20\sin^3x +16\sin ^5x
\end{align}

so you got the coefficients right, but the notation reversed.

In the second task, is the star connected to convolution, cepstrum transforms or similar?

No, it is just $g(x)=\exp(α x)=e^{α x}$ defined on the interval $[-\pi, \pi]$. Then

\begin{align}
c_k
&=\frac1{2\pi}\int_{-\pi}^\pi g(x)e^{-ikx}dx
 =\frac1{2\pi}\int_{-\pi}^\pi e^{(α-ik)x}dx
\\
&=\frac1{2\pi}\left[\frac{e^{(α-ik)x}}{α-ik}\right]_{x=-\pi}^{x=\pi}\\
\end{align}


Can I get more information from this signal other than magnitude? First a bit of background. I'm working on a metal detector that produces a pulse across a coil. The pulse signal is influenced if metal is present within the coil aperture.This link provides more information on Pulse Induction Metal Detectors. The ideal signal looks like this:

Detection of metal relies on back EMF and the signal recovers after a pulse. I'm working on a design that removes all aspects of the signal except for this "recovery curve" aspect. In MATLAB there is a clear definition between the recovery curve when there is metal and no metal:

At the moment I am calculating the average of the entire "recovery curve" over a number of curves (reading in 50 curves and then calculating the average of each point in the curve) to create an average representative curve and then turning the curve in to just a single number i.e. the average single value of 50 curves. Here is a number of recovery curves of 4 different metal pieces.

And here is when I calculate the average of 50 for each metal piece to create one representative curve:

I am then turning the curve into a number and if it exceeds a set threshold then metal is considered present. Effectively I am turning the average of 50 curves into a single magnitude and then comparing it with a no metal reference value.
While this works, I feel like a better comparison can be made. I want to be able to generate a profile of what a no metal present "recovery curve" looks like and then compare subsequent curves to my reference curve. At the moment I am just using single magnitudes and a comparator to my reference curve. Are there any other techniques I can use to make a more detailed comparison?
Because the shape of the curve changes, i.e. the intial slope can be steeper, the knee of the signal can become rounder etc. is there any kind of technique or filter that I can use to more accurately profile the shape of the curve and compare the shapes rather just the average magnitude? I'm a complete newbie in regards to DSP so If I can be given some direction it would be appreciated.
thanks.

There's a lot of answers -- look at what's different between the metal and no-metal case, and try to make some signal processing that'll do the best job of differentiating between them.
Something that I suspect that will work well would be to start a clock as soon as the pulse goes negative (possibly according to the transmitter timing, if that will be reliable).  Capture the height of the peak, then integrate the area of the whole response over time.  That peak to area ratio should be very indicative of the presence of metal.
My gut feel is that the ratio of the peak to the integral of the overall response should be more indicative than the actual value of the integral -- but if the pulses are repetitive enough, just that integral may be enough.  If you want to go down that road, be sure to try to detect every different thing you may be concerned with, at every different range.
It may be more reliable to capture the peak by integrating for a short time, and comparing that number to the integral of a longer time span.
Continuing in this direction, if you were to integrate the pulse over time and look at the shape of the resulting integrated pulse it should be very indicative of what's going on.  Once again my gut feel is that you want to look at ratios -- I'd normalize the integrated pulse to its final value.  The points at which it reaches 1/4, 1/2, 3/4, etc., of its final value should be very telling.


How to find fundamental frequency of Fourier Series The original problem is from the Problem Set 7 of MIT OpenCourseware: Find the Fourier series coefficients for
$$
x(t)=sin(10\pi t+\frac{\pi}{6})
$$
What I did is to rewrite it in exponential form $\frac{1}{2j}e^{j\frac{\pi}{6}}e^{j10\pi t}-\frac{1}{2j}e^{-j\frac{\pi}{6}}e^{-j10\pi t}$, and take $\omega_0=10\pi$ as the fundamental frequency. The non-zero coefficients I got are $a_1=\frac{1}{2j}e^{j\frac{\pi}{6}}$ and $a_{-1}=-\frac{1}{2j}e^{-j\frac{\pi}{6}}$.

However, the solution takes $\omega_0=2\pi$, which gives $a_5=\frac{1}{2j}e^{j\frac{\pi}{6}}$ and $a_{-5}=-\frac{1}{2j}e^{-j\frac{\pi}{6}}$.

I can't find a proper explanation for this problem. What should I take as fundamental frequency? Can the fundamental frequency be arbitary number for continuous time Fourier series?

Here is the link of the problem set: Problem set 7

The link of the solution: Solution

Your solution is correct. Obviously, any periodic function with period $T$ is also periodic with period $kT$, $k\in\mathbb{Z}^+$. In the given example, the smallest possible period is $T=1/5$, which is equivalent to $\omega_0=10\pi$. There is no obvious reason why one should choose $T^\prime=kT$ and $\omega_0^{\prime}=\omega_0/k$, with $k>1$. Nevertheless, such a choice, even though unmotivated, is equally correct, because the resulting Fourier series are identical:

$$\begin{align}x(t)&=\sum_nc_ne^{jn\omega_0t}\\&=\sum_nc_n^{\prime}e^{jn\omega_0^{\prime}t}\\&=\sum_nc_n^{\prime}e^{j\frac{n}{k}\omega_0t}\\&=\sum_mc_{km}^{\prime}e^{jm\omega_0t}\end{align}\tag{1}$$

With $\omega_0^{\prime}=\omega_0/k$ we get

$$c^{\prime}_n=\begin{cases}c_{n/k},&n=mk\\0,&\textrm{otherwise}\end{cases}$$


After advice about detecting focus quality of objects in a photo detected using YoloV3 I've spent the last couple of days playing with YoloV3, and have had very good results.  My use case is sports photography, and the object detection for people/bikes etc is very very good, I'm very impressed.  In the future I think I'll train it with my own dataset to improve it further, but out of the box it does a fantastic job already.

What I want to improve:

Once an object has been detected, how can I generate some sort of metric to quantify how well-focused it is?

Past/current approaches

1) My first thought/Google, was "variance".  First I turn the image grayscale, then use a Laplacian convolution to highlight the edges.  Then simply compute the variance of the pixels in the bounding box.  High numbers "probably" mean high contrast i.e. high focus, whereas lower numbers would tend to mean low contrast and probably poor focus.

It works pretty well, but not 100% reliable.  Imagine you have a person stood up with their arms spread wide, against a skilled bokeh background.  Most of the bounding box is blurred background, so the variance ends up being low.

2) To improve upon this, I came to the conclusion that there will almost always there be a significant portion of the object somewhere around the center of the bounding box.  I set my criteria to center square, 20% of the width of the BB and 20% of the height.

Unfortunately this throws up scenarios where that square happens to be 10% background, between someone's arm and body when running around a tight corner facing the camera, etc.

3) "OK, a cross".  Thickness equal to 20% of the BB, up and down through the middle and left and right through the middle.

Not bad, not bad.  Still getting a lot of background on some images though, as the edges of the box are where background is going to tend to reside.

4) "OK, a reduced cross".  Same as above, but only extending from the center 2/3rds of the way out to the edges.

Almost fantastic.  With the caveat that sometimes you end up with just a competitor's chest, and if they're wearing a single-colour top.... the variance isn't all that.

Examples:

In this photo, the motorbike (close enough...) apparently has great focus, while the person isn't so good.  Mainly due to the near-uniformity of his central cross.

Here's a more troubling example.  Look at that variance, 5 FFS!

So I think that's the end of that approach.

The future...

I could go on and on with this, and I'm ALWAYS going to end up with some photos that it just doesn't work well for.

I think a different approach is needed.

One thought is just to take the largest variance over a small region, say a 10% width/height square that roams across the bounding box.

But then if the foreground is completely out of focus, and the background is sharp, we'll get a false positive from the background.

Anyone cleverer/more experienced than me have a fantastic solution for this?

It's clearly possible, not least because http://remove.bg and PhotoShop already do a fanastic job of separating foreground from background.  But how?

EDIT: I completely neglected to mention that I'm using a Laplacian convolution on a grayscale version of the photos before computing the variance, to detect the edges.

Here is what I would try on the source image:

Split your image into 5x5 pixel blocks (maybe 3 maybe 7, who knows?)

Create output image one fifth (third, seventh, ??) size

For each block
  For each color channel

    Find best fit plane
    Measure RMS of (pixel value-plane value)

  Next

  Set output pixel to  RMS(R,G,B)

Next

In blurry/plain areas the RMS should be low.  With sharp lines it should be high.

Ced

So, I thought I'd give this a go with 5x5 block.

Here are some results.

Clearly (pun intended), not only is this algorithm an edge detector, but it also a focus metric at those locations.

Those are kind of dim, so I multiplied by 5.

If nothing else, they are kind of cool looking, but I think you can tell where the well focused areas are.

The 5x5 takes a while to process, so I decided to try a 3x3 fit on a simple 4 point gradient estimate (Down from a 16 point gradient estimate).

Still (pun intended again), the blurry lady pic proves this is a focus metric.

Here is the relevant code called for each color channel:

'==================================================================
Public Sub FindFocus(ArgV As Float[]) As Float[]

        Dim w, h, x, y As Integer

        w = ArgV.Bounds[0]
        h = ArgV.Bounds[1]

        Dim theFocus As New Float[w, h]

        For x = 1 To w - 2
          For y = 1 To h - 2
            GoSub CalculateFocusAtPoint
          Next
        Next

        Return theFocus

'-------------------------------------------------------------------
CalculateFocusAtPoint:

        Dim dx, dy As Integer
        Dim a, b, c As Float

'---- Plane Estimate:  z = ax + by + c

        a = (ArgV[x + 1, y] - ArgV[x - 1, y]) * 0.5
        b = (ArgV[x, y + 1] - ArgV[x, y - 1]) * 0.5
        c = ArgV[x, y]

'---- Calculate the RMS of the NonPlanar

        Dim v, e, s2 As Float

        s2 = 0

        For dx = -1 To 1
          For dy = -1 To 1
            v = a * dx + b * dy + c
            e = ArgV[x, y] - v
            s2 += e * e
          Next
        Next

'---- Set the Value

        theFocus[x, y] = Sqr(s2 / 9)

        Return
End
'==================================================================

Finally, just because I could, I turned your pic into a "drawing".


how can we tell that OTSU binarization is feasible or not. OTSU is a very popular image binarization method, and Matlab provides a function to implement this method, see graythresh for details. In this function, it has two outputs, one is the threshold calculated with OTSU algorithm, the other is called effectiveness metric EM, and it claims:

    The effectiveness metric is a value in the range [0 1] that indicates the
effectiveness of the thresholding of the input image. The lower bound is
attainable only by images having a single gray level, and the upper bound is
attainable only by two-valued images.

I think it would be nice if we can use effectiveness metric EM to tell whether binarization is fine or not. If EM value is slow, then we will use different binarization method. It would be very useful in the context of document image binarization. Then, a natural question is how can we set this threshold reasonablly. Any ideas on how EM is calculated?

Assuming you have MATLAB, you can just look at the source code, where you will see that the EM is defined as

  em = maxval/(sum(p.*((1:num_bins).^2)') - mu_t^2);

where

I = im2uint8(I(:)); # the input image
num_bins = 256;
counts = imhist(I,num_bins);
p = counts / sum(counts);
omega = cumsum(p)
mu = cumsum(p .* (1:num_bins)');
mu_t = mu(end);
sigma_b_squared = (mu_t * omega - mu).^2 ./ (omega .* (1 - omega));
maxval = max(sigma_b_squared);


Why Cramér spectral representation and not DTFT for stochastic process In a lot of time-series analysis references I find (written by mathematicians or statisticians rather than engineers), I find the following signal decomposition for a stochastic process, termed the "Cramér representation" (e.g. eqn 8.11 of this reference):
$$
X[n]  = \int_{\langle 2\pi \rangle} e^{-j\omega n} d Z(\omega)
$$
The factor $dZ(\omega)$ is referred to as a spectral increment. I found another reference (ref, eqn 77) that said that the spectral increments are orthogonal (w.r.t. the expectation operator) if the process is stationary.
Compare this to the inverse discrete-time Fourier Transform (IDTFT), non-normalized, angular frequency convention (eqn 4.2.28 of Proakis & Manolakis, Digital Signal Processing, 4th ed):
$$
X[n] = \frac{1}{2\pi} \int_{\langle 2\pi \rangle} e^{j\omega n} X(\omega) d\omega
$$
With the exception of trivial differences in convention (minus sign in the exponent, normalization factor), the two representations appear to be the same. Ignoring the minus sign convention on $\omega$ for now, I am tempted to just conclude:
$$
dZ(\omega)  = \frac{1}{2\pi} X(\omega) d\omega
$$
but I suspect there is a deeper mathematical reason why this would be wrong and that the statistics literature uses spectral increments instead.
Why do statisticians prefer the Cramér representation? Are there any computational or theoretical advantages to using it?
Does it have something to do with the convergence (or existence) of some type of integral? Or some issue relating to the fact that $X[n]$ is explicitly a stochastic process in the Cramér representation whereas the DTFT might rely on the signal being deterministic.
I wonder this because engineering education (at least mine was this way) tends to abuse notation or gloss over certain mathematical difficulties because those nuances wouldn't matter for the situations in which an engineer would be using said mathematical tools.
For instance, as an undergrad I never had to learn what a Lebesgue integral was, even though I was implicitly computing Lebesgue integrals in my probability course.

I will introduce some terminology and intuition that will be helpful when reading other references. It will be neither complete nor completely rigorous.

The measures that we first encounter in real analysis assign sizes (non-negative real numbers) to measurable subsets of $\mathbb{R}$; Lebesgue measure is the measure that agrees with the intuition we build in calculus (the measure of the interval $[a,b]$ is $b-a$, etc).

$Z$ is a measure, but it is a stochastic measure†. It does not assign numbers to measurable subsets of $[0,2\pi]$. Rather, it assigns a random variable to each such subset:
\begin{equation}
X_A = \int_{A}dZ(\omega).
\end{equation}
The convergence of the integral on the right-hand side is an issue that I would rather not try to explain (to you or to myself).

In particular, the $Z$ used for WSS processes is an orthogonal stochastic measure. One result is that random variables assigned to non-overlapping sets are independent of one another.

If $A$ is a Lebesgue-measurable set, then $Z(A)$ is a random variable, and the expectation of $\left|Z(A)\right|^2$ is
\begin{equation}
\mathsf{E}[\left|Z(A)\right|^2] = \textrm{Lebesgue measure of $A$}.
\end{equation}
Hence, Lebesgue measure is "under the hood" even if we stick with the notation $dZ(\omega)$.

Just as we can use Lebesgue measure to integrate functions over subsets of $\mathbb{R}$, we can use $Z$ to integrate functions over subsets of $[0,2\pi]$ (such as all of $[0,2\pi]$).

Let $\mu$ be Lebesgue measure on $\mathbb{R}$, and let $\nu$ be another measure on $\mathbb{R}$. $\nu$ is said to be absolutely continuous with respect to Lebesgue measure if there is a function $f$ such that $d\nu = fd\mu$, or the measure $\nu(A)$ of $A$ is equal to
\begin{equation}
\nu(A) = \int_{A}f(x)d\mu(x).
\end{equation}
The function $f$ is called the Radon–Nikodym derivative of $\nu$ with respect to $\mu$.

Not all measures are absolutely continuous with respect to Lebesgue measure. The example most familiar to electrical engineers is Dirac measure. Lebesgue measure assigns measure zero to any set consisting of a single point, and a measure that is absolutely continuous with respect to Lebesgue measure must do the same. But the Dirac measure $\delta_0$ assigns measure 1 to the set $\{0\}$ and to any set that contains $0$. Since $\delta_0$ is not absolutely continuous with respect to Lebesgue measure, $d\delta_0$ cannot be written as $fd\mu$.

There are also more exotic measures that are not absolutely continuous with respect to the Lebesgue measure.

I have found no evidence of the notion of absolute continuity of stochastic measures.

EDIT: While theoretical results about spectral representations of WSS processes are crucial for applications, the $dZ$ notation may be off-putting and perhaps even doubt-inducing. I suspect that writing $Y(\omega)d\omega$ for $dZ(\omega)$ is a useful abuse of notation that allows the user to manipulate symbols as though some analogue of the Radon-Nikodym derivative existed. Rigor can be added after the fact.

Note that rigor might arrive decades after the fact. Plenty of ideas seem to work just fine without complete mathematical rigor.


How to remove weak edges in a binary photo? See attached photos for reference.

I tried watershedding, it breaks apart the main body and doesnt even touch the part I want disconnected!

morphology is no good, the holes in the main body easily break in two before the edges I want disconnected come off.

imfill doesn't solve this, it fills the hole between the two trouble areas as well as the holes in the main body.

I have thousands of similar photos, so I need to come up with something robust.

I'd appreciate any help I could get!

When you are dealing with a binary image (black and white with 1 bit per pixel intensity) you loose all the continuous gradient information, and therefore cannot make use of edge detection (or removal) algortihms that rely on it.

Instead you can return to geometric features such as line length, edge thickness (indeed line thikness), surface area etc to classify features found in the image.

A weak edge should have a corresponding definition in terms of feature geometries so that it can be detected and removed (replaced) if necessary.


Does down-mixing a signal before taking the fourier transform aid in signal detection? Let us assume that we have a narrowband signal $s[n]$, at frequency $f_1$, and let us further assume that this $f_1$ is just below $\frac{f_s}{2}$, where $f_s$ is the sampling rate.

Furthermore, let us also assume that we have a fixed number of time-domain samples of this signal, $N$, and that our FFT length is also $N$.

The objective is to detect this signal $s[n]$'s presence in the frequency domain. (Further assume that we know to expect this signal in or around that frequency $f_1$).

My question is simple: With all things being equal, which is better from a detection perspective:

Simply do the FFT of this signal and look at the peak at that frequency bin

Or

BPF a band of where we expect the signals' frequency $f_1$ to be, (eg, BPF this band $f_1 - \epsilon \text{ to} f_1 + \epsilon$). Then, mix this band down such that $f_1 + \epsilon$ is guaranteed to be just above DC, and LPF that entire band. Now, take the FFT, and detect energy.

Would the second method be better from a detection point of view because we have more samples per cycle, whereas the first one has just over 2 samples per cycle?

Thanks,

TLDR: Is it better from a detection standpoint to FFT a signal or down-mix it and then FFT it (with appropriate filtering)?

If you are doing are doing this in both the analog and digital domain, then there may be significant differences due to trade-offs between the noise sources in the analog stage(s) in the filters and the mixer oscillator frequency, and the corresponding sampling jitter and quantization noise in converting to digital.

Even strictly in the digital domain, the bandpass and lowpass filter responses may have different effects on how potential stop-band noise might affect your detection accuracy, compared to how windowing (rectangular or otherwise) affects the continuous frequency response of each FFT bin.


Determine modulation function from spectrum? I have a two periodic signals $x(t)$ and $y(t)$ which are related to each other by a periodic modulation function $\alpha(t)$:
$$y(t) = \alpha(t) x(t)$$
If $x(t)$ is a sinusoid, is it possible to determine $\alpha(t)$ from the spectrum of $y(t)$? How might I do this?

The Fourier transform is reversible, so the same amount of information is contained in the spectrum (not in the power spectrum, that doesn't contain the phase, and hence, only half of the information)! Just convert your spectrum back to time domain.

However, if you just got the power spectrum plot of a signal, you'll find that many modulations have "characteristic" spectral shapes. With a bit of practice, you can tell binary GMSK from an 4-FSK pretty easily.

Now, in general, the problem of detecting a signal and its modulation is pretty hard – governments pay loads of money for well-working signal classificators, because in many cases, you're looking at a very weak $y(t)$ within a mixture of hundreds, and it's still hard to guess whether you're looking at an 16-QAM or an 8-PSK if you don't have timing information.

Thus: without a lot of additional info on what you're looking at, guessing the modulation class of $\alpha$ from $|Y(f)|$ is very hard. And it's harder to actually estimate $\alpha(t)$ from that.


Relationship between DCT and PCA I have a basic implementation knowledge of the 2D 8x8 DCT used in image & video compression. Whilst reading about Principle Component Analysis, I can see a lot of similarity, albeit PCA is clearly more generic. When I've read about DCT previously it was always presented in relation to DFT. So my question is how can the DCT be derived from a PCA perspective? (even a hand-wavey explanation is sufficient)

Many thanks

The main difference between DCT and PCA (more precisely, representing a dataset in the basis formed by the eigenvectors of its correlation matrix - also known as the Karhunen Loeve Transform) is that the PCA must be defined with respect to a given dataset (from which the correlation matrix is estimated), while the DCT is "absolute" and is only defined by the input size. This makes the PCA an "adaptive" transform, while the DCT is data-independent.

One might wonder why the PCA is not used more often in image or audio compression, because of its adaptivity. There are two reasons :

Imagine an encoder computing a PCA of a dataset and encoding the coefficients. To reconstruct the dataset, the decoder will need not only the coefficients themselves, but also the transform matrix (it depends on the data, which it does not have access to!). The DCT or any other data-independent transform might be less efficient in removing statistical dependencies in the input data, but the transform matrix is known in advance by both the coder and decoder without the need for transmitting it. A "good enough" transform which requires little side information is sometimes better than an optimal transform which requires an extra load of side information...
Take a large collection of $N$ 8x8 tiles extracted from photos. Form a $N \times 64$ matrix with the luminosity of these tiles. Compute a PCA on this data, and plot the principal components that will be estimated. This is a very enlightening experiment! There is a very good chance that most of the higher-ranked eigenvectors will actually look like the kind of modulated sine-wave patterns of the DCT basis. This means that for a sufficiently large and generic set of image tiles, the DCT is a very good approximation of the eigenbasis. The same thing has also been verified for audio, where the eigenbasis for log- signal energy in mel-spaced frequency bands, estimated on a large volume of audio recordings, is close to the DCT basis (hence the use of DCT as a decorrelation transform when computing MFCC).


Identify whether to have unique output in this ARMA system my friends and I have no idea whether to diagnose to have the specific output from  the formula below:

The question is that Assume we are given a specific input signal x[k]x[k]. Explain and justify whether or not there is a unique output y[k] satisfying this difference equation.
Would somebody please answer this question asap??
Thank you

As I mentioned in my comment, I think there is a typographical error in the equation. I think it is supposed to be
\begin{equation}
y[k] = \sum_{i=1}^{M}a_iy[k-i] + \sum_{j=1}^{N}b_jx[k-j].
\end{equation}
This is an example of a linear constant-coefficient difference equation. These are introduced very early in most textbooks about digital signal processing. For example, they are introduced in Section 2.5 (page 33) in the 1989 edition of Discrete-Time Signal Processing by Oppenheim and Schafer.
Given an input signal $x$, there will not be a unique solution signal $y$. This is because a specific solution requires specifying the input signal as well as initial conditions for $y$. Different initial conditions lead to different solutions.
You may gain some insight into the subject from Oppenheim himself in lecture 3 of a DSP course recorded in 1975 and now posted online.


How long does it take for snd_pcm_writei() (ALSA) to send data? Using a pretty typical call to Linux's ALSA, snd_pcm_writei(pcm_handle, pcm_buffer, pcm_buffer_size) how do you know how long it should take to send this data? What is the formula?

The sound card should send frames (a frame is two samples for stereo, in ALSA terminology) at a certain rate, however, what is this rate? Sample rate? I don't understand how the sample rate relates to this.

If I have a sample rate of 44100, does it mean it sends 44100 samples per second? Given that my buffer in my program contains 8820 samples, I should see it drain this buffer 5 times every second (44100 / 8820 = 5). However, in reality I see it drain the buffer something like 10 times per second.

I could provide code, but I think it only would make the question less clear.

I've read https://stackoverflow.com/questions/24040672/the-meaning-of-period-in-alsa and www.linuxjournal.com/article/6735?page=0,1#N0x19ab2890.0x19ba78d8 as well as some other guides. However, it is still unclear to me how long it actually takes for the sound card to send data that it has received from snd_pcm_writei().

Please correct me if I'm wrong.

snd_pcm_hw_params_get_period_size() gives the number of times per second a chunk of data is sent, i.e. in Hz.

snd_pcm_hw_params_get_period_time() gives the number of microseconds to send such a chunk.

E.g. for my example, it takes in reality around 10 times per second, and these calls return to me 4410 Hz and 10000 usec, respectively. Hence, this appears correct.

However, it still eludes me how the buffer size plays into this. If anyone who happens to be experienced with ALSA would like to add that, that would be awesome.

Edit:

I believe it works like this:

In the same example, 4410 frames will be required for each period. This makes it 8820 samples. If you're using S16_LE, i.e. short data types, then each frame is 2 byte, hence the buffer should be 17640 bytes.

So, that is the buffer which is treated during each period.

However, there is another important concept in ALSA, the hardware buffer, which can be quite large and introduce a latency when copying it---Hence why ALSA uses the concept of periods at all.

Reference: https://www.linuxjournal.com/article/6735

  " The size of the buffer can be programmed by ALSA library calls. The
  buffer can be quite large, and transferring it in one operation could
  result in unacceptable delays, called latency. To solve this, ALSA
  splits the buffer up into a series of periods (called fragments in
  OSS/Free) and transfers the data in units of a period.

  A period stores frames, each of which contains the samples captured at
  one point in time. For a stereo device, the frame would contain
  samples for two channels. Figure 1 illustrates the breakdown of a
  buffer into periods, frames and samples with some hypothetical values.
  Here, left and right channel information is stored alternately within
  a frame; this is called interleaved mode. A non-interleaved mode,
  where all the sample data for one channel is stored followed by the
  data for the next channel, also is supported. "


Element-wise powering vs. discrete version of binomial theorem on a sum of Gaussians I had recently posted a question on applying powers on a sum of Gaussians (here) to enhance signal resolution artificially. The discussion has given rise to another query. Consider this problem. We have two Gaussian time functions $G_1$ and $G_2$. Suppose we wish sum these Gaussians and raise it to a power of 2.

Situation no. 1: $(G_1+G_2)^2$, we should get another function $G$= $G_1^2$ + $2G_1G_2$ + $G_2^2$. The product of two $G_1G_2$ is another Gaussian. We should see three peaks in this case.

Situation no. 2: We sum the two Gaussians element wise i.e., $S_i$ = $g(t)_{1i}$ + $g(t)_{2i}$ where $g(t)_i$ 's are the corresponding elements of $G_1$ + $G_2$. Now each summed element $S_i$ is raised to power 2. We get only two peaks as a result.

I discussed this issue with some of my colleagues including the person who had proposed this method; they seek agree that both operations are equivalent. It seems they are not, because in situation 2, we will never see $G_1G_2$. What would the discrete version of situation 1?

I feel there is some fallacy here and some mathematical insight should be able to settle this problem.

Thanks.

For Situation no. 2: Here is the MATLAB code.

t=[0:1/200:60]';

u1= 19; % mu the mean time peak i

u2= 22; % mu the mean time peak j

A= 250; %Area

s_ij=[0.84 0.87];; %s_ij represent the sigma of i,j pair; Choose elements by indexing

G_ij=A*normpdf(t, u1, s_ij(1,1))+ A*normpdf(t, u2, s_ij(1,2)); %Sum normpdf is a built in function for a unit area Gaussian from statistical package.

Power= (G_ij).^2;

subplot(2,1,1); plot(t,G_ij);title('Unresolved Gaussians');

subplot(2,1,2); plot (t, Power);title ('Gaussian Sum Raised to Power 2');

([a b c] .+ [d e f]).^2 = [a+d b+e c+f].^2
                        = [(a+d)^2 (b+e)^2 (c+f)^2]
                        = [a^2+2ad+d^2 b^2+2be+e^2 c^2+2cf+f^2]

I suspect that:

The problem may be in your assumptions: "We should see three peaks in this case" sounds like a hypothesis that may be false.
You may have a bug in your code, which might be why "We get only two peaks as a result".

Here is some code to illustrate that (A.+B).^2 = A.^2 .+ 2A.*B .+ B.^2:
t=0:1/200:60;
gauss1 = 250*normpdf(t,19,0.84);
gauss2 = 250*normpdf(t,22,0.87);
m1 = gauss1.^2+2*gauss1.*gauss2+gauss2.^2;
m2 = (gauss1+gauss2).^2;
fprintf('Difference: %1.2f\n', sum(abs(m1-m2)));
subplot(2,1,1);plot(m1);title('Method 1');
subplot(2,1,2);plot(m2);title('Method 2');

I get:
Difference: 0.00

and

So, the math is correct. You may still be wondering, how come the third Gaussian $G_1G_2$ is not apparent in the plot? This is the point where it's useful to challenge your assumption that the third Gaussian should be visible. If we plot $G_1G_2$ by itself we obtain:
plot(2*gauss1.*gauss2);

Look at the amplitude and the position in the x-axis: this third Gaussian is not apparent in the plot because it is completely swamped by the other two, much larger, surrounding peaks.


The transform function of a digital FIR filter The transform function of a digital FIR filter can be characterized by
$$
H(z) = \sum_{n=0}^{2N} h(n)z^{-n}\quad\text{where $2N$ is the order of the filter.}
$$

But when we use $2N$ or $N-1$ , $2N+1, \ldots$, on any basis chosen, is it related to linear phase or  filter characters , or related to various configurations (low pass, high pass, band pass and band cut).
On the other hand , is $h(n)$ the coefficients or weights ?
Can you give me a realistic example of bringing the idea closer?
Are they converted to electrical signals later? Or what?

I Need Clarification Please.

First, we call $H(z)$ the transfer function (not transform).

Next, I believe you should have a look at Wikipedia first, many of your questions shows no preliminary research.

Still, I'll try to clarify few things for you.

I am not sure why your summation uses $2N$, but it should be $N$.

$$H(z)=\frac{Y(z)}{X(z)}=\sum_{n=0}^{N}{h(n)z^{-n}}$$
Where N is the order of the filter.

Looking at a block diagram from Wikipedia should help you understand what this equation means: (Note that $b(n)$ is the same as $h(n)$ here)

To your questions:

$N$ represents the "oldest" delay, which is the same as the filter order. 5 delays makes a fifth order filter. With bigger order, you can achieve more steep rolloff (better attenuation between passband and stopband). All other characteristics are defined by the value of $h(n)$
Coefficients or weights? Aren't they the same thing? Weighting a signal comes down to multiplying it with a coefficients. $h(n)$ are constants value with which you multiply the delayed version of your signal.

Then, your last question is hard to understand, what would be converted to electric signal? In any case, what we consider, is a system with an input ($x[n]$) and an output ($y[n]$). The signal that is fed to that system could be anything, including a electric signal that is converted to digital values.

The coefficients $h(n)$ are merely constants that have no units. DSP theory tells us how to properly select these constants in order to achieve a given behaviour.

Let's make a very easy example with a moving average. The moving average is easy to understand. We could define the moving average like this:

$$y[n]=\frac{x[n]+x[n-1]}{2}$$

If we use this moving average on a noisy signal, we can obtain the following image. Yellow is the noisy signal, blue is the averaged signal.

If we take our moving average equation from above and do a little math we find.

$$Y(z)=\frac{X(z)+X(z)z^{-1}}{2}$$
$$Y(z)=\frac{X(z)(1+z^{-1})}{2}$$
$$\frac{Y(z)}{X(x)}=H(z)=\frac{1z^0}{2}+\frac{z^{-1}}{2}$$
So we got our transfer function. The $h(n)$ coefficients are
$$h(n)=[\frac{1}{2}, \frac{1}{2}]$$

These constants values are what define the moving averager, which is itself a type of low-pass filter. Select other constants and you will find other type of filters.

Last point, phase linearity will happen when the $h(n)$ coefficients are symmetrical. For example : $h(n) = [0, 0.2, 0.4, 0.2, 0]$


How does Correlation for 1D Signals work? On page 2 of the slide, I don't understand the example. Like $j$ and $s$ are both on the x-axis, what does j+s mean? The filter is $(f(−1), f(0), f(1)) = (−0.5, 0, 0.5)$. How can we use it? Could anyone explain how we get $r(j)$ with numbers?

For the simple $3$-tap filter given in the example, each output value $r(j)$ is computed as
$$\begin{align}r(j)&=f(-1)I(j-1)+f(0)I(j)+f(1)I(j+1)\\&=-\frac12 I(j-1)+\frac12 I(j+1)\tag{1}\end{align}$$
So each input value is replaced by the (scaled) difference between its neighboring values.


upsampling in the frequency domain I have the Fourier transform $S(\nu)$. If I oversample by a factor $R$ in the frequency domain. Does the duration of the signal $s(t)=IFFT(S(\nu))$ change?

If by "oversampling", you mean interpolating more frequency samples between your existing frequency samples, then in this case this would have the effect of zero padding your time domain waveform. Therefore the duration of your signal itself doesn't change, but the duration of time over which the DFT is performed would, with the addition of zero padding to fill out the added time duration (which could be before or after the waveform, changing the phase component of the resulting frequency spectrum accordingly).

To understand zero padding it is helpful to understand the difference between the Discrete Fourier Transform (DFT) and the Discrete Time Fourier Transform (DTFT). The Fast Fourier Transform (FFT) is an efficient algorithm for computing the DFT.

Zero padding a discrete waveform in time results in interpolation of DFT result in frequency. The DFT result is discrete, given by samples of the DTFT. The DTFT in contrast is a continuous function in frequency. By zero padding, we introduce more samples on this same DTFT.

Comparing the formulas for the DFT and DTFT will give you further insight into why this is so:

Here the plots show how the DFT is samples of the DTFT, and by adding zeros to the time-domain waveform prior to taking the DFT (zero-padding), we are approximating the continuous DTFT by filling in more samples (interpolation)-- the more zeros we add the closer we can get to the continuous waveform DTFT. So zero padding and interpolating the frequency spectrum are Fourier Transform pairs.

Further notes:  In the plots I say that the DFT is repeating in time, but in actuality it is only given over N samples. The repetition in time is a mathematical equivalence: similar to the Fourier Series Expansion which is defined over a finite time T, it's frequency components only exist at integer multiples of 1/T (discrete in frequency). If you allowed the time domain waveforms of those components to extend to $\pm \infty$, the base waveform that was defined over t=0 to T would repeat in time. The DFT as given above would be the same as the DFT of the same waveform repeating in time, and seeing that helps provide a lot of insight--- anything that repeats in time MUST be discrete in frequency, anything that is discrete in time MUST repeat in frequency (A/D sampling is an example). The DFT does both. Since the DTFT is performed on a single time domain waveform over all time, it never repeats in time and for this reason has a continuous frequency response.


Spectrogram Time Resolution with Overlap I understand that the time resolution of a spectrogram with no overlap is the sampling frequency divided by the window length. How does this change with, say, 50% overlap? Does this double the time resolution? If so, why can't I just have a large overlap and window length to achieve both high time and frequency resolution?

A sliding FFT (maxed-out overlap of N-1) used for a spectrogram will turn a single sample event into a fat wide hump (shaped like the window applied) of width 2N-1.  You might call the fat wife blurry peak of a single hump “resolution”, but multiple overlapped events of different heights plus a bit of noise might make it impossible to “resolve” (as in separate) adjacent event peaks in time.

To make the blurry humps for short events appear narrower you have to use shorter FFT windows, no matter what the overlap.

A 50% overlap does provide more time resolution, as the typical quantized Von Hann window eats edge events.  The extra windows in a 50% overlap brings that info back up into the picture, as the sum of 50% overlapped Hamming or Von Hann windows has a constant gain.


Why integrate over $2\pi$ in inverse DTFT? In DTFT of a signal, the spectrum of a sequence is periodic with period $2\pi$ and all the information needed for derivation of the original signal from its spectrum is contained in $\pi <\omega <\pi$.

But , why do they integrate over $2\pi$ and not from $-\infty$ to $+\infty$ , as in continuous Fourier transform in inverse DTFT? After all, when you derive a transform of a signal, the inverse of it is not arbitrary and to recover the original signal, you can't integrate over an other  interval because all the information is contained in there.

There are two reasons:

As you say, the sequence is periodic so all of the information is contained in the $-\pi < w <\pi$ region.  Thus, integrating over that region tells you everything that there is to be learned from the inverse transform.
Integrating from $-\infty < w< \infty$ would only give answers of $0$, $\infty$, or $-\infty$.  The reason for that is that integrating over $-\pi < w <\pi$ will generally give a finite answer which can be something negative (e.g. -4), 0, or something positive (e.g. 2.5).  If you integrate over $-\infty < w< \infty$ each $2\pi$ region will give the exact same answer, so the negative results will accumulate to $-\infty$, the zeros will sum to $0$, and the positive results will sum to $\infty$.  (Note that the results are often complex, but that doesn't change the primary point.  It just increases the dimensionality of the infinities.)