How do optical illusions work (photo mosaic) from a signal processing perspective? Hello fellow investigators

I have two question about optical illusions

1) A photo mosaic is something like this:

What are the signal processing principles behind our eye merging the many tiny images into a large coherent one?

2) Are there etiological explanations of optical illusions in general - by this I mean can it be explained why image inputs get parsed in 'wrong' ways by the eye/brain, why does the human parser get confused; what transfer function is in play?

Thanks for your time

What are the signal processing principles behind our eye merging the many tiny images into a large coherent one?

From what I've read (Source 1: Front-End Vision and Multi-Scale Image Analysis, Source 2: Information Visualization), our eyes apply something very similar to a set of gabor filters to the visual inputs we get. So if you, i.e. focus on Mona Lisa's eye, your brain processes gabor filters at the scale roughly the size of the eye, and such a filter removes details at smaller scales.

  Are there etiological explanations of optical illusions in general - by this I mean can it be explained why image inputs get parsed in 'wrong' ways by the eye/brain, why does the human parser get confused; what transfer function is in play?

Regarding optical illusions in general: There is tons of research. The two books I've mentioned above quote some of it, but I don't know enough about it to compress it into a good answer here.

But regarding this special case, there's actually a very elegant explanation in Front-End Vision and Multi-Scale Image Analysis, p208, which I'll try to repeat here: Take a "natural" sample image. Split the image into small blocks. Multiply each block with a window function (e.g. a gaussian).

For example, I'll take the "Lena" image in Mathematica:

sample = ColorConvert[ExampleData[{"TestImage", "Lena"}], "Grayscale"];

split it into 16x16 pixel blocks:

blocks = Flatten[Partition[ImageData[sample], {16, 16}, {5, 5}], 1];

multiply each block with a gaussian window:

window = Table[Exp[-(x^2 + y^2)/25], {x, -7.5, 7.5}, {y, -7.5, 7.5}];
blockTimesWindow = Map[window * # &, blocks];

the blocks look something like this:

Next I'll use SVD to find the principal components of these blocks:

{u, s, v} = SingularValueDecomposition[Flatten /@ blockTimesWindow, 10];
Image[0.5 + Partition[#, 16], ImageSize -> 100] & /@ Transpose[v]

and visualize them:

These are the principal components of small parts of my test image. Or in other words: any small part of my test image is (approximately) a weighted sum of these principal components, plus noise. And it's not just this test image. Take any photograph of a natural scene, the result will look very similar.

As you can see, the principal components look suspiciously like gabor filter kernels or gaussian derivatives. That explains why filters like these are a good choice to analyze natural images: in a natural image, any small part is (almost) the weighted sum of a few gabor filter kernels plus noise. And it explains why evolution would choose this kind of low-level processing.

Now your image is not a natural image: Neighboring pixels are not correlated the way they would be in a natural image. But our eyes still use processes that were optimized for natural scenes, where differences at smaller scales don't carry as much important information, so they are filtered out. That explains why the illusion works.


Why analog anti aliasing filter is used before analog to digital converter when there is already a digital filter after ADC? Normal data acquisition consist of:

Analog anti aliasing filter( Sampling frequency : $5\textrm{ kHz}$)
ADC - Digital Filter - (Sampling : 200K samples /sec)
Digital low pass filters Filters
DAC

Questions:

My question is why analog anti-aliasing filter is used when their is already a digital low pass filters after ADC to prevent anti aliasing.
If analog anti-aliasing filter have sampling frequency $5\textrm{ kHz}$ , the system will not take frequency greater than $2.5\textrm{ kHZ}$,

then why ADC frequency is $200\textrm{ kHz}$?
Doesn't analog anti-aliasing filter  sampling frequency limit overall system frequency ?

Question 1:
The anti-aliasing filter before the ADC is exactly for the purpose of rejecting high frequencies, that will become lower frequencies (i.e. aliasing) after the ADC. The digital lowpass after the ADC cannot help here, as the aliasing has already happened. Consider this example:

Your ADC has a sampling frequency of Fs=100kHz.
Your input signal is a sum of two sine waves, with frequencies 10kHz and 220kHz.
After ADC you would find two sine waves: one at 10kHz, one at 20kHz (220kHz-2*Fs).
Hence, you have aliasing occured, and no digital lowpass can remove this aliasing.

Question 2:
Without more information on the system this cannot be answered. However, here are some thoughts:

filter of 5kHz requires a sampling frequency of at least 10kHz (ideally). You state you only need 2.5kHz. I think you mix something here.
in reality, no anti-aliasing filter is a perfect low-pass, hence its cutoff-frequency does not mean, that higher frequencies are perfectly blocked. Instead, they are more and more attenuated. To cope with non-ideal anti-aliasing filters, the sampling frequency should be higher than 2 times the cutoff. However, 200kHz from your example occurs still quite high for me.


Designing High Pass Filter to Remove DC Offset from a Signal I want to design a high pass filter (f0*S/1+f0*S) in order to remove the dc offset from a signal ? (i.e. how to choose the frequency f0).

Thanks,
Hacene

Ps : the filter is for MATLAB implementation, this is why i would like to have it in continuous time.

If the DFT is the Uniform sampling from $ 0 $ to $ 2 \pi $ then the first bin is given by:

$$ x[k] = \sum_{n = 0}^{N - 1} x[n] $$

Namely it is the sum of all the samples.
Hence in order to remove the DC (Mean) all you need is a filter which has zero in its DC bin.
Since, the filtered signal, which is a convolution (Circular) of the DFT of the input signal and the filter will have zero at the first bin which means the sum of the output is zero which means its mean is also zero, as wanted.

Simple and intuitive FIT would be an FIR with the length of signal which removes from the current sample the mean of all samples:

$$ y[m] = x[m] - \frac{1}{N} \sum_{n = 0}^{N - 1} x[n] $$

This is a simple FIR.

Pay attention I assumed you have all the samples.
If you don't, you need to make this FIR casual by only "Calculating" the sum of given samples.


Downsampling impact on complex phase For my application, I have to downsample a bandpass complex signal which spectrum is located on the second Nyquist zone.
Knowing that this processing will cause a spectrum inversion, what would be the additional side effects (SNR reduction, less amplitude ..) ?
Will the phase's linearity be conserved (supposing that the pass-band anti aliasing filter has a linear phase) ? Does anyone have a reference detailing those side effects of the downsampling ?
Thanks for any help.
Mourad.F

Decimating a signal (selecting every Dth sample and discarding the rest) does not distort the signal within the passband in any way other than to cause aliases from higher frequencies to fold into the signal bandwidth.  Depending on how we model the system the phase may be effected since $z^{-n}$ is replaced with $z^{-n/D}$, but the phase will still be linear, so the signal itself will not be distorted. (In multi-rate control loop modelling applications for example I had to pay attention to this change in phase slope). For further information on this side detail see the last paragraph.

Therefore an "anti-alias" filter is required, no different than what is done in and A/D application, just in the digital domain (decimation is a digital-to-digital converter).

The ideal decimation filter would pass the passband with no distortion while completely rejecting all the energy in the image bands. If that could feasibly be done then the decimation operation itself would be perfect (no effect on the signal in the passband). Actual performance depends on how well we design the decimation filter.  So in your case, of a signal in the second Nyquist zone (that once sampled will be in the first Nyquist zone with spectrum inverted as you describe) will not go through any further distortion with a subsequent decimation (if linear phase filter is used the result will be linear phase).

Since the folding areas are in distinct bands, a multi-band filter is a good choice for decimation filter design (firls and firpm in Matlab support the design of multiband filters).

As far as the impact on noise due to decimation itself, see this post that further details how the variance of the noise is not effected by the decimation: Properties of Up- and Downsampling

However the devil is in the details, as with any digital design you must pay attention to quantization effects, filter realization and passband distortion etc but this is not due to decimation specifically.

To further explain what I touched on in the first paragraph regarding the change in phase slope, consider a unit delay prior to a decimate by 2. A unit delay (with transfer function $z^{-1}$) has a constant magnitude over all frequencies, and a phase shift that varies from 0 to $2\pi$ as the frequency goes from DC to the sampling rate. After the decimate by two, the transfer function of the unit delay becomes $z^{1/2}$ which has a constant magnitude over all frequencies with a phase shift that varies from 0 to $\pi$ as the frequency goes from DC to the new sampling rate. So in the digital domain where we work with normalized frequencies, the slope of the phase has changed, but with regards to the related analog frequency in Hz, the slope is the same: At our new sampling frequency in Hz we see the same phase shift we got at that particular frequency prior to decimation.


A question about Gibbs phenomenon example in Richard Lyons book If anyone has a copy of this book could you please shed some light here: I was trying to reproduce the example shown in Chapter 5, section 5.3.1 Window Design Method's Figure 5-19 which illustrates the ripples in the frequency domain in octave. Just as the book suggested, when you have any dropped tap in time domain, you will see Gibbs phenomenon in the frequency domain. However, my question is that, why don't you use all 32 taps in this case (the book is dropping 1 tap from 32 taps to get 5-19(c) )?  When I use all 32 taps, I can see the un-distorted frequency domain here since the iDFT and DFT pretty much cancelled out each other here...I think I am missing some points here the book is trying to make...

Thanks,

Will

I don't have the book, but from what I understand from your question, you take a discrete desired frequency response with $32$ elements, apply an IDFT to get a $32$-tap filter, and then you remove taps to see the effects of truncation in the time domain. Well, this is not exactly what the window design method is about, and it is no surprise that if you don't remove filter taps in the time domain and apply a DFT that you get back your desired $32$-tap frequency response (because, as you've figured out yourself, IDFT and DFT are inverse operations of each other).

However, your actual frequency response does not have $32$ values in the frequency domain but it is continuous in frequency. It is given by the discrete-time Fourier transform (DTFT):

$$H(e^{j\omega})=\sum_{n=0}^{N-1}h[n]e^{-jn\omega}\tag{1}$$

where $N$ is the number of time domain taps (filter length), and $h[n]$ are the filter coefficients. You can compute a sampled version of $(1)$ by using the DFT, and for example, if $N=32$, you can (and should) compute more than $32$ frequency domain samples by using zero-padding. In Matlab/Octave, you can use the commands fft or freqz to compute sampled versions of $(1)$ with an arbitrary number of frequency domain samples.

To help you improve your understanding of the window method, and of DSP in general, I would like to recommend to you a really good (and free) book: Introduction to Signal Processing by S.J. Orfanidis.

Also take a look at this answer, this answer and this answer to related questions.


How to implement filter poles in noise shaping for ΔΣ modulator? I am currently stuck trying to implement filter poles to stabilize a ΔΣ modulator. What I think I know is:

a noise shaping filter with only complex zeros along the unit circle destabilizes a 1 bit ΔΣ modulator when the order is beyond 2.
filter poles can be used to reduce the out-of-band noise, making the modulator more stable

For testing this approach, I am trying to realize a simple 2nd order modulator with 2 zeros and 2 poles.
Initially, I tried it without the poles and only with the two zeros which works and I implement it in the following way for testing (using python/numpy):
import numpy as np
roots = np.exp(1j*np.array([0.05, -0.05]))                  # two exemplary complex roots
b = np.poly(roots)                                          # and corresponding filter coeffs
errors = np.zeros_like(b)
signal = np.random.rand(100) - 0.5                          # any signal smaller than [-1,1]
samps = np.zeros_like(signal)
for iSAM in range(len(samps)):
    desired = signal[iSAM] + np.dot(b[1:],errors2[:-1])     # add the filtered errors to the input
    samps[iSAM] = np.sign(desired)                          # discretize to 1-bit
    errors = np.roll(errors, 1)
    errors[0] = samps[iSAM] - desired                       # calculate new error

As expected when FFT'ing the results, I see a frequency response like that of the desired filter.
Also as expected when adding another root at exp(i*0)=1 (creating a 3rd order modulator), the whole thing becomes unstable, even for small inputs.
I have no issue making filters with poles, and defining their frequency response. I am also able to take a time domain signal and filter it, with the designed biquad filter (2 zeros, 2 poles) and I obtain correct results.
The question is: how do I use the poles in the noise shaping algorithm. How would I even use a simple biquad filter for noise shaping ? I have watched many videos and papers today, but still don't understand it. Everyone draws these diagrams like this here from a video series i watched: https://youtu.be/IE8tU_10Hpg?t=63
A more detailed flow chart is in Stanley Lipshitz paper from 1991.
However, even if i follow this flow chart literally, I get no results. The state variables just explode. The guy in the video series also implemented a 3rd order filter with 3 zeros and 3 poles to define a stable 3rd order modulator. While I can reproduce the filter coefficients, and the Noise transfer function (NTF), I don't understand how to actually use the filter. In one video he makes a quick remark about testing the filter performance in a simulation. He says that performing the test is "trivial" or something like this. So obviously I am missing something very fundamental here :)
UPDATE:
After more days of trying, I am really at a loss. Can anyone please post an example of noise shaping with even a very simple IIR filter, even with only 2 zeros and 2 poles. The shaped noise spectrum should reproduce the magnitude response of the designed filter.
I am putting some bounty on it as much as I can. Because all the guides on the internet haven't helped me. I don't care about which programming language you use for it but the shaped spectrum should demonstrably produce the correct magnitude.
Best regards

I believe I have actually found an answer myself after even more research ;)
The crucial piece of information is the following: The filter function H(z) shown in many models is not the actual Noise Filter Function. Usually the noise filter should be some kind of high pass and I used this directly as H(z) erroneously.
Instead, the noise is shaped by the noise transfer function NTF(z) = 1-H(z).
Therefore, after designing the highpass noise transfer filter NTF(z), one needs to calculate the filter model via H(z) = 1-NTF(z).
It turns out that, if the NTF is a "poleless" filter (FIR filter) of form NTF(z) = 1 + Sum(b_n z^-n), then the filter function will be very similar: H(z) = 0 - Sum(b_n z^-n) where n starts from 1. The first coefficient will be 0.
I actually implemented this noise shaping directly in the past, by ignoring the leading coefficient 1 and inverting the sign of the filter result at the summing node with the signal. I came to this conclusion, because I picked it up somewhere but didn't understand its origins.
When the NTF contains poles, however, the H(z) will deviate more significantly. The denominator of NTF(z) and H(z) will be in fact identical. But the numerator will be X(z) = Sum(a_n z^-n) - Sum(b_n z^-n), where a_n are the pole coefficients, and b_n are the zero coefficients. As typically, the leading coefficients of the NTF(z) will be both 1 ( it is a minimum phase filter after all), the leading coefficient of H(z) will become 0. This produces the missing first coefficient naturally.
With this properly defined H(z), following the flow chart as given e.g. in Lipshitz 1991 paper, gives correct results. I used +/- 0.5 of the quantizer range as uniform white dither.
Update:
after correcting some minor mistakes, it all works as it should. In Python I use the Inverse Chebyshev Filter model to make the Noise Transfer Function. Then I correct the coefficients, to turn them into the H(z) filter model:
fDAC=20e6
b, a = signal.cheby2(5, 160, 20000, 'hp', fs=fDAC)
b /= b[0] # normalize according to Gerzon/Craven
b = a - b # turn NTF(z) into H(z) for noise shaping

Then the actual noise shaping (left some redundancy in the code, to keep it closer to Lipshitz 1991):
l = 200000
ts = np.arange(l)/fDAC
wav = 0.55 * np.sin(2*np.pi*3200*ts)
errbuf = np.zeros_like(b)
samps = np.zeros_like(wav)
for iSAM in range(l):
    desired = wav[iSAM] - np.dot(b, errbuf)
    samps[iSAM] = np.sign( desired )
    error = samps[iSAM] - desired
    errbuf[0] = - error
    errbuf[0] = - np.dot(a, errbuf)
    errbuf = np.roll(errbuf, 1)

Note that no additive dither is being used. Even though the quantizer noise is not white itself, the noise shaping is efficient enough to remove the quantizer distortion. Undithered Delta Sigma modulators do have drawbacks tough, e.g. on stationary input, as can be read in numerous sources. If you want dither, it should be added right before quantization, i.e. inside the sign() function

Figure: Example of a 5th order delta sigma modulator. Vertical scale unit is full scale. Horizontal is Frequency in Hz. The orange curve is a designed 5th order inverse Chebyshev Highpass filter as Noise transfer function. Blue is the result with a ~-5dBFS sine wave.
PS: Since I originally set a bounty on this question, which was keeping me up at night: Please feel free to correct me and claim the bounty with a better answer of your own.


Frequency response of a rolling linear regression I am looking for a way to characterise the frequency response of the slope from a linear regression. We are exploring the effect of window length of the regression to the magnitude of the slope of the regression and find that for our data (discrete time series data), the slope is positively related to the window length.

The slope can be thought of as the expected change of our output variable for a given input variable over a time period equal to the window length. We find that for a short window length, this slope is very low whereas for a long window length the slope is high. We interpret this to mean that the system responds to changes in the input in the long term but not in the short term. To me, that sounds a lot like a low pass filter killing off the fast changes in the input but passing the slower (or rather, more sustained) changes.

The slope, or beta, of the regression for a rolling window is, in a way, similar to a moving average and it is possible to characterise the frequency response of a moving average filter as a function of its window length as shown here. I was wondering if anyone knew of a similar way to characterise the frequency response of the beta of a rolling regression as a function of its window length in a similar manner?

Note that I have read this paper which explores using regression beta as a filter but I was unable to get an intuitive understanding of why the larger the window length is, the lower the cutoff frequency appears to be.

Here are a few descriptive charts of our data that might show why we are asking this question:

In the chart above, the thick grey line is the regression over all the data (i.e. a long window length) and has a slope of 1.79. The coloured lines (and points) show the regressions for single years of data (the betas for these lines are in the brackets in the legend of the chart). The beta for any individual year is significantly lower than the overall beta. In other words the grey line is much steeper than the other lines. This is true for any arbitrary 1 year window on the data (i.e. there is no year-long window of data that has a slope as steep as the overall data).

This next chart shows the rolling betas for different window lengths:

So the question we are trying to answer (and hoping to do so by characterising regression as a LPF) is why do the lines in this chart have higher biases (higher on average) as the window length increases (e.g. the yellow line is higher than the green line). How can we interpret this phenomenon?

This is a statistics question and should not be asked here. But I will try to answer the question as best as I can (it won't be the best answer). First off, I don't see what the phenomenon you are observing has to do with frequency response or LPFs. Low pass filters makes data smoother, which is not what you have. Instead you have a situation, where estimate for one set of data differs greatly from a other set of data. The OLS estimator also is not an MA process, the formula being cov(x,y)/var(x). The filter described in the paper is not a linear time invariant filter.

Indeed, if you regress a line on white noise process or linearly growing process, the regression line should not change slope significantly as more data is gathered. Instead as the amount of data grows you should be more certain of that the coefficient is right, generating greater variance for the rolling process, but on average the slope should be right.

What you have here is not a process that satisfies Gauss Markov assumptions. As can easily be seen the expected value of the error is not zero, and there is heavy auto correlation. This means that the error term is correlated with the independent variable, producing a biased estimate. Besides that the data is not randomly selected.

Most importantly, the data is increasing as a function of time. This generates a so called spurious correlation. Since Y increases faster than X as time passes, the beta becomes upwardly biased. You can model autocorrelation with an MA process, but this is not really important...

I am sure statistics stackexchange (crossvalidated) would provide you with much more mathematically precise answer and a better explanation of what the phenomenon. One thing fix that might help is to difference both of the variables, which changes the autocorrelation structure of the error. This can be equally viewed as a (1-1) MA process applied to both of the variables, which is a comb filter at a high frequency, producing a high pass filter below the nyquist frequency. But the data has such a non-standard structure that I wouldn't be so sure that the fix will help you.


How to estimate the mean of a clipped signal? An accelerometer is used on a high vibration environment (a bicycle) to measure  gravity. The signal clips, so data higher than a clipping value is clamped to that value.
So, simply filtering the signal gives the wrong median mean value when the vibration is high.

I want the mean value over a time much longer than the vibration frequency, so a few second filter is fine. I do not care about the vibration waveform itself. Just want the average value. How to remove the error caused by the clipping?

Bad idea one:
I had suggested that the signal's noise distribution be modeled, the amount of clipping predicted based on measured vibration amplitude, and that used to correct the data. Maybe the distribution is Rayliegh, maybe Gaussian, maybe... So maybe that first idea is not a good one.

Possibly brilliant idea two:
But how about just measuring the distribution on the spot? Make a histogram with several bins, accumulating a count in a bin each time the sampled data is in a bin's range.

Then, fit a curve over the top 3 or 5 highest bins, and the position of the peak is the signal median the modal mean, where the signal mean value would be, if it were not clipped.

Is this crazy? Good? Done before and has a name?

Is this crazy? Good? Done before and has a name?

I wouldn't say that it is "Good" because clipping is not noise that can be modeled as an additive component with some distribution.

  fit a curve over the top 3 or 5 highest bins, and the position of the
  peak is the signal median

The position of that peak would be the modal mean, not the median

The effect of clipping on the signal's distribution would be to a) truncate the real distribution (whatever that is) at the maximum (which depends on the sensor & measurement (or "environment" as you put it) and b) distort the distribution towards this maximum depending on how much clipping was present in the signal.

b) means that whatever you do, your median / mean will be biased by the amount of clipping.

Ideally, you would have to change the sensor to one that can respond to higher accelerations.

If you absolutely must use the current sensor and deal with the clipping then prior to any further processing of its output you could try removing the clipped values (with a simple thresholding operation) and fitting a curve to the data. If you were to fit an spline for example, then the curve would depart at the beginning of the clipped part with the slope of the usable signal prior the clipping and arrive at the end of the clipped part with a slope of the usable signal after the clipping. In this way, all clipped parts will be converted to some "soft" peak (because of the splines "smoothing" effect).

Of course, this is a compromise. It could well be that the clipping is so much that very little "usable" signal remains on either side to fit a "good" spline. Furthermore, depending on your sampling frequency ($Fs$), the signal might be full of impulse-like data which would lead to "very tall" splines as the algorithm would try to fit the sharp slopes. More information would be needed to take a decision on that but in general a higher $Fs$ would help.


Band-limited random signal with arbitrary distribution? I'd like to generate a random discrete-time signal that is band-limited to some bandwidth B (by means of a digital filter, ie in MATLAB). The catch is that I'd like this signal to have an arbitrary PDF, for example, a uniform distribution.

If I start with a uniform distribution, and then filter to band-limit the signal to B, I end up with something that's basically gaussian. Since my band-limiting FIR filter has a large number of taps (say 128), this makes sense given the Central Limit Theorem, as I am basically summing a number of IID random variables.

What I don't understand, is if it is possible to now transform this signal to a new distribution, such as uniform, while maintaining the band-limitedness. I understand that I can transform a RV from one PDF to another by basically integrating/mapping over the CDF (hand-wavy, it's been a while!), but I believe this significantly changes the frequency content of the band-limited signal. For example, if I have a gaussian distribution, and I want to transform to a uniform, I basically map values "near the mean" of the gaussian to something "closer to the max" in the uniform RV. But now the frequency content has changed, and my spectrum will be smeared.

Does what I'm describing make sense, and has anyone looked at something similar? Perhaps what I need to do is look at transforming the signal as a vector, not a scalar (ie treat it as a random process, not individual random variables). But this is where I get a little out of my depth.

Consider this approach:

General a white Gaussian random sequence.
Filter the white Gaussian sequence.  The output, which we'll call $z$, will be Gaussian because a linear combination of Gaussians is a Gaussian.
Pass the output $z$ from stage 2 through the transformation $g(z)=F_Y^{-1}(F_z(z))$, where $F_Z(z)$ is the Cumulative Distribution Function (CDF) of $z$ (the Gaussian CDF) and $F_Y^{-1}(y)$ is the inverse CDF of your target non-Gaussian distribution.  The output $y$ will now have the desired non-Gaussian distribution.

The only question is how to choose the filter taps in stage 2 so that the power spectral density at the output of the nonlinearity in stage 3 is what you want.  The sloppy but simple and often effective answer is that stage 3 typically doesn't change the power spectral density much, and you can often get away with just filtering $z$ so it has the desired power spectrum.  The more rigorous approach is to note that the autocorrelation of $y$ at lag $\tau$ $R_y(\tau)$ depends only on the autocorrelation of $z$ at the same lag $R_z(\tau)$.  Thus the problem of finding the needed $R_z(\tau)$ to achieve the target $R_y(\tau)$ reduces to separate optimization problems for each lag $\tau$.  It also turns out that $R_y(\tau)$ is a non-decreasing function of $R_z(\tau)$, so you can solve for $R_z(\tau)$ numerically fairly easily using very crude techniques.  Then, once you have the needed $R_z(\tau)$, you can design a filter to operate on the Gaussian sequence to achieve it.

Justification for the above statements can be found by searching for the following in Google Scholar: Cario & Nelson, Modeling and Generating Random Vectors with Arbitrary Marginal Distributions and Correlation Matrix


Best way of segmenting veins from arm? I have done fews researches and found out approach such as

multiscale and
region growing,

that could be used of extracting forearm veins.

However, I tried to use frangi filter and Figure 1 is the segmentation result using the best parameter set I could found;

Figure 1

Could anyone suggest what approaches as well the technique to adopt while considering this problem?

I give extra image for further analysis:

So one good step to enhance the vein-like structures is coherence enhancing diffusion:

  Weickert, Joachim. "Coherence-enhancing diffusion filtering." International Journal of Computer Vision 31.2-3 (1999): 111-127.

So I first apply this algorithm to your image, aggressively. The next step is to identify the curvilinear structures, which would in this case extract the veins. I already explained how to do this robustly in a previous post. So please check this.

The result is the following:

One could apply all the desired noise cleaning to this image. Note that the edges appear very weak, so it is better to tune the parameters, probably even more than I did.


Is multirate signal processing involved in ADSL? Telecom companies offer different data rates for home/office purposes depending on the money they pay.

How are the telecom people able to control the flow of data of different cables?
Say supposing a cable can be used till $20\textrm{ MBps}$. But he paid only for $5\textrm{ MBps}$ scheme, why is the user not able to steal the remaining bandwidth?
Is there a multi rate theory behind this?
Or it is a simple filtering?

"The telecom people" are the ones providing the internet data to you. Hence, they can easily throttle data at their end, making sure that you do not get faster access than you pay for. If the actual signal processing is adapted based on this I would say depends on the solution.


Is there a way to break down a signal into its respective phases (4 in this case) to be later analysed individually Is there a way to extract amplitude information but over 4 phases of the same frequency. That is, if the frequency is known could I extract amplitude peaks (between given values)in the time domain in lots of four, then separate them into "bins" or "arrays" ready for four individual FFTs?

This would be trivial to do as simple decimation where each block increments its samples by $n+4$ with each block starting at sample 0,1,2,3 respectively. This is common with polyphase filter implementation and similar techniques to reduce the overall clock rate requirement for the processing (parallel processing).

For more details on both of those approaches see these posts:

How to implement Polyphase filter?

4-phase FIR LPF in MATLAB

Why does this parallel IIR filter have such a weird form?

Another straightforward and interesting approach specific to DFT processing is to use the Decimation In Time binning that is inherent in the (DIT) FFT algorithm itself.

See the block diagram below for further detail on how an 8 point DFT can be binned in time into 4 2 point DFTs, and specifically how to combine the resulting DFT outputs for any size data block.

Expanding this to any general block length of radix 4 length (a length that is a multiple of $2^4$ so that it can be divided into 4 identical blocks. For blocks not in increments of $2^4$ simply pad with zeros to get to this length), the input sequence is binned by two decimation in time operations each with a factor of 2, best demonstrated by example:

Consider a sequence $x[n]$ numbered 0, 1, 2, 3, 4, 5, 6, 7, ...

After one decimation in time by 2 we divide this into two groups, noted here as Group one and Group two:

Group one (even samples): 0, 2, 4, 6, 8, 10 ....

Group two (odd samples): 1, 3, 5, 7, 9, 11 ...

Each of these is then decimated by 2 to form each of the four separate sample blocks, as follows:

Group $x_1[n]$ (even samples of Group 1): 0, 4, 8, 12, ...

Group $x_2[n]$ (odd samples of Group 1): 2, 6, 8, 10, ...

Group $x_3[n]$ (even samples of Group 2): 1, 5, 7, 11, ...

Group $x_4[n]$ (odd samples of Group 2): 3, 7, 11, 13, ...

The reason for binning in this fashion is for the efficient recombining of the DFT results for each block using 2 point DFT "butterflies" as shown in the graphic above and further detailed below.

The four DFT results are combined to be equal to the DFT of the whole sequence as follows:

Given each $N/4$ length DFT result of the four groups as $X_1[k]$, $X_2[k]$, $X_3[k]$, and $X_4[k]$, where k is the integers from $0$ to $N/4-1$ representing the $N/4$ samples in each group:

Combine $X_1$ with $X_2$ to form a $N/2$ length DFT of the combined result. $k$ is incremented from $0$ to $N/2-1$ and due to the cyclical nature of the DFT, $X_1[N/4]=X_1[0]$, $X_1[N/4+1]=X_1[1]$ etc
$$X_{1,2}[k] = X_1[k]+ W_{N/2}^{k/2} X_2[k]$$

Similarly combine $X_3$ with $X_4$ to form a $N/2$ length DFT of the combined result.

$$X_{3,4}[k] = X_3[k]+ W_{N/2}^{k/2} X_4[k]$$

Finally combine the two resulting $N/2$ length DFT's to have the full $N$ length DFT result, here k is incremented from $0$ to $N-1$:

$$X[k] = X_{1,2}[k]+ W_{N}^{k} X_{3,4}[k]$$

The notation above $W_n^k$ are the "roots of unity", coined as "Twiddle Factors" by W.M Gentleman and G. Sande in their 1966 paper "Fast Fourier Transform - for fun and profit". The twiddle factors are simply convenient shorthand for the exponential expression below and represent phase rotations:

$$ W_N^k = e^{j2\pi/N}$$

Below shows a graphic example for N = 8


Is the filter described by specific difference equation time independent Is the filter described by
$$y[n] + 2 y[n-1] + y[n-2] = x[n] - 2 x[n-1]$$
time invariant?
($x[n]$ - input, $y[n]$ - output)

Theorem of Foures-Segal says no (?)

When I asked about Foures-Segal I meant the following.
F-S theorem says that if transfer function $T$:

$$y=Tu$$

is time invariant then there exists a function $G(z)$
which is analytic and bounded in exterior of unit disk
on complex plane,
such that

$$T = Z^{-1} G(z) Z$$,

where $Z$ is the Z-transform.

In my case the function $G(z) = \frac{z(z-2)}{(z+1)^2}$
having pole at $z=-1$, is not bounded. I would say
it means the system is not time invariant?

From the given difference equation, nothing can be said about linearity or time-invariance. Note that the difference equation (DE) doesn't even uniquely specify the output for a given input $x[n]$. To see this, assume that $y_1[n]$ is a solution to the DE. Let $y_h[n]$ be a solution to the corresponding homogeneous equation

$$y[n]+2y[n-1]+y[n-2]=0\tag{1}$$

Then any sequence $y_2[n]=y_1[n]+y_h[n]$ must also be a solution to the original DE. Note that if $y_h[n]\neq 0$, the corresponding system is neither linear nor time-invariant because that part of the output signal does not depend on the input $x[n]$. So no matter if you shift the input signal, or if you multiply it with a constant, that term remains unchanged.

In order to uniquely specify $y[n]$ for a given input $x[n]$, we need auxiliary conditions. There are several options, but the auxiliary conditions that make sure that the system described by the DE is linear and time-invariant (and causal) are initial rest conditions. If we assume that $x[n]=0$ for $n<0$, then the initial rest conditions are $y[-1]=y[-2]=\ldots=0$. This makes sure that $y_h[n]=0$ for $n\ge 0$, and consequently, the output $y[n]$ contains no term that is independent of the input $x[n]$.

A system described by a linear DE with constant coefficients (such as the one in your question) describes a linear and time-invariant (and causal) system only in combination with the aforementioned initial rest conditions.


Why does spectral leakage arise in an FFT? I've tried googling and wikipedia-ing it, but I haven't gotten any answers beyond 'it's because the frequency of the input signal is sitting between two bins'.

I understand that this is the reason, but what I can't understand is why the leakage seems to extend to several adjacent bins rather than only one adjacent bin.

To illustrate what I'm talking about, here's some simulated data (code at the end of the post) :

Above is the FFT spectrum (plotted on a log scale) of a sine wave of frequency 10. The sampling rate is one, and the number of samples is 100. The graph has been FFT-shifted. There's clearly only a peak at bin 10, and the rest is on the order of numerical error, or there about.

This is the frequency spectrum at a generated frequency of 10.1.
Clearly there is 'leakage' into more bins than just the immediately adjacent bin.

This is the plot for a frequency of 10.5.

Question: Why is there this leakage, and why does it extend to all the other bins, rather than the immediate adjacent bin?

Code, for anyone who's interested (Python code)

import numpy as np
import matplotlib.pyplot as plt

xFreq = 10.5
xSize = 100.0
xPeriod = xSize/xFreq
x = np.linspace(1,xSize,xSize)

data = np.sin(2*np.pi*x/xPeriod)
fft = np.fft.fft(data)
fft = np.fft.fftshift(fft)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(abs(fft), "o")
ax.set_yscale('log')
plt.show()

I changed the xFreq value from 10.0 to 10.5, etc.

An FFT has finite length, and thus constitutes a default rectangular window on a data stream.  A window in the time domain results in a convolution in the frequency domain with the transform of the window.  Note that the transform of a rectangular window is a Sinc function (sin(x)/x), which has infinite width.  It's not just 2 bins in width.  Thus the ripples of the Sinc function will show up as "leakage" far from any spectral peak that is not perfectly periodic in the FFT's length.

The picture below shows part of the frequency response of the sinc function.  When the tone is centered on one of the bins all of the other points are lined up with the nulls in the frequency response.  If it isn't centered on a bin then it's like shifting the whole frequency response, which causes the other bins to fall on non-null portions of the frequency respoonse.

Another way to look at it is that an FFT is just a filter bank, where each filter stop-band floor has lots of ripples, and is certainly not infinite in attenuation more than 1 bin away from center frequency.  Some windows (von Hann,etc.) other than rectangular have lower stop bands, which is one reason for their popular use.


Derivation of Bessel filters I was reading Winder's Analog and Digital Filter Design and the section on Bessel filter.

I was hoping to see a complete derivation of the Bessel filter theory, but Winder's book gives only

\begin{align}
H(s)&=e^{-sT} \tag{$\scriptstyle\text{pure delay}$}\\
H(s)&=e^{-s}=\frac{1}{\sinh(s)+\cosh(s)} \tag{$\scriptstyle\text{normalized pure delay, T=1}$}
\end{align}

and then using series expansion of $\sinh$ and $\cosh$ we're left with

$$H(s)=\frac{a_0}{B_n(s)}$$

where $B_n(s)$ is a Bessel polynomial.

Anyone know where can I get a more exhaustive derivation of the Bessel filters?

The optimality criterion for Bessel (low pass) filters is that they have a maximally flat group delay at $\omega=0$, as correctly pointed out in a comment by Peter K.. The derivation is slightly tortuous, but you can find it for the continuous-time case in these lecture notes, which are taken from the book Introduction to Circuit Synthesis and Design by Temes and La Patra. I can recommend this book to anybody who wants to gain a deep understanding of electrical network theory.

The interesting fact is that unlike for other standard filters (Butterworth, Chebyshev, Cauer), the discrete-time version of a Bessel filter cannot be derived from its continuous-time counterpart via the bilinear transform. The reason is the frequency warping inherent in the bilinear transform, which doesn't affect the optimality conditions of the other standard filters (maximal flatness of the magnitude, or equiripple behavior of the magnitude, or a combination thereof), but it does not preserve the maximally flat group delay. A discrete-time filter derived from an analog Bessel filter via the bilinear transfer will not be optimal in any way, and it will not have a maximally flat group delay. This means that a discrete-time Bessel filter must be derived directly in the discrete-time domain. This has been done in this famous paper by Thiran. Again, the derivation is slightly painful but I'm afraid there is no easier route.

What you can take away:

Bessel filters have a maximally flat group delay.
Discrete-time Bessel filters cannot be obtained from continuous-time Bessel filters via the bilinear (or any other) transform.


Circular Convolution Matrix of $ {H}^{H} {H} $ We all know that Discrete Fourier Transform (DFT) corresponds to circular (not linear) convolution. That is to say, if $x(n),h(n)$ and $y(n)$ is the original signal, the filter and output signal in domain and $X(\omega),H(\omega),Y(\omega)$ are their counterpart in Fourier domain, then if we want $Y(\omega) = H(\omega)X(\omega)$ we must have $y(n) = h(n) \circledast x(n)$, where $\circledast$ is circular convolution (modulo the length of x).

If $h$ is the linear operator that stands for circular convolution with $h(n)$, then I think $h^T$, the transpose of $h$, is a circular convolution plus a circular shift. My question is, is $h^Th$ also a circular convolution of some kernel?

The motivation for thinking about this is to look for a fast method to compute $(I+h^Th)^{-1}$ through FFT. I know it seems common to introduce FFT to compute this, but how exactly when h is a convolution kernel?

Thanks in advance,
Bob

If $ H $ is a matrix form of Circular Convolution then it is a Circulant Matrix.
Being a Circulant Matrix means it can be diagonalized by the Fourier Matrix $ {F} $:

$$ H = {F}^{H} D F $$

Where the matrix $ D $ id a Diagonal Matrix with the Fourier Coefficients of $ \mathcal{F} \left( h \right) $ on it main diagonal.
Also pay attention that we use the Conjugate Transpose $ {\left( \cdot \right)}^{H} $ operator as $ F $ and $ D $ are built by complex numbers (The space is over the complex numbers, hence the Adjoint Operator is the Complex Conjugate).

Now, it is easy to see that:

$$ {H}^{H} = {F}^{H} {D}^{H} F $$

Which is again an operator diagonlaized with the Fourier Matrix hence it is a Circulant Matrix -> Cicrcular Operation.
Basically it is a multiplication in Fourier Domain (DFT) by the conjugate of the DFT of the vector $ h $. By the DFT properties it means it is the flipped version of $ h $.

By the way:

$$ {H}^{H} H = {F}^{H} {D}^{H} F {F}^{H} D F = {F}^{H} {D}^{H} D F $$

Now, pay attention to $ {D}^{H} D $. It is a diagonal matrix multiplied by its conjugate. Namely it is, in the Fourier Domain, the Squared Magnitude.
In time domain it means it can be done by circular convolution with the correlation function of $ h $ (Or just do circular convolution in one direction and then in the other which is basically convolution with the flipped version).


How to Obtain Mathematically High Frequency and Low Frequency Component Separately Using Bilateral Filter? I have asked this question before in the sense that what does a high frequency and low frequency component signify in a image and i got satisfactory answers now i want to know that how i can get high frequency and low frequency component separately from the image.

I mean what changes i have to make in the equations of bilateral filters and its implementation in matlab so that i can get both high frequency and low frequency component of that image.

Because in my experiment i need both the high frequency and low frequency component of a image.

Bilateral filters takes a weighted sum of the pixels in a local neighborhood; the weights
depend on both the spatial distance and the intensity distance.

The value of a pixel assigned is given as

$$
      BF[I_p] = \frac{1}{W_p} \sum_{\substack{ q\in S }} G_{\sigma_s}(\parallel p-q \parallel) G_{\sigma_r}(\mid I_p - I_q\mid) I_q
      $$

Do I have to make changes to this equation of bilateral filter to get high frequency and low frequency components separately .

What component does this equation signify is it the high frequency component or low frequency component ?

One you have convolved your original image with the gaussian kernel (of appropriately chosen variance), then simply subtract the original image from the result of the convolution. This will give you the image with only the high frequency components of your original.


Seamless looping of a signal without pops I'm trying to seamlessly loop an audio signal. When I loop the current audio signal I have, I hear audible clicks or pops right as it loops back to the beginning. I want to be able to loop the audio signal without having any pops.
I have tried changing the signal near the end by writing an interpolation which follows this formula
$$
B_n=\frac{(A_n d_1) + (Cd_2)}{d_1+d_2}
$$
Where $B$ is the output, $A$ is the audio signal itself, and $C$ is the first signal value within the audio signal. $d_1$ equals to how many signal values until the first sample. $d_2$ equals how many signal values since the specified start. And I choose when the specified start is based on how many milliseconds I wish for the interpolation to be. In this case I gave 15ms a try. So, the specified start will be 15ms before the end of the signal.
And in simple terms; it's a crossfader that is linearly transitioning from the signal source $A$ onto the single value $C$.
The end result is this.

The transition points are highlighted via the cursor. Also, do not mind that both signals are not exactly the same. For the upper signal, I cut a portion of the start so that it would start at a non-0 value.
Since the crossfade is fairly quick for human ears (15ms), I did not manage to hear any disintegration of audio. But is this good practice? I am new to signal processing and cannot tell whether or not this will work for all types of signals I might have. Is there any other way this could be done in a more "proper" manner?

But is this good practice?

Yes, that's a perfectly reasonable approach.

Will this will work for all types of signals I might have ?

There are a LOT of signals out there and there is always an outlier. In your case, I think the most "vulnerable" would be low frequency sine waves. A 40 Hz sine wave has a period of 25 ms and looping one at an inconvenient point without audible artifact is tricky. So this all boils down to: how likely are these signal for your application, what level audibility is acceptable and what are the constraints around resources,  looping points & length, etc.

Is there any other way this could be done in a more "proper" manner?

If you have some leeway with the exact loop length and loop points, try to put the loop points at the nearest zero crossing. Ideally both zero crossing go in the same direction.
You can do an energy cross fade instead of a linear cross fade. This will minimize the volume dip in the middle of the cross fade, but at 15ms it's not going to make much of an audible difference.


What are the differences between fft and dct? what are the differences  between Fourier transform and cosine transform.....
i try to find the typical answer
the difference ,from where ,the size of input and output form,in time complexity,the domain ,the effect on image itself .

First, you should compare the DCT to the DFT, not the FFT. FFT is just a family of algorithms to approximate the DFT.

Also:

  [...]
  In particular, a DCT is a Fourier-related transform similar to the discrete Fourier transform (DFT), but using only real numbers. DCTs are equivalent to DFTs of roughly twice the length, operating on real data with even symmetry (since the Fourier transform of a real and even function is real and even), where in some variants the input and/or output data are shifted by half a sample. There are eight standard DCT variants, of which four are common.

Source: https://en.wikipedia.org/wiki/Discrete_cosine_transform


Detecting automobile honking I'm detecting vehicle honking from environmental sounds (engine noise, speech music, siren etc.) based on binary classification. Common audio features (spectral flux, centroid, mfcc, harmonicity etc.) seemed not doing well --- we get model with nice sensitivity but the problem is that the false-alarm rate is high. Failure cases includes sirens and speech.

Any features that describes most common vehicle honking sounds with discrimination from similar sounds for example, siren or speech? Thanks.

if you wanted to make use of the fact that the horn is normally two horns at different pitches, you might consider a pair of tone detectors (using either BPF or parametric EQ or a Goertzel filter and squaring and LPFing the squared output) similar to DTMF.  you have to detect a minimum portion of the whole audio energy in the two tones and they have to be reasonably equal energy for the two tones.  that will identify one model of horn.  you would need other pairs for other horn models.  but you should be able to detect a single horn honk that is louder than the other din.

detecting a cacophony of horns is a bigger problem.

you might want to gather a collection of sampled horn sounds (as Sagie suggested) and run a MATLAB spectrogram on each to see what you might have.


Difference between the average operation and the application of the equivalent LPF I would be glad if someone could answer my question. I have a bit stream of data and I want to extract the low frequency component. When post processing the data I can calculate the average value over a certain number of samples with the MATLAB operation mean(data(1:N))
Then I tired to build the equivalent LPF of this operation:
$$V_o(k) = \frac 1N \big[V_i(k-1) + V_i(k-2) + \ldots + V_i(k-N)\big]$$
Plotted the equivalent filter transfer function: freqz(ones(1,N),N,1e3).
And it looks like a sort of first order LPF. Now if I filter the data I don't get the same results of the average operation, filter(ones(1,N),N,data).
So what am I missing? The purpose of this investigation is to look for a better filter that will be implemented on FPGA to do real time processing.

You are missing $1/N$ factor in your filter if you want to do a moving average.
Moreover, FIR filter has a group delay of $\frac{N-1}{2}$. So if N is odd the value of mean(data(1:N)) should be at the $\frac{N-1}{2}$ sample of the filter output


Finding a signals complimentary sound / frequency I know about complimentary colors and how to get them, since colors are basically frequency / wavelengths.  Is it possible to find a signals complimentary, additive or subtractive frequencies in an audio file / from a signal.  I know how to use FFT to get the phase, amplitude, and frequency of a signal.  I'm just looking for a general idea of how I would go about creating similar color charts but for sound, we have a color charts of complimentary, additive or subtractive colors do we have one for the sound / audio range of humans?

Yes and no. Only monochromatic colors consist of a single wavelength, and those do not really "wrap around" in the form you describe. Color addition and subtraction have more to do with the three types of cone cells in our eyes than with wavelengths. Since there is no audio analogue to these cone cells, the concept of complimentary colors doesn't exist in the same way.

However, if you consider notes instead of frequencies, in the sense that two frequencies that are one or more octaves apart constitute the same note, then you can say that two notes are complimentary if they are a tritone apart.


Modeling a signal through the environment I have a signal that is composed of a sinusoid.  I will send it over the air through radios.  Part of the signal will go through processing on a computer, which means that it will be converted from analog to digital, and back again.

I have two questions.  How can I model any frequency shift that will occur in this system?  And, how can I model how the time of the signal may change?  That is, if the signal occurs in bursts (e.g. 200 ms of tone followed by 100 ms of absence of tone), how can I determine how the time may decrease?  Is the time changed only by the digital part of the system (due to sampling), or might the radio path also affect the duration?

Thank you

You can model frequency shift this way:

$s = s \circ exp(j \cdot 2\pi \cdot (f-\hat{f}) \cdot t)$,

where $\circ$ means Hadamard product and $s$ is input signal vector (e.g. frame or part of data stream of interest). You can implement this operation in the baseband or in carrier frequency. The effect will be generally the same.

As you've told your signal time can increas if multipath exists. You can model it as FIR filter of the order of your impulse response length. Refer to the channel modeling techniques for details.

Also your signal time can vary (increase and decrease) if relative motion of receiver and transmitter exists. This effect appears in the case of so called wideband model, e.g. sonic channel in the underwater communications. For such a case signal bandwidth is comparable to carrier frequency so Doppler shift is very different for low frequencies and high ones in the band of interest. It causes varying of signal duration in the receiver end. If you deal with radio signals your signal model is probably narrowband so this effect is neglectable.


Normalize data for transfer function determination? I have the input and the output of a system and I want to determine the transfer function of it. I am using Ident Tool from Matlab.
How is it more ok: to normalize the output and the input between 0 and 1 or to leave it as they are?

For what I know there is no reason to normalize data in general, if you use least squares like methods, all you need is that your input is persistently exciting, in the sense that the Hankel matrix $H$ of the observations must be such that $H^T H$ is nonsingular.


Parameters of standard audio channel filters I want to do some simulations of speech recorded on microphones. In order to simulate the channel response, are there any standard (FIR) filter parameters available that I can use?

Your question is not very specific, nor the answer will be. For example what do you mean by 'standard filter parameters'? I assume that you want to make given speech signals sound like ones being recorded with a specific microphone. Solutions would be

Hard way. You must obtain the impulse response of your microphone, for example by using Maximum Length Sequence or Sweep Sine method. In theory such measurement should be performed in anechoic environment (free space outside, well damped room, anechoic chamber), otherwise you will add even more reverberation. When you have your impulse response, simply convolve it with your signal and you will get very good approximation of how your microphone modifies the sound. Obviously speaker characteristics are also part of this impulse response, so you might want to use one with most linear amplitude response. In theory speaker response can be deconvolved to some extend by using reference microphone.
Easy way. Measure the frequency response of your microphone by using signals described above. You also measure the frequency response of the reference (maximally flat) microphone. Just subtract two amplitude characteristics and you get how your microphone distorts signals. Having that, use it as the input for the Parks–McClellan algorithm and design FIR filter that resembles this response in a best way.
Dummy way. Take the frequency characteristic provided by manufacturer in the manual, trust him, sample it, and use as the input for algorithm mentioned above. In MATLAB that would be firpm function.


FFT - preferred way of tone detection I am using FFT to detect certain frequencies (DTMF). I can clearly see the frequencies in the spectrum. Is there a preferred / optimal way to decide if a certain frequency is present or not by looking at the spectrum?

I can think at least of:

Bin absolute value (ie. bin amplitude)
Difference between maximum bin and second highest
Ratio between maximum bin and second highest
Ratio between maximum bin and sum of all others

Number 4 : Comparing narrow-band energy to the energy of the total spectrum is a commonly described method of DTMF detection.  This may involve using more than one FFT result bin and/or interpolation between bins.

The other methods depend on whether one has an accurate estimate of the signal level, noise floor, and how well the tone frequency is centered in a single FFT bin (related to frequency, sample rate, FFT length and window).


Calculate IFFT using only FFT? Say you start out with $A=\left\{3,−1,0,0\right\}$, run fft() and end up with $B=\left\{2,3+j,4,3-j\right\}$.

What specific steps would need to be done to obtain $A=\left\{3,−1,0,0\right\}$ again using some pre-calculation on $B$, then fft(), then maybe post-calculation steps before obtaining $A$? I have heard that this is possible but have not been able to work out how.

We want to do inverse DFT:
$$
x(k) = \frac{1}{\sqrt{N}}\sum_n y(n) e^{j 2\pi kn/N},
$$
but can only do forward DFT. So let's try to make it "look" like a forward DFT. Note that we want a minus sign in the complex exponential. We can introduce it via complex conjugation:
$$
\overline{x(k)} = \frac{1}{\sqrt{N}}\sum_n \overline{y(n)} e^{-j 2\pi kn/N},
$$
That is starting to look like a forward DFT now! So the recipe is:

Complex conjugate the given sequence that we want to inverse DFT
Calculate its forward DFT
Calculate complex conjugate of the result.

That gives you the inverse DFT of the original sequence.

(Note that I cheated a little bit on the scaling factor to make it easier!)


Why is it assumed that $x[n]$ is limited from $0$ to $N-1$ while evaluating DFT? I am a total beginner in this topic of DFT. I get that the series must be finite for DFT calculation. But everywhere we are assuming that this series must be limited from $0$ to $N-1$. How to evaluate DFT for a sequence extending from $n=2$ to say $n=N+1$?

I am following this video

https://www.youtube.com/watch?v=GDFTb-BwA0o

Here in the minute 21:20 to 21:50, Professor seems to be really strict about this particular range. Why is this so?

How to evaluate DFT for a sequence extending from n=2 to say n=N+1?

Formularistic: trivial variable substitution of the running variable $n$ by $\tilde n = n+2$, done.

Looking at it as a signal: When doing an DFT, it helps to think of it simply as mapping of $\mathbb C^N\mapsto\mathbb C^N$, i.e. a complex vector goes in, same size complex vector comes out.

Your "shifted" input hence is simply a different input value for the DFT.

Now, when we interpret the DFT as a mapping between discrete time- and discrete frequency-domain, we have to assume that the time signal is periodic with the DFT length. That means the two "new" values are just the oldest two values. I.e., your reindexing just leads to a circular shift of the values in the vector.


Is the system $y(n)=x(n^2)$ time invariant? According to my solution,

$$y(n) = x(n^2) \tag{1}$$

for $y(n,k)$ i.e output for delayed input $x(n-k)$

$$x_2(n) = x(n-k)$$

so $$y_2(n) = x_2(n) = x\big((n-k)^2\big)$$

and for delayed output signal $y_1(n)$, replace $n$ by $n-k$ in equation (1), so we get,

$$y_1(n) = x\big((n-k)^2\big)$$ and therefore system is time invariant

But in the answers to the book in which this question it says the system is time variant. Can anyone point out the mistake in my steps, and give correct way of solving this sum?

The system is obviously time-varying. However you make the following usual mistake in the first case of showing the response to the delayed input $x_2[n] = x[n-k]$ as follows.

Writing the output $y[n]$ in the formal manner: $$ y[n] = \mathcal{T}\{x[n]\} = x[n^2] $$ and denoting the shifted input as $$x_2[n] = x[n-k]$$ you would have the corresponding output $y_k[n]$ as:

$$y_k[n] = \mathcal{T}\{ x_2[n]\} = x_2[n^2] = x[n^2 -k] $$ which is defintely not equal to $y[n-k] = x[(n-k)^2]$ hence the system is time-varying.

You should have defined the output as $$ x_2[n] = x[n-k] \implies x_2[n^2] = x[n^2 -k] $$ but you didn't. See that the equality sign between $x[n]$ and $x_2[n]$ enables you to operate on $n$ only, not the whole argument of both sides.


Hydrophone kit to transmit (and receive) underwater May I ask you if there exists a complete kit (e.g., by including the amplifier, etc.) of professional hydrophones to transmit and receive signals in underwater environments? If not, can I ask you for help about the configuration and which are the best instruments that can I buy? Scope: Research and Development.

Yes, may be you are confused between Hydrophone and Modem.Hydrophone is just a passive device that records sound underwater. And what you just describe in your question is modem, which emit and receive sound waves underwater. there are definitely so many options you have. There are number of factors you need to consider before buying, mostly the operating frequency as most of the equipment in underwater works with KHz frequency. There are some amazing Software Defined Modems are available in the market but definitely those are expensive. If you want to just design and develop some communication protocols then you can use the UnetStack, basically convert your laptop sound card into a modem and have pretty good functionality. If you’re planning to buy some cheaper uw-modem, you can try buying popoto development kit, which includes transmitter and receiver. If you’re interested in some of the hydrophones then it worth give it a shot to Oceansonic icListen SC2, i must say this is an impressive military grade performance Hydrophone. This guys also had icTalker, but this can not generate custom signals, but definitely generate some chirp, continues sweeps, tones etc.
Also there are good market players in underwater domain such as,

Subnero
Benthos
EvoLogics GmbH
WHOI Micromodem
Popoto


Change of bin magnitude when phase of signal changes when doing FFT? I am working on an implementation in which I need to apply an FFT to determine which of 4 frequencies are present in a signal.

I am doing a 64 point FFT on a buffer which only is partly filled (the remainder of the buffer is zero-ed).

When looking at the screenshots of the signal in the time domain, you will notice that the signal is roughly the same in each screenshot, but it has different phase in each screenshot.

When looking at the screenshots of the FFT of the two versions of the signal, you will notice that the magnitude for bin 4, which corresponds to the signal frequency, is 874.83 in one of the screenshots but 366.94 in the other, even though the only difference between the two versions of the signal is a phase difference (?)

I don't understands why this happens and how to deal with this. I am generating the signal myself and want to make sure that phase doesn't play a role here. When the signals are more or less the same in amplitude, I would expect their magnitude after calculating the FFT to stay the same.

In case you need this information, the sampling frequency is 7800 Hz and the actual frequency involved is 487.5 Hz which is exactly 7800/64*4, or the start of bin 4.

In the screenshots of the FFT, the magnitude data is plotted upside down (i.e. zero is the top of the screenshot).

Can anyone please shed some light on this? Thanks for helping!
B

The FFT magnitude, for non-zero sinusoidal-like signals, will only be the same, given a change in phase relative to the window, for signals that are exactly periodic in the window width (the non-zero padded portion of the FFT aperture).


Measuring vibrato rate of a sung note I'm doing a project where I need to measure the rate of vibrato cycles in a sung note under various conditions. The stimulus will be a several-second-long audio sample of a held a capella note with vibrato; the output should be the rate of the vibrato repetition (e.g. 4Hz, 10.5Hz, etc). The audio will be captured specifically for this study; I won't be sampling existing audio, and don't need to deal with interfering signals. The underlying pitch of the note can be whatever is needed, although precision probably won't be possible.

Here's an example sonogram:

Ideally, this would be done with a Linux command line chain or shell script, but if I need to write some Python that would be fine, too.

This publication may be of interest for you:

Jonathan Driedger, Stefan Balke, Sebastian Ewert, and Meinard Müller: Template-Based Vibrato Analysis of Music Signals, in: Proceedings of the International Conference on Music Information Retrieval, New York, USA, 2016, pp. 239–245. [PDF]

Furthermore, there is an accompanying website for this paper.


Applying Viterbi algorithm to compensate for ISI in PSK31 I am implementing a receiver for the amateur radio mode PSK31. It is effectively a BPSK modulation with a pulse shape of a raised cosine twice the width of one symbol. Here are four consecutive pulses, with 4 samples per symbol, as they would be transmitted:

The trouble is matched filtering introduces ISI. 4 pulses again, after a matched filter in the receiver:

If my understanding is correct, the issue here is one pulse is non-zero at the ideal decision point of the adjacent pulses.

However, my thinking is this is predictable, and the interference never extends beyond the adjacent symbols. So if the previous and next symbols are known the ISI can be cancelled by subtraction of the interfering adjacent pulses.

The trouble is the previous symbol isn't known with certainty. I think the Viterbi algorithm is a possible solution. And I've watched some videos on the algorithm and think I understand how the algorithm works, given a particular trellis.

What I don't understand is how to apply the Viterbi algorithm to this particular problem. The material I've found is either theoretical with math that blows my head off, or some other application of the algorithm. Could someone please provide an intuitive, practical explanation in this context?

(After some deliberation, I have a guess at answering my own question. Someone tell me if I'm wrong.)

The first step is to look at the contribution of each pulse at the decision point. The filter is normalized so a single pulse results in 1. And in this particular case, the tails of adjacent pulses work out to 1/6.

We can then look at any possible combination of three pulses, summing their contributions together. If we say the three pulses are $x_{n-1},\ x_{n},\ x_{n+1}$, then the output $y_n$ at the decision point is:

$$ y_n = 1/6\:x_{n-1} + x_{n} + 1/6\:x_{n+1}$$

To give some examples:

111: 1/6 + 1 + 1/6 = 1.33
010: -1/6 + 1 + -1/6 = 0.67
110: 1/6 + 1 + -1/6 = 1

Say state 01 means the previous bit was a 0, and the one before a 1. We can then ask, "if we add another bit, what does the output become?" For example, if the next bit is a zero that yields 001, with an output of -1 according to the formulation above.

We can then build a state table:

╔═══════╦═══════╦════════╗
║ input ║ state ║ output ║
╠═══════╬═══════╬════════╣
║     0 ║    00 ║ -1.33  ║
║     1 ║    00 ║ -1     ║
║     0 ║    01 ║ -1     ║
║     1 ║    01 ║ -0.67  ║
║     0 ║    10 ║ 0.67   ║
║     1 ║    10 ║ 1      ║
║     0 ║    11 ║ 1      ║
║     1 ║    11 ║ 1.33   ║
╚═══════╩═══════╩════════╝

And the associated trellis diagram:

The notation "0/-1.33" means "to send a 0 from this state, the output is -1.33".

Either the magnitude or the square of the difference between the observed values and the expected values can be used as the metric. The Viterbi algorithm will find the path through the trellis with the lowest metric regardless of what metric is used, as long as it's something that monotonically increases with the distance between observation and expected.

The Viterbi algorithm can then be followed in the usual way to determine the most likely sequence that explains the received signal, knowing how adjacent pulses will influence the received values.

I've found these videos especially helpful in understanding the evaluation of the Viterbi algorithm:

Viterbi Algorithm, by AppBooke Xperience
Viterbi Algorithm, by Keith Chugg


variance of the product of two samples with awgn problem solved itself, sorry for your inconvenience. I'll try to post better questions next time.

There is a paper, which is more general than Gaussian distributions.

  Bohrnstedt, George W., and Arthur S. Goldberger. "On the exact
  covariance of products of random variables." Journal of the American
  Statistical Association 64.328 (1969): 1439-1442.

https://www.jstor.org/stable/pdf/2286081.pdf?casa_token=m1_rSG4hiwIAAAAA:sGjisD6dekanCat5EBGl0yvdn1EJ5LY783fq4VKdxNhd0RA4z4izkcyDJ6yU7MylDke_qD335SWUNu1gJ5opTzXFWIVQWDrjoMvxcFgvMBEgG6r8En0W

Their equation (9)
$$
 V(xy) = E^2(x)V(y) + E^2(y)V(x) + V(x)V(y). \quad (9)
$$
reduces to your result


PSD looks to clean/continuous I generated an signal with a base frequency of $50\textrm{ Hz}$ and its harmonics and quantization noise in it (as an array in MATLAB). It's sampled fast enough to make sure there is no aliasing etc.

I calculate the PSD with a rectangular window in MATLAB. Due to the fact that the signal is synthesized in MATLAB I know where my harmonics are, so I chose length of my FFT to have a frequency resolution of $5\textrm{ Hz}$. I did this to be able to see my hamonics without the spectral leakage effect.

My problem is, that the PSD estimates differ from my expectations. I expected to see spikes caused by my harmonics and a kind of random noise floor. But what I get is kind of a continuous looking PSD with spikes as you can see in the image below:

The $x$-axis is the frequency in $\textrm{Hz}$ and the $y$-axis is in $\textrm{dBc/Hz}$ normalised to my base frequency power.

My Matlab code is:

NFFT = length(bitstream)-1; %df = fs / NFFT = 5 Hz
fVals = fs*(0:NFFT/2)/NFFT;

% calc fft
X_f = abs(fft(bitstream,NFFT));
% calculate periodogram (one sided)
pxx = ((X_f(1:NFFT/2+1))/NFFT).^2;
% plot
plot(fVals,10*log10(pxx/max(pxx)),'black');

The MATLAB function periodogram() as shown below (with rectangular window and same FFT length) produces the same output, so I thought my code is correct.

pxx = periodogram(bitstream,rectwin(length(bitstream)),NFFT);

Can you tell my why I see that?
Is the periodogram method not suitable for that case?
Is it caused by the rectangular window and/or the frequency resolution? If I use the built-in non-parametric estimation functions like pwelch(), periodogram(), pmtm() etc. without any further input it looks more noisy (as expected).
Or is it even better than expected because one of the frequencies in my simulation is the greatest common divider of the sampling frequency (using e.g. a signal with a different frequency like $49\textrm{ Hz}$ also leads to more noisy results)?

I believe this is due to your choice of using a rectangular window, which has the best frequency resolution (1 bin) but worst spectral leakage (dynamic range). The spectral leakage is of the form of a $\textrm{sinc}$ function (with nulls at $1/T$ where $T$ is your window length in time), which has an envelope that rolls off very slowly ($1/f$) in frequency. What I believe you see above is the cumulative effect of the $\textrm{sinc}$ function in frequency due to your window convolving with each tone.

My favorite window choice is the Kaiser window kaiser(N,beta), as you can control the precision and dynamic range trade directly through the use of beta. Beta factors I typically use range from 4 to 12 but you can experiment with it and see what works best for you.

Try this to see the window performance (MATLAB / Octave) in the frequency domain (showing the main lobe width (precision) and sidelobe level (dynamic range).

wind = kaiser(N, beta)
freqz(wind)


What is databending and how does it actually work? http://www.google.com/search?q=databending

It is, from what I understand, opening images as audio, applying effects, exporting back as image and voila! glitch-art

Example:

credit

But how does it really work? How/what change in frequencies(?) affect the pixels? It can't be totally random since some of them aren't completely distorted..

How does time(in audio) relate to the dimensions and colors of the image?

It's pretty simple, really.

Basically, this relies on the fact that the images must be saved in an uncompressed file.

If you think about how an image is stored, it's functionally a bunch of sequential numbers:

(PX 0,0 - R:x, G:y, B:z)
(PX 1,0 - R:x, G:y, B:z)
(PX 2,0 - R:x, G:y, B:z)
(PX 3,0 - R:x, G:y, B:z)
....
(PX 639,0 - R:x, G:y, B:z)
(PX 0,1 - R:x, G:y, B:z)
(PX 1,1 - R:x, G:y, B:z)
....
(PX 639,479 - R:x, G:y, B:z)

For most bitmaps, R, G, B are either 8 or 16 bit values (1 or 2 bytes). Therefore, the "image" is really just a very long sequence of numbers, with some agreed-upon way of mapping those numbers into pixel-positions and colors.

Now, how is an uncompressed audio file represented? Very, very similarly:

(Sample 0 - R:x, L:y)
(Sample 1 - R:x, L:y)
(Sample 2 - R:x, L:y)
(Sample 3 - R:x, L:y)
....
(Sample z-1 - R:x, L:y)
(Sample z - R:x, L:y)

Most audio samples are 16 bit (some are 24, it can vary) per sample. However, again, the audio file is really just a whole bunch of numbers, with a agreed-upon method for translating the sequences of numbers back into an audio stream.

Now that you understand how the files are structured, you can start to see how the audio effects relate to the images. First, since uncompressed images are generally stored in a raster-scanned topology, audio effects are primarily going to have a horizontal "smearing" effect.

Additionally, you can see that basically the effects are going to permute the images based on the effective horizontal frequency components of each horizontal row of pixels.

The "Random EQ" setting boosts some the frequency components of the horizontal edges in the image, and the rippling that surrounds each edge is a function of the impulse-response of the filter.

The "Distortion" image looks to be adding some high-frequency noise, in addition to some apparent posterization, which makes me think it may decimate the internal sample-values in some manner.

The "Bass-Boot" is clearly exaggerating the swings between wide areas of low-values (e.g. black) and high-values (e.g. light colours).

I'm not sure what exactly the phazer is doing (I don't know what it would sound like, so I can't extrapolate to the DSP it involves).


What's the Q of a pole at the origin of the s-plane? I'm writing a simple function in Python to return the Q of a circuit, given its poles.

The Q of a system with pole p, can be found using $Q=-\frac{\left|p\right|}{2*Re(p)}$.  Obviously, in a natural system, such a pole must come as part of a pair, and both would yield the same result.

For coincident poles on the real axis, this gives a Q of 1/2 (critically damped).  For poles on the imaginary axis, that gives a Q of $\infty$ (undamped).

Writing the function left me with the question, as a matter of completeness, of what to do if the pole is at the origin.  Is there a customary value to assign to poles at the origin?

Background material:

For anyone interested in where the equation relating pole location to Q comes from,   It derives from the 2nd order transfer function:

(source)

Graphically described here:

(source)

A decent discussion of Q and damping in 2nd order systems can also be found here:

(source)

which gives an excellent review of how the pole placement corresponds to various damping characteristics.

Note that the formula given in your question is valid for a system with a complex-conjugate pole pair $p$ and $p^*$ with $\text{Re}(p)\neq 0$. As you've correctly pointed out, if for $|p|>0$ the real part $\text{Re}(p)$ approaches zero, the $Q$ factor approaches infinity. This is the case for a pole pair on the $j\omega$-axis with $|p|>0$, i.e. not at $s=0$.

Also note that you can't use that formula for two distinct real-valued poles. That situation corresponds to an overdamped system with $Q<\frac12$.

The problem with a double pole at $s=0$ is that there is no reasonable way to evaluate the formula given in your question. The corresponding limit does not exist.

Having said that, we can think of a reasonable definition of $Q$ for the case of a double pole at $s=0$.
For $p=p^*=0$ the transfer function of the corresponding system is (assuming no finite zeros)

$$H(s)=\frac{A}{s^2}\tag{1}$$

with some constant $A$. The impulse response is

$$h(t)=\mathcal{L}^{-1}\{H(s)\}=A\,t\,u(t)\tag{2}$$

We have to ask ourselves what we mean by $Q$ in this case. Since the damping is zero, you could say that $Q=\infty$. I haven't come across a definition of $Q$ in the case of a double real pole at $s=0$, and I think that for this degenerate case there cannot be any useful definition. If I had to decide on a value, I'd say $Q=\infty$, in analogy with all other pole pairs on the $j\omega$-axis.


Direction of arrival and distance to source Suppose we have $3$ receivers $a,b,c$ and one wave source which are all located on the same $XY$ plane and we do not know the position of only $c$ on the plane.

Assume a plane wave is transmitted $s(t)$ from an aribitrary point source in far field.

Using the received signal $s_a(t)$ and $s_b(t)$ by pattern matching would it be possible to infer the direction of arrival at $a$ and $b$ and time delay of signal to reach $a$ and $b$?
Suppose the received signal at $c$ is $s_c(t)$ would it be possible to guess the distance and direction of source to $c$ and also $(x,y)$ coordinates of $c$ with $s_a(t),s_b(t)$ and $s_c(t)$? If not what other minimal information do we need to find these?
Is there a reference for these kind of problems for multiple sources and $XYZ$ coordinates?

Your statements are a little confusing. First you say that the sensors and the source all are on the same XY plane and then you mention that the source signal $s(t)$ is transmitted from a point source in far field. I will go with the latter.

In array signal processing, the source's or the sensor's position is usually estimated using the concept of relative time delay. For a given source at location $\boldsymbol{\beta} = \left[\phi, \theta \right]$ where $\phi$ is the azimuth angle and $\theta$ is the elevation angle, the signal received by the $k^{\text{th}}$ sensor at position $\boldsymbol{p}^{\text{T}} = \left[x_k,y_k,z_k=0 \right]$ will experience a time delay of

$$ \tau_k(\boldsymbol{\beta}) = \frac{1}{c}\boldsymbol{p}^{\text{T}} \begin{bmatrix}{\cos\phi \sin\theta \\ \sin\phi \sin\theta \\ \cos\theta} \end{bmatrix} $$

To simplify this setting, let us only work with the elevation angles. Let the sensor $a$ is at location $\left[x_a, y_a \right]$ and the sensor $b$ is at location $\left[x_b, y_b \right]$. So the distance between them will be $d = \sqrt{\left(x_a-x_b\right)^2 + \left(y_a-y_b\right)^2}$. Let the plane wave reaches sensor $a$ first. Solving a basic geometry will tell you that the extra distance the plane wave has to travel to reach $b$ is $\delta = d\tan\theta$. So the time delay will be

$$\tau = \frac{\delta}{c} = \frac{d\tan\theta}{c} $$

So to answer your first question, Yes! It would be possible to estimate the direction of arrival at sensor $a$ and/or $b$. You do not need to get it for both the sensors as they are receiving the signal from the same source. You would need to know the time delay which you can get by cross-correlating the signals received by the two sensors. Once you know the time delay, you can use it to get the direction of source i.e., $\theta$.
To answer your next question, Yes! Once you find out the position of the source. Find out the time delay between the sensors $a$ and $c$ and the sensors $b$ and $c$ and then apply the same equation. You should pay attention to the sign of the delay i.e., between sensors $a$ and $c$ whether is reaches $a$ first or $c$ first and same with $b$ and $c$. Since the formulation only has the distance $d$ and not the specific coordinates, the absolute value of the delay will be the same in the circle of radius $d$ with origin being the sensor $c$.
I have included the formulation of the delay in terms of the 3-D coordinates. There are a lot of good resources available online to help you understand this in detail. A basic google search will go a long way.

There are a lot of other ways to solve the direction of arrival estimation problem. You should find out more and select the one based on your needs.

Hope this helps!


what is nyquist rate of $h(t)\cdot h(t)$ and $h(t)\circledast h(t)$ Let's say we have $h_c(t)$ as a continuous-time signal with bandwidth $B$ and we would like to sample it.
To be able to reconstruct it correctly, the sampling rate must be greater than $2B$.
Now assume we have $h_{c1}(t)=h_c(t)\cdot h_c(t)$ and $h_{c2}(t)=h_c(t)\circledast h_c(t)$ what is Nyquist rate of these two signals? Is it $4B$ for both of them because each $h_{c}(t)$ needs $2B$?

squaring a signal effectively doubles the bandwidth of the signal (and doubles the Nyquist rate, which is the lower limit of the sampling rate).  convolving a signal with anything, even itself, does not in general increase the bandwidth.


Evaluating discrete spectral density at only a few frequencies I'm trying to obtain the spectral density at three particular frequencies for a computational chemistry problem that I'm working on (if you are curious, it has to do with the estimation of Nuclear Overhauser Effects from molecular simulations). Three is a small enough number, especially given the length of the signals being sampled, that it seemed helpful to use a pruned DFT to do this (a Goertzel-like algorithm). Now I'm wondering if I can get away with that, or am somehow thinking about the problem incorrectly. I'm not a signal processing expert by any stretch of the imagination, so I am worried about missing something obvious.

In the chemistry literature, these effects are always functions of the spectral density of the autocorrelation function of some order parameter (in this case, if it matters, the dipole interaction tensor between an internuclear vector and the laboratory magnetic field, which I'm just going to be treating as 'the signal' below). From wikipedia's article on autocorrelation, for a signal $X(t)$:
$$
\begin{align}
F_R(f) &= \mathrm{DFT}[X(t)] \\
S(f) &= F_R(f)F_R^*(f) \\
R(\tau) &= \mathrm{IDFT}[S(f)]
\end{align}
$$

Since the spectral densities are given as the $\mathrm{DFT}[R(\tau)]$ in the textbooks/papers I've found discussing this physical effect, and since I only need the spectral density at three frequencies (0, 600, and 1200 MHz, if that matters) I thought the most straightforward solution was to use a Goertzel-like algorithm (this especially because much computation and disk IO not related to the DFT needs to be done to obtain each sample, so the single pass character of Goertzel-like algorithms is nice for this application).

$$
\begin{align}
S(f) &= \mathrm{DFT}[R(\tau)] \\
F_R(f_0) &= \mathrm{Goertzel}_{f_0}[X(t)] \\
S(f_0) &= F_R(f_0)F_R^*(f_0)
\end{align}
$$

To reiterate the question, is the above reasoning valid? I've gotten funny results that have been hard to troubleshoot so I'm going back to questioning my underlying assumptions. It's also quite possible that all this works, but I'm applying it incorrectly, so if there are restrictions on which frequencies can be monitored in this way (or other things that would be obvious to someone in the field) that'd be helpful to hear.

The algorithm I'm using is sourced from Clay Turner's 'Oscillator Theory' approach to single frequency DFTs: here's C-like sample code he provided for it on a comp.dsp post from a while back:

// The input data is in x[], the data has N samples. And the bin number is w.

y1=0;
y2=0;
k=2*sin(pi*w/N);             // not 2 pi !!

for (j=0;j<N;j++) {
      y2=y2-k*y1+x[j];
      y1=y1+k*y2;
    }

// And the energy is simply

E = y1*y1 + y2*y2 - k*y1*y2;

Note that I chose this after doing some looking when I read Gentleman's 1969 paper suggesting that Goertzel's algorithm performs poorly for low frequencies (I need zero, so that seemed bad). I am working in C++, but because my 'signal' is matrix-valued my code does not look exactly like this. (Yes I can post it, but it might be TMI so I won't unless someone feels that might help).

I believe the OP's question is simplified to the following (confirming I didn't actually miss the salient question):

  Given we can compute a power spectral density from the DFT using the
  conjugate product as:

  $$S(k) = |X(k)|^2 = X(k)X^*(k)$$

  Can we use a more efficient algorithm, such as the Goertzel, to
  compute a subset of $S(f)$ when only a few points are needed. The OP
  also wanted to understand if there was a limitation with using the
  Goertzel specifically, especially when only the lowest frequencies are of interest.

The answer is yes we can use a more efficient algorithm and specifically we can approximate DFT bins directly using the Goertzel. A windowing function is also recommended in either case (full DFT) or with the Goertzel to minimize aliasing issues (spectral leakage between bins).

The Goertzel is more efficient than the FFT when $M<log_2(N)$ where N is the total number of samples and $M$ is the total number of bins to compute. Other advantages are we can center each bin of interest at exact frequency values.

The Goertzel algorithm does have an error bound proportional to $N^2$ that is most pronounced near $\omega = 0$ and $\omega= 2\pi$. An alternate approach that would be identical to FFT results is to compute the DFT bins directly using the DFT formula, with the windowing suggested above also shown:

With windowing applicable to any approach: $x_{win}(n) = w(n)x(n)$, where $w(n)$ is a window, such as Kaiser, Hamming, Blackman, etc.

$$X(k) = \sum_{n=0}^{N-1}x_{win}(n)e^{-j2\pi nk/N} \tag{1}\label{1}$$

Which assuming a complex $x_w(n)$ requires a total of $4N$ real multiplies and $4N$ real adds for each bin computed.  In comparison, the FFT would compute all $N$ bins with $2Nlog_2(N)$ real multiplies and $2Nlog_2(N)$ additions.

Therefore the total number of bins where it there is savings in computing directly using the DFT is given by $M$ for:

$$4NM < 2Nlog_2(N)$$

Resulting in:

$$M < \frac{log2(N)}{2}$$

So for example, if we had a 1024 point DFT, if less than $log_2(1024)/2 = 5$ bins were needed, it would be more efficient to compute them with the DFT equation directly rather than use the FFT. The Goertzel extends this threshold to 10 bins.

Further, as we can in the Goertzel, we can modify the DFT equation to center the bin on any frequency $\omega_o$ for $\omega_o \in [0, 2\pi)$ as follows:

$$X(\omega_o) = \sum_{n=0}^{N-1}x_w(n)e^{-j\omega_o n} \tag{2} \label{2}$$

In the DFT equation, $\omega$ is limited to discrete frequencies given by $k\omega_o$ where $\omega_o = 2\pi n/N$, thus it should be clear how $\ref{2}$ and $\ref{1}$ are related and how the form in $\ref{2}$ allows us to place the bin on any frequency along the continuous $\omega$ axis rather than the discrete axis for the DFT.

If windowing, please also refer to the last paragraph in this post to properly compensate noise and signal terms after using a window:

Blackman-Tukey Autopower equation

The FFTW page also provides a useful reference on Pruned FFTs and specifically an approach when the interest is in only the first K outputs which may be applicable to the OP's challenge:

http://www.fftw.org/pruned.html


Dynamic Time Warp through SPTK tool I'm an IT student and got an assignment to do about Dynamic Time Warping(DTW) using the Speech Signal Processing Toolkit (SPTK) and comparing some words spoken by 2 speakers and finding the similarities. I managed to get the SPTK working, collected 8 people(4 females, 4 males) who recorded 8 words each for me(same words for every person) and saved them as files with a .wav extension.

My .wav files are in the following format:

RIFF (little-endian) data ;
WAVE audio ;
mono ;
sampled at $16\text{ kHz}$.

I transfered every .wav file into .raw data files with the SPTK conversion function :

wav2raw +f source_maleA.wav

The .raw file I get is:

Unsigned 8 bit ;
Sampled at $8 \text{ kHz}$ ;
Mono.

I transfered every .raw file to a .mcep file with this line of code:

x2x +Sf < source_maleA.raw | frame -l 320 -p 160 -n | window -l 320 -L 512 -n 1 | mcep -l 512 -m 13 -a 0.42 > source_maleA.mcep

After that, I went to compare the .mcep files with this line of code:

dtw -m 13 -s Scorefile target_maleB.mcep < source_maleA.mcep > source_maleA_target_maleB.dtw

The output of that command line (the "Scorefile") is a binary file.

I programmed a code that decodes the binary sequence written inside of the Scorefile and the output is some kind of decimal number: $0.678382$, which could imply a percentage of something.

Since by the SPTK documentation they are using some .short files of their own, I tried to remove the header of the .wav files and use them as .raw format and then convert it to .mcep file for the further use.

Since my understanding of DTW is limited by itself, plus my understanding of the SPTK tool is also limited, could anyone tell me if what I did is good?

I divided the files into frames of $320$ samples and periods of 160 and I used $x[0]$ as the first point of the first frame. I applied a Blackman's window of 512 samples length with the normalization computed by this function:

$$\sum_{n=0}^{L-1} w^2 (n) = 1 $$

After that I used the mel cepstral analysis to get the coefficients with $alpha = 0.42$ and I set the mc order to $13$.

Since when changing the length of the window, frames, period etc. etc. I don't get any better result, could anyone explain me what would be a good setup for all those factors?

Link to the SPTK (for the documentation and examples)

The documentation clearly says (emphasis mine):

  If –s option is specified, the score calculated by
  dynamic time warping, that is, the distance between the test data and the reference data is output and sent to Scorefile.

That means that your decimal number, $0.678382$, is the actual DTW distance between your 2 speech samples.

Now DTW is a distance measurement, and as any distance it only gives you relative information. For example, you could set a treshold and says that under this treshold, the 2 signals correspond to the same word.  If you apply DTW on a pitch contour, you could determine whether the 2 melodies are similar (that's used in query by humming).

Of course, the physical meaning of the distance depends completely on the one of the time series compared : if you compare voltage values, or frequency values, etc.
In your case, you compare normalized volumes mapped to unsigned 8-bit data.

I advise you to learn more about DTW and to try to analyse what the distance actually tells you, and how that relates to your application


Discrete Fourier Series of square of a signal I am having trouble in deriving a result, which asks me to find a relation between $X_{1}(k)$ and $X_{2}(k)$ where $X_{1}(k)$ are the DFS coefficients of single period of signal $x(n)$ and $X_{2}(k)$ are the DFS coefficients of two periods of the same $x(n)$. My approach is this.

$$X_{1}(k) = \sum_{n=0}^{N-1}x(n)W_{N}^{-kn}$$
and $$X_{2}(k) = \sum_{n=0}^{2N-1}x(n)W_{2N}^{-kn}$$
Since we know that $W_{N} = W_{2N}^{2}$. That means that, if I square $X_{2}(k)$, then I get this
$$X_{2}^{2}(k)=\sum_{n=0}^{2N-1}x^{2}(n)W_{2N}^{-2kn}=\sum_{n=0}^{2N-1}x^{2}(n)W_{N}^{-kn}$$
Now this is where I can't move further. What should I do with the last expression. I mean what are the DFS coefficients of the square of a sequence? And also note that now we are over two periods ($2N-1$).

Here is an outline (in case it's homework)

Use a different frequency index for $X_2(l)$
Split the $X_2(l)$ into two sums: from $0:N-1$ and from $N-1:2N$
Use the periodicity of x[n], i.e.  $x(n) = x(n+N)$
Substitute $m = n-N$ in the second sum
Recombine the two sums again
Stare what happens to even $l = 2k$ and odd $l = 2k+1$ terms