best window for close frequency components I'm currently designing a receiver which has to determine whether a signal contains specific frequencies or not. The frequencies are at constant 215Hz difference:

{16,351  16,566  16,781  16,996  17,211   17,426  17,641  17,856  18,071}

I sample at 44100kHz and I have 1024 samples each time. Therefore my resolution is 44100/1024=43.06Hz, and the bins centers are multiples of 43.06:

Bin[1]=43.06Hz,... BIN[385]=16578.10    BIN[386]= 16621.16
BIN[387]=16664.22    BIN[388]=16707.28     IN[389]=16750.34
BIN[390]=16793.40    BIN[391]=16836.46    BIN[392]=16879.52
BIN[393]=16922.58    BIN[394]=16965.64    BIN[395]=17008.70
BIN[396]=17051.76    BIN[397]=17094.82,....

My algorithm is based on comparing the values of the desired bins (or the nearest ones) to the values of the bins in between.
for example, if I wish to know if 16,566 exists or not, I compare BIN[385] to BIN[386],BIN[387],BIN[384],BIN[383] and if it's value is significantly larger (say by a hundred) then the average value of the side ones then I conclude that this frequency does exist.

My question is as follows: considering the fact that I work in high noise environment, and the algorithm I use is the one above, what would be the best window for me to apply before the FFT, other than the obvious rectangular window I currently use.

Thanks

I'm presuming that you have some or all 9 of the input frequencies available at one time, and you just want to find out which of them is present (you didn't specify in your problem statement).

Two methods come to mind: 1) compute DFTs at the frequencies of interest, or 2) use the Gaussian ratio technique, as described in reference_1 and reference _2.

Regarding 1) above: When Moses came down from the mountain, he did not carry tablets of stone on which it was commanded that in the equation: F = k*(df) = k*sample_rate/N, that k has to be an integer.  In other words, you can compute a DFT bin centered at ANY frequency, not just those that correspond to the usual k = integer of your FFT.  For example, given your sample rate and N, a frequency of F = 16351 would correspond to a DFT bin computed at k = (N*F)/sample_rate = 1024*16351/44100 = 379.669 (ie: between bins 379 and 380 of your FFT).

So for the DFT, X(F) = (sum over n){x[n]*e[-j*twopi*n*k/N]}, where X[F] is the complex result, N = 1024, F = 16351, k = 379.669, and n = 0 to N-1.

Some of the other values for your problem are:

Frequency   and value of   k

16,566  ==> 384.662

16,781  ==> 389.659

16,996  ==> 394.646

I'll leave it to you to figure out the rest.  And, as you can see, you've got about 5 bins of separation between the input frequencies.  That may, or may not be enough based on your SNR and method used.  And you might want to use N = 1025 instead – it makes your input frequencies much more 'bin-centered.'

Method 2) above uses an FFT, which, given the right circumstances, can be extraordinarily accurate when computing frequencies.  Note that in the first reference above, by using a Gaussian window, and then taking the log ratio of amplitudes from adjacent bins, you are actually solving for frequency, subject to the requirements of the technique.  Some code (DevC++ compiler) for your specific problem, followed by the programs' output is shown below (explanation to follow):

#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <cmath>
using namespace std;
void fft_recur(long, double* r, double* i);
int main (int nNumberofArgs, char* pszArgs[ ] ) { // arguments needed for Dev C++ I/O
const long N = 1024;   double sample_rate = 44100., r[N] = {0.}, i[N] = {0.}, w[N] = {0.} ;
long n;      double  t, twopi = 2*3.141592653589793238;
double amp_sq[N], F, fbin, amp, power, c = 2.3;

for (n = 0; n < N; n++)  {  // generate test data
    t = twopi*n/sample_rate;
    r[n] = 1.*cos(16351.*t)+1.*cos(16996.*t)+1.*cos(17211.*t)+1.*cos(17426.*t) ;
} // end for

//Gaussian window as computed in [MCEACH94] program
double z;
for (n = 0; n < N/2; n++) {
    z = (.5+(N/2)-1.-n)*(twopi/2.)/N;
    w[n] = exp(-c*z*z);
    w[N-1-n] = w[n];
}

//multiply input data by window points
for (n = 0; n < N; n++)  {
    r[n] = w[n]*r[n];    // we're overwriting r[n] with the windowed r[n]
} //end for

fft_recur(N, r, i);       // note: fft is not scaled

// compute amplitude squared
for (n = 0; n < N/2+1; n++)  {
    amp_sq[n] = r[n]*r[n]+i[n]*i[n];
} //end for

// compute natural log of FFT outputs
for (n = 0; n < N/2+1; n++)  {
    r[n] = .5*log( amp_sq[n] ) ; // r[n] is now the natural log of the amplitude of FFT bin 'n'
} // end for                     // .5*log(A squared) = log(A)

cout<<"\n\nfrequency and amplitude estimation results\n\n n        magnitude      frequency       amplitude\n\n";
for (n = 370; n < 415; n++)  {
    fbin = n+.5+.5*c*(r[n+1]-r[n]);   // r[n+1] and r[n] are ln(b) and ln(a)
    F = (sample_rate/N)*fbin;
    amp = sqrt(c)*exp(r[n]+((fbin-n)*(fbin-n))/c)*(2.0*sqrt(twopi/2.)/N) ;
    printf("%2d\t%9.5f\t%9.5f\t%9.5f\n",n,amp_sq[n],F,amp);
} //end for

system ("PAUSE");
return 0;
} // end main
//******************** fft_recur ***********************
void fft_recur(long N, double *r, double *i)  {
long h, i1, j = 0, k, i2 = N/2, l1, l2 = 1;
double c = -1.0, s = 0.0, t1, t2, u1, u2;

for (h = 0; h < N-2; h++) {    // ***** bit reverse starts here ***
    if (h < j) {
       t1 = r[h]; r[h] = r[j]; r[j] = t1;
       t2 = i[h]; i[h] = i[j]; i[j] = t2;
    }
    k = i2;
    while (k <= j) {
          j = j - k;      k = k/2;
    }
    j = j + k;
}   //****** bit reverse done ******

for (k = 1; k < N; k = k*2) {
    l1 = l2;     l2 = l2*2;
    u1 = 1.0;    u2 = 0.0;
    for (j = 0; j < l1; j++) {
        for (h = j; h < N; h = h + l2) {
            i1 = h + l1;
            t2 = (r[i1] - i[i1])*u2 ;
            t1 = t2 + r[i1]*(u1 - u2) ;
            t2 = t2 + i[i1]*(u1 + u2) ;
            r[i1] = r[h] - t1;
            i[i1] = i[h] - t2;
            r[h]  = r[h] + t1;
            i[h]  = i[h] + t2;
        } // end for over h
        t1 = u1 * c - u2 * s;
        u2 = u1 * s + u2 * c;
        u1 = t1;  //x = u1 - u2;    y = u1 + u2;
    } // end for over j
    s = - sqrt((1.0 - c) / 2.0);
    c =   sqrt((1.0 + c) / 2.0);
} // end for over k

} // end fft

frequency and amplitude estimation results

 n        magnitude      frequency       amplitude

370       0.01156       15959.53142       0.00065
371       0.01328       16002.86947       0.00070
372       0.01542       16046.20441       0.00076
373       0.01810       16089.42371       0.00082
374       0.02138       16136.92108       0.00095
375       0.03020       16209.62311       0.00211
376       0.14116       16370.55303       3.20360
377      77.00512       16349.81449       0.95808
378     3193.91932      16351.21160       1.00396
379     24622.71009     16350.90700       0.99982
380     32939.04014     16351.10559       0.99850
381     7803.22610      16350.67285       1.01020
382     319.11953       16353.81879       0.87280
383       2.60322       16324.16039       7.57917
384       0.00113       16605.10228       0.00051
385       0.00724       16606.09917       0.00052
386       0.00851       16650.89175       0.00058
387       0.01072       16694.51462       0.00065
388       0.01381       16737.66957       0.00074
389       0.01787       16785.61197       0.00090
390       0.02814       16862.78585       0.00251
391       0.17568       17013.74935       2.83597
392      85.58946       16994.84074       0.95935
393     3413.62331      16996.22071       1.00408
394     25288.20406     16995.88948       0.99974
395     32472.56144     16995.79357       1.00043
396     7296.89833      16957.17339       4.09340
397      60.55471       17218.13259       1.24775
398     3330.74781      17212.33274       1.01810
399     25463.82887     17210.95691       0.99978
400     32350.54927     17210.69761       1.00169
401     7144.91966      17170.91626       4.35710
402      55.60789       17436.25628       1.40301
403     3423.65695      17427.15709       1.01474
404     25642.98140     17426.03425       1.00002
405     32244.80235     17425.93874       1.00073
406     7095.46971      17426.35995       0.98918
407     278.99401       17421.20777       1.27242
408       1.56515       17494.25360       0.02622
409       0.02946       17619.31213       0.00091
410       0.01521       17677.62446       0.00071
411       0.01452       17719.83432       0.00069
412       0.01340       17762.64014       0.00066
413       0.01224       17805.61294       0.00063
414       0.01113       17848.67454       0.00060

First, 4 tones of the input are generated (I use twopi when computing time because it's easier than writing it 4 times into the cosines).  Then, the Gaussian window is generated, after which it is applied to the input (by the way, the window shown here, which is from reference_1, is slightly wrong – I don't want to use the 'corrected' version of it because I don't want to write a page describing how I derived it).  Then we take a 1024 point FFT, followed by the (non-normalized) calculation of amplitude squared.  Then comes the calculation of the natural log of each FFT bin.  Finally, I compute the frequencies and amplitudes via the equations in reference 1 and reference 2above.

As you can see from the results following the code, we have an isolated peak between bin 379 and 380 from the first input frequency, which is estimated as 16350.90700 and an amplitude of .99982.  We can also see that the other 3 inputs exhibit a little 'peak merging' from bins 392 to 407, but the peaks and associated estimated frequencies and amplitudes are quite accurate.

The results could be even more accurate if I had: 1) adjusted parameters, 2) used the 'correct' window, 2) used a non-recursive FFT, or 3) used one of the variations of the technique (which I am not going to show here).

You might also notice that, as your value of the parameter 'c' gets smaller ('c' is related to the 'beamwidth' of the Gaussian window), you can actually improve some of your estimates, but it can cause some results to blow up due to numerical sensitivities.  The effective 'Q' or 'gain' of the Gaussian ratio technique is on the order of 5,000 – 10,000, and it can be much higher if you know how to adjust the parameters properly.

Of course, if your SNR is too low, you may have to increase N to get more integration time, subject to the limits of signal stability.


Defining Model in Kalman for Getting Position Using I am stuck at modeling a system model, i.e. getting my state vector and input vector for navigating just using navaid and ins (tactical). My guess is that position is my only state vector and INS output (acceleration or velocity) is input vector. What is state vector and what is input vector in my case? As I just need 2 position which sensor must be used for update and which must be used for predict. Should I use 1d 2 kalman implementation for lat and lon or no etc?

In general, for position-related Kalman filters, you want your state vector to contain $x$, and $\dot{x}$ (location and velocity) components. See, for example, the Wikipedia page or this question and answer here.

If you have a measure of velocity, then it can certainly also be an input to the Kalman filter.


Why discretize a continuous transition matrix in Kalman Filter? In the Kalman filter toolbox at http://becs.aalto.fi/en/research/bayes/ekfukf/cwpa_demo.html the example code shows that a function lti_disc is called, which is essentially a matrix exponential function. What is the purpose of such matrix exponentiation? Given a really simple transition model such as, e.g.

F = [ 1 0 ;
      0 1 ] ;

and doing the matrix exponentiation with lti_disc will always give a different transition matrix A, depending on the values of the Q process noise covariance, which is also an input to the lti_disc function.

Have you taken a look at the documentation - section 2.2.4 discusses a Linear Kalman filter model that is very similar to the one you described. From that example you see that:

The resulting $\bf{A}$ matrix does not depend on the process noise
The $\bf{A}$ matrix only depends on the size of the time step.
The $\bf{Q}$ matrix only depends on the size of the time step, and the value of the continuous noise covariance, ie. diag(q) in the example

You can use discrete noise kinematic models, they are probably a bit easier to understand. The continuous noise model assume that the noise is applied throughout the time step interval i.e. from time $t_0$ to $t_1$. The discrete model assumes the noise is applied instantaneously at times $t_0$ and $t_1$ (and so on).

See chapter 6 -"Estimation for Kinematic Models" in the book - "Estimation with Applications to Tracking and Navigation" (Bar-Shalom, Rong Li, Kirubarajan) for a derivation of the continuous and discrete time versions of the state transition models.


Using Spreading spectral with Maximal likelihood As asked here, HERE  if we have the signal

x =
[0.7 + 0.7i;
0.7 - 0.7i;
-0.7 + 0.7i;
-0.7 -0.7i];

Which was spread over code c and transmitted over channel H whose dimension is [4x4], so the convolution of signal after spreading will become:

r = reshape(H*reshape(x,4,[]),[],1);

comparing it with signal without spreading, it was

 r = H*x;

That was explained well in the above link.

My question, suppose I am using Maximum likelihood estimation, in case if we didn't spread the signal, we will check the likelhood compared with the channel H, but what's about after using spreading ? how will become the channel ? It supposed to be a vector of [16x1], is that right ? but how will it be ?

thank you!

H was not reshaped, as you see in your command r = reshape(H*reshape(x,4,[]),[],1);, you reshaped the data itself.

In that case, you are going to add noise, then using ML estimation based on the received data. So that, H will be H without changing, what will be changed is the received data, you can reshape it similar to that way in the transmitter reshape(y,4,[]), where y is the received data, then reshape the results again into [16 x 1].

Good luck


Output signal as convolution of impulse response and step signal Hi All: This question is kind of an add on question to the gorgeous answer by Juancho provided at this link:

How does this "simple filter" work?

The answer explains how the limit of exponential smoothing given an impulse response $x_0$, (so that new $x_i$'s are not arriving. only $x_0$ arrives) can be computed by convolving the impulse  with a step signal.

The difference between the framework at that link and my framework is that  actual series is not driven by an impulse response. By this I mean that the same process, occurs over and over because a "new" $x_i$ is arriving during each period so one could think of there being an impulse response as described but occurring in each period.  So, there is a step signal convolved with the first impulse response. Then, in the next period, another step signal is convolved with a second impulse, etc, etc. In other words, now there is $x_1$, $x_2$, $x_3$, ..... $x_n$, where each $x_i$ gets operated on in the same manner as described by Juancho.  Assuming that my setup is clear, my question is: is there a simple recursion for this non-impulse response case. Does my case have a name ?  It seems like there should be a way to use recursion  but I can't figure it out. Thank you very much for any references, answers, books, whatever.

It is a normal IIR filter, and a very simple one at that.

Everything that can be described with a linear time invariant difference equation is an IIR filter and can be designed and analyzed using well established mathematical methods.

You calculate the output of the filter for ANY given input (impulse, step, the Dow Jones, Wagner's Nibelungen, etc.) with the equation
$$y[n] =\sum_{i=0}^{N}b_i\cdot x[n-i]-\sum_{i=1}^{N}a_i\cdot y[n-i]$$

In your case you have $b_0=0.2$ and $a_1 = -0.8$ with all the other a and b coefficients being zero, so the equation simplifies to
$$y[n] = 0.2 \cdot x[n] + 0.8 \cdot y[n-1] $$

That is indeed a simple first order low pass filter with a -3dB point at a normalized frequency of 0.0355 or thereabouts. It's not a very good low pass filter since it doesn't go to zero at the Nyquist frequency. A somewhat better choice would be

$$y[n] = 0.1 \cdot x[n] + 0.1 \cdot x[n-1] + 0.8 \cdot y[n-1] $$

There are plenty of smoothing algorithms that are not IIR filters but your specific example is just a simple first order low pass.

EDIT: based on the question in the comment
The equation above describes the operation for ALL input signals, not just an impulse.
Let's consider the input signal x = [1.0 1.1 1.05]; So we have x[0] = 1.0, x[1] = 1.1, x[2] = 1.05, and x[3...] = 0; You would calculate the output as follows

y[0] = 0.2*x[0] + 0.8*y[-1] = 0.2; // since y[-1] = 0;
y[1] = 0.2*x[1] + 0.8*y[0]  = 0.38;
y[2] = 0.2*x[2] + 0.8*y[1]  = 0.514;
y[3] =            0.8*y[2]  = 0.4112;     // since x[3] = 0;
y[4] =            0.8*y[3]  = 0.32896;
etc.


Zero Centering Audio for Forward FFT I'm trying to work out the best way to zero centre an audio clip for a forward FFT.   I've been trying rotate my audio so that a peak exists mid way (N/2).   This seam a little naive and I thought I had read that you can zero pad the centre of the audio clip in order to achieve the same result.  Am I mad?  what is the best way to zero centre an audio clip for Forward FFT

Eureka! You can get a zero centered waveform of an audio clip by applying a window, i used hann (cosine) window and then inverting the left and right sides.


how can i get to know the frequency of a speech signal? What is the actual meaning of frequency? what if the signal is aperiodic??
In general, we say that frequency means the repeatition of an event in unit time. But what about a speech signal? how can we say that a speech signal is having 4KHz Bandwidth(say)?

A speech signal can contain a lot more than 4kHz bandwidth, however experiments showed that if you removed all but the spectrum below 3kHz (thus allowing the use of cheaper circuitry), it was still reasonably understandable by another human.


How to decide what subcarrier allocation should I use in OFDM? For those familiar with the liquid-dsp library, or with OFDM in general, I am hesitating what subcarrier allocation should I use in a particular wireless system.

By subcarrier allocation I am refering to how many null subcarriers, pilot subcarriers and data subcarriers should I use, and in what order. Let's say that the total number of subcarriers is forced to be 64, due to channel bandwidth limitations.

How is this related to the perfomance of the system?

The liquid-dsp library provides a default subcarrier allocation, but I would like to use an allocation tailored to my system's requirements. I am interested in making my system more robust to Doppler effect (and thus to frequency shifting).

Is it right to assume that more null subcarriers or more pilot subcarriers will improve the OFDM performance under Doppler effect?

As is typical, the answer is that it depends. In most cases, you're implementing a system that someone else has designed already, and you just need to match their design. However, if you're in the situation of designing a system of your own, then there are a few high-level considerations:

Data subcarriers are the ones that actually communicate information through the channel. Ultimately, the goal of any communications system is to transfer information with as little energy and time expended as possible, so maximizing the number of data subcarriers used in each symbol helps to reach the second of the two goals.
Null subcarriers don't do anything for information transmission; instead, they are used to shape the spectrum of the transmitted signal. Therefore, we would want to minimize the number of null subcarriers that are in use, as they reduce the information rate of the system. They are often used for two reasons:

Allowing some room for analog anti-aliasing filters. In a practical system, you'll need to leave some portion of the upper end of the frequency band near the Nyquist rate empty to allow for rolloff of pre-ADC and/or post-DAC filtering. Typically, you'll see the upper ~10-25% of the band unused.
Exclusion of the center bin in the band. Common, low-cost direct conversion receivers usually exhibit a high level of interference at band center (due to ADC DC offset after the conversion step). Therefore, it's common in wireless standards that are aimed at low-cost markets to just place a null subcarrier in that slot to avoid the potential interference.

Pilot subcarriers are used in the equalization process. Again, since their contents are known to the receiver a priori (or else they wouldn't be useful as pilots), they don't convey any information through the channel. So, one would want to minimize the number of pilot subcarriers in use. However, you want to have enough of them so that you can effectively estimate the response of the channel across the entire frequency band. This is why you'll often see them spaced across the OFDM symbol. Also in many cases, the pilot positions will rotate over time, allowing the receiver's equalizer to periodically get a view of the entire channel response.

Note that the number of subcarriers doesn't define the channel bandwidth. OFDM transmitters are often implemented using an IFFT structure so you can lay out the subcarriers in the frequency domain first, then transform to the time domain to get a waveform that you can apply to the channel. The sample rate that is used on that waveform, along with the subcarrier allocation and FFT size, is what defines the bandwidth of the transmitted signal. You can increase the FFT size arbitrarily while still maintaining the same signal bandwidth (you will just end up with more lower-rate subcarriers that are spaced more finely in the frequency domain).


Do Spread Spectrum technologies combat time-selective fading of a channel? The Spread Spectrum technologies such as DSSS and FHSS are successful against frequency-selective fading because the original signal is spread to a larger bandwidth and different frequencies are used a different times.

Do these technologies also mitigate the problems imposed by time-selective fading?

I hadn't heard the term "time-selective fading" before, and after a bit of research it appears that it refers to the effects of doppler spreading.  Given that, the answer is "not really".

DSSS is typically extremely sensitive to carrier offset since the pseudo-random chipping sequence is usually quite long.  The long time means that a little carrier offset can lead to destructive interference.  This is dealt with by transmitting the data with a pilot.  It is easy to search for the pilot across multiple carrier offsets by correlating using FFT techniques.

FHSS does not really have any particular advantages or disadvantages that I am aware of when it comes to handling doppler.


Scipy fourier transform zero frequency spike (from DC offset) - de-meaning and hanning window have no effect I am trying to plot the FFT of essentially a random signal that has a non-zero mean shown below.

The FFT of the signal is peaked over the zero frequency which usually indicates a DC offset. Although I have already demeaned and applied a hanning window to the data.

I am not sure what I am doing wrong here. I have also tried plotting from [1:] to ignore the zero frequency and applying a highpass filter which have had no effect. My code is below.
# a = signal.detrend(hhe_trim[0].data, type='constant')
a = hhe_trim[0].data.astype(np.float64)
a -= np.mean(a)

win = signal.hann(174001)
dt = 0.01
n = 174001

X = fftpack.fft(a*win, n=n)
freqs = fftpack.fftfreq(len(a), dt)
plt.plot(freqs, np.abs(X))
# plt.yscale("log")
plt.show()

My data has 174 001 samples and a sampling rate of 100Hz. Any help would be appreciated!

It doesn't happen with a random signal.  Your signal must have low frequency content around 0 Hz that shows up even after you've nulled out 0 Hz itself?
import numpy as np
from scipy import signal
from scipy import fft
import matplotlib.pyplot as plt

a = np.random.randn(174001)
a -= np.mean(a)

win = signal.hann(174001)
dt = 0.01
n = 174001

X = fft.fft(a*win, n=n)
freqs = fft.fftfreq(len(a), dt)
plt.plot(freqs, np.abs(X))
# plt.yscale("log")
plt.show()

If you zoom in on 0 Hz, I bet you can see that the 0 Hz spike itself goes away when you null out the mean, but the content around 0 Hz stays the same.


Phase locked loop for three phase grid I have been implementing a control software where one of its core functionalities is the ability to synchronize with the three phase grid. Basically I need to implement some form of the phase locked loop (PLL) suitable for three phases. I have found that there are several posibilies how to do that. One of the simplest I have found is basically based on calculation of the time integral of the phase voltages. My first idea how to do that was to use the trapezoidal integration rule in recurent form
$$I(k) = I(k-1) + \frac{T}{2}\cdot\left[x(k-1) + x(k)\right].$$
I have implemented the above written difference equation in C and I have tried to pass a sinewave into it
    #define PI 3.14
    #define N 64

    float x, y;
    for(int k = 0; k < N; k++){
        x = sin(k*2*PI/N);
        y = integrate(1.0, x);
    }

    float integrate(float T, float xk)
    {
        static float intk_1 = 0; // I(k-1)
        static float xk_1 = 0; // x(k-1)

        float intk = intk_1 + T/2*(xk_1 + xk);
        intk_1 = intk;
        xk_1 = xk;

        return intk;
    }

Based on the outcomes of my experiment I have realized that in case the integration process starts from the point where the sinewave has zero value the resulted integral has large offset

In case the integration process starts from point where the sinewave has its positive maximum the integral has zero offset.
I need to have the integral with zero offset to be able to estimate the phase and frequency of the three phase grid. The obvious solution could be to somehow detect the positive maximum of the sinewaves and then start integration process. The problem is that the frequency of the three phase grid fluctuates in some tolerance band. Can anybody give me an advice how to do the integration process so that the resulted integrals have zero offsets?

I recommend a standard dq-PLL, check this reference https://vbn.aau.dk/ws/portalfiles/portal/273236528/PLL_Review_RSER.pdf

There are 2 required integrators. For algebraic loop reasons, the last integrator, the one that converts the angular frequency to the phase, should be a forward-Euler integrator. The integrator in the PI controller, can be a trapezoidal integrator if you want.
Edit : Ts = Sampling Period


What assumptions should be used to invert spectrum into time domain data? This question comes from a realistic problem that how to obtain a time domain wind pressure from known wind spectrum. Or similarly, gets rail surface curve from spatial rail spectrum.

Obviously, this question does not have unique solution, since phase information is not given. Usually we assume that phase of the time/spatial domain data in uniformly distributed in ($-\pi$, $\pi$], and using a random distributed phase to do inverse Fourier transform to obtain the time/spatial domain data. However, I feel not so comfortable to assign such rapid changing phase to the spectrum. I wonder whether there are some more reasonable assumptions to get these jobs done?

BTW, question ifft of signal without phase information discussed some strategies on how to do this job. @LutzL suggested Human perception can be used to help determine the phase information. But I cannot understand how Human perception is working. @Hilmar said minimum phase may be a choice, but I do not know whether there is a physical meaning of minimum phase. Is there any discussion to show minimum phase is a good choice for physical phenomena, like wind.

If you know anything at all about the time-domain signal, for example that it's piecewise linear, or should never be negative, or should look like a series of spikes with exponentially decaying tails, you could define a cost function that measures how well a proposed signal meets those criterea.  How to do that is hard to describe in a generic way.  Then feed that cost function into your favorite multi-variable optimizer.  The phases of the Fourier amplitudes are the variables to be optimized.   Try a variety of initial values, beware local minima, and all the other usual cautions for optimization apply.


How does sampling and subsequent resampling affect aliased harmonics? If I sample a signal with a harmonic higher than half my sample frequency I can predict the aliased frequency easily as

$$|F_h-kF_s|$$

But what if after sampling the signal, I re-sample at a lower rate (decimate)?

Does the lower sampling frequency simply trump the original higher frequency or will there be mixing or modulation effects due to the two sampling frequencies that change or create new aliased harmonics?

If so, is there a way to predict what the aliased harmonics will be?

As a simple example assume I sample a signal containing a $7$ kHz harmonic at $3$ kHz, then decimate by taking every third sample, so effectively re-sampling at $1$ kHz.

The alias behaves exactly the same, in every sense, as if it were an actual sinusoidal component at the alias frequency.

In your example, the 7 kHz harmonic sampled at 3 kHz is actually a 1 kHz sinusoid, and decimating it is exactly the same as decimating a 1 kHz signal. In this case, you end up with a DC value.


What is the inverse Laplace transform of squared denominator term? Referring to the image below, what would the inverse Laplace transform be? I can't seem to find any tables that include this case.

This is quite straightforward to solve. Either just use this table where you can directly find the corresponding result, or "derive" it yourself with very basic knowledge of the Laplace transform.

You should know that

$$\mathcal{L}\{u(t)\}=\frac{1}{s}\tag{1}$$

Integration is equivalent to multiplication with $1/s$, so integrating $(1)$ gives

$$\mathcal{L}\{t\cdot u(t)\}=\frac{1}{s^2}\tag{2}$$

Replacing $s$ by $s+a$ is equivalent to multiplication with $e^{-at}$. Consequently, we have

$$\mathcal{L}^{-1}\left\{\frac{2}{(s+2)^2}\right\}=2\,t\,e^{-2t} u(t)\tag{3}$$


Plot the frequency response of MTI delay line canceller using Octave I try to plot the frequency response of the delay line canceller (FIR filter which used for MTI).

The delay line canceller has the following structure:

              +-----+
              |     |
x [k] >---+---|  T  |
          |   |     |
          |   +-----+
          |      |
          |   +-----+
          |   |     |
          |   |x -1 |
          |   |     |
          |   +-----+
          |      |
          |   +-----+
          |   | ___ |
          \---| \   |---> y [k]
              | /__ |
              +-----+

Its frequency response is well known (Pg 20) and equal:

$$H(\omega) = 2 \cdot \left|\sin (\frac{\omega \cdot T} {2})\right|
\tag{1}$$

I also try obtain frequency response of the canceller using algorithm which described on this question:

$$y [k] = x [k] - x [k - 1] \tag{2}$$

$$Y (z) = X (z) - X (z) z ^{-1} \tag {3}$$

$$ H(z) = 1 - z^{-1} \tag{4}$$

Now I try to plot (using GNU Octave) both of responses (1) and (4).

w = linspace (0, 2 * pi, 100);
z = exp (-j .* w);

H_z = 1 - z .^ -1;
H_w = 2.0 * abs (sin (w / 2) );

hold ("on");
plot (w, H_z, "1", "linewidth", 2);
plot (w, H_w, "2", "linewidth", 2);
title ("Frequency response of the delay line canceler.");
set (gca, 'XTick',     0: pi / 2: 2 * pi)
set (gca, 'XTickLabel',{'0', 'pi / 2', 'pi', '3 pi / 2','2 pi'})
xlabel ("Angular frequency.");
ylabel ("Magnitude.");
legend ("H (z)", 'H (\omega)');

I expect that they will be the same, but they are different.

Where is my mistake?

P.S. If I add modulus for (4) (like: H_z = abs (1 - z .^ -1);) they will became the same.

Magnitude |H(z)| is always defined as an absolute value of your transfer function H(z) values which are complex. What you are plotting is the vector that contains complex numbers. You need to take their absolute value and that will give you magnitude.

If you call the angle function then you are going to get phase. Transfer function H(z) is always complex. Try to print values of real(H_z) and imag(H_z), then you will see real and imaginary values of $H(z)$. Now magnitude |H(z)| is sqrt(real * real + imag * imag) for every sample...


RMS from calculation from PSD using Welch's method I've used Welche's method on a vibration signal to get the averaged FFT which I want to used to calculate the PSD from which I can calculate the RMS.

Within the Welche method i've use Hann window with zero overlapping. The following post states, that using
windowing couses distortion in the energy, which has effect on the RMS:
https://community.sw.siemens.com/s/article/window-correction-factors .

However not using windowing might lead to spectral leakege.

My question is, that what is the correct way of dealing with this problem? Shall I use uniform windowing, or shall I use correction factors on the FFT?

Also, you can use the PSD for estimation of the RMS, the easier way for estimation of the RMS is a direct calculation.

In case, you still want to use the PSD for estimation of the RMS, you should normalize the values by the square of the average value of the window. For the Hann window, the normalization factor will be $\frac{1}{4}$.

update
The question is a duplicate from RMS calculation in frequency domain after windowing


pre-emphasis with scaling? I am working on speech synthesis with LPC and I originally implemented a pre-emphasis filter:

+(void)processBuffer:(Buffer *)buffer {
    float alpha = -0.9375f;
    for (int i = 1; i < buffer.size; i++) {
        buffer.samples[i] += buffer.samples[i - 1] * alpha;
    }
}

And found that my output audio was distorting, so I adding some code to scale the values, and so now it looks like this:

+(void)processBuffer:(Buffer *)buffer {
    float preEnergy = [self energyFor:buffer];

    float alpha = -0.9375f;
    for (int i = 1; i < buffer.size; i++) {
        buffer.samples[i] += buffer.samples[i - 1] * alpha;
    }

    [self scaleBuffer:buffer preEnergy:preEnergy postEnergy:[self energyFor:buffer]];
}

+(void)scaleBuffer:(Buffer *)buffer preEnergy:(double)preEnergy postEnergy:(double)postEnergy {
    float scale = sqrt(preEnergy / postEnergy);

    for (int i = 0; i < buffer.size; i++) {
        buffer.samples[i] *= scale;
    }
}

+(float)energyFor:(Buffer *)buffer {
    double sum = 0.0;

    for (int i = 0; i < buffer.size; i++) {
        sum += buffer.samples[i] * buffer.samples[i];
    }
    return sum;
}

I am finding that the speech I am generating still doesn't sound very good, and so I am just wondering if my implementation for this filter is ok?  Does anyone see anything I am doing wrong or does this look ok?

The most glaring problem is that buffer.samples[0] is never filtered because the for loop goes from 1 to size-1.  This is an inherent problem in the implementation of FIR filters in software where you are dealing with buffers of data rather than an infinite stream. In the more general case, where your filter has N states you would need to store N-1 states so that they can be used on the next invocation of the filter. Since your filter only has two taps and the implementation is rather more simple than a general solution would be.

// state is a member variable of the class initialized to zero.
float alpha = -0.9375f;
buffer.samples[0] += state * alpha;
for (int i = 1; i < buffer.size; i++) {
    buffer.samples[i] += buffer.samples[i - 1] * alpha;
}
state = buffer.samples[buffer.size-1]; // store for the next time.

Scaling

Your filter has close to 2x gain in the passband that might be causing you problems. I suspect samples need to stay in the range of -1 to 1. Imagine samples was a stream of [...,1,-1,1,-1,...] (this happens to be a sine wave at the sampleRate/2).  Your filter output is going to be [...,-1+1*-0.98,1+-1*-0.98,...] which is approximately [...,-2,2,-2,2,...]

The strategy of scaling each consecutive buffer by the filter gain for that buffer is a nice idea but it has the unfortunate property of producing a different scaling factor every time you call it. This would likely introduce a transient between each buffer.

It is better to analyze the filter and choose a constant at design time to scale by.

In your case the would be 1/2.

// state is a member variable of the class initialized to zero.
float alpha = -0.9375f;
buffer.samples[0] += state * alpha;
buffer.samples[0] *= 0.5;
for (int i = 1; i < buffer.size; i++) {
    buffer.samples[i] += buffer.samples[i - 1] * alpha;
    buffer.samples[i] *= 0.5;
}
state = buffer.samples[buffer.size-1]; // store for the next time.

Of course if you wanted the gain then disregard this part.


How do I understand Fourier descriptors more visually and intuitively? I read the book Image Processing, Vision and Machine Vision and find the concept Fourier descriptors hard to understand, although literally its derivation is somewhat reasonable. Can anyone give me a detailed picture about what's really going on about it? Thx.

Fourier Descriptors (FD) are a way of representing the shape of a closed curve at varying forms of detail. In the following discussion, I am assuming that you are familiar with the one-dimensional Discrete Fourier Transform (DFT).

FD are essentially the application of the Fourier Transform to decompose the form of the curve into a sum of shapes that look roughly like circles whose radius has been modulated. In fact, the equivalent of the zero frequency component (or DC) is a perfect circle with a radius equal to the average distance to the centroid of the curve.

We can continue thinking in terms of superimposed "squiggly" circles that result to a complex curve, or, we can "transfer" the curve to a more familiar setting of a "waveform".

To do this, first estimate the centroid of the curve by calculating the arithmetic mean of the curve's coordinates.

Now imagine that you have the familiar $x,y$ coordinate system but "wrapped" around the curve. The $x$ axis now corresponds to some angle ($\theta$) around the curve and the $y$ axis now corresponds to the distance between a point on the curve and the centroid at angle $\theta$.

If you apply a polar to cartesian coordinate transform that takes the pair of $\theta, r$ coordinates to the familiar $x,y$ space, you will have unfolded that axis to its more familiar "rectangular" shape.

So, now, you have a very simple one-dimensional waveform whose overall shape you can analyse by using (usually) just a few FDs.

Just two further points that might not be immediately obvious:

The curve has to be closed because the Fourier Transform is defined over periodic waveforms.
The way to "input" the curve to the Fourier Transform is by taking its $x,y$ points to form a complex number (e.g. $z= x \cdot i +y$). After this, you apply the Fourier Transform on that set of complex numbers (i.e. a series of $z_k$ where $k \in \left[0 .. N-1\right]$ and $N$ being the total number of $x,y$ points that the curve is composed of.

Hope this helps.


Converting Google's text-to-speech voice to, say, Justin Bieber's voice We are trying to convert text to an artist's voice. For example, we want input text in Justin Bieber's voice. We are converting text to mp3 files in Python using gTTS (Google Text-to-Speech). We can then change the pitch or speed to some extent.

How can we convert gTTS's mp3 file to Justin Bieber's voice? Any pointers will be helpful. Thanks :)

How we should go about from converting Google Text to Speech's m3 file to Justin Bieber's voice?

Ask Mr. Bieber to repeat what your MP3 says.

I think letting one speech synthesizer first speak something and then hoping to convert it is a bad approach (as you can see if you consider how complicated that would be – first, you need to recognize in some low-level way the things that the speech synth said, then re-synthesize them with a Bieber voice synth).

Instead, "just" implement a Bieber voice synth. That will "only" require you to have a model of the target voice, instead of another model for utterance recognition.

I'm not deep enough into voice synthesis, but common voice models have a metric ton of parameters that you can tweak – maybe you can tweak a male voice model until it resembles Mr. Bieber's voice?

I honestly think you might be widely underestimating the nuances and number of parameters that human voice has.

For example, record you saying something, and then record someone from two generations older than you (or younger than you, whichever is easier to find) say the same thing. Compare rhythm, compare amount of pauses, compare stress on syllables and words – you can't just go ahead and convert your voicing into someone else's voicing just by means of simple transformation.


How do Anti-Aliasing-Filters in audio signal processing work? I am a little confused about the use of anti-aliasing filters. As on how I am confused, please consider the following task:

Consider a simple sampling rate conversion system with a conversion rate of 4/3. The system consists of two upsampling blocks, each by 2, and one downsampling block of 3.

The frequency spectrum has periodic repetitions at integer multiples of the sampling frequency. Therefore, upsampling creates additional - unwanted - spectral images. These can be cancelled by an anti-imaging filter:

The spectrum within the blue dotted rectangles should be my output Y_1. This output will be upsampled a second time, which I think results just having 4 spectrum repetitions until 2*pi.

But, how does downsampling work? I am given the following figure from my lecture notes:

So, apparently before downsampling we apply the anti-aliasing filter, which will result in the part of the spectrum within the blue-dotted shapes. Now, we increase Omega by a factor of 2, such that.. I don't know. What exactly happens here ?

I think you make this needlessly complicated. Let's do a specific example: let's say we want to go from 48 kHz up by 4/3 to 64 kHz and your signal bandwidth is 20 kHz.

In the first step, you would up-sample by a factor of 4 by inserting 3 zeros between each sample. Your new sample rate is 192 kHz. Your original spectrum is preserved and you get mirror spectra centered around 48 kHz & 96 kHz.

To eliminate this you need a lowpass filter with cutoff of 20 kHz at 192 kHz sample rate and sufficient attenuation at 24 kHz.

Finally you just down sample the result by factor of three by simply throwing away every second and third sample. This works without aliasing, since the lowpass filter has already eliminated all content that would alias.

Regarding your graphs: in general I would avoid drawing anything past the Nyqusit frequency. It's just periodic repetition and doesn't add any useful information and it can get confusing.


FEC decoder processing delays I am trying to implement a burst demodulator. I have catered bits in premable for carrier and timing recovery and for fec also
in the recieve side if i am using viterbi decoder/TPC coder will the processing add any extra processing time. So if my burst time is say 1 ms. Is it possible to decode inside the 1ms or will there be any processing delays?

Well, yes, everything that you can calculate needs some time to calculate. Viterbi is quite well-behaved with respect to that, since you can decode up to the last couple of bits (only the last constraint length bits are relevant to the current decision), and you can already prepare the path metrics while receiving sequentially, and only need to do the backtracking after reception of the last code bit.
However, in most practical systems, you'd want some sort of interleaving. Also, you'll have filtering and all sorts of operations with inherent delay before you even decode. So, there's group delay and processing latency in your system anyway.
That, together with the fact that Viterbi is not an overly complex method and typically not used with convolutional codes with very large constraint lengths, and that 1 ms is a rather long time for a modern processing system:
You are probably barking up the wrong tree. You haven't defined acceptable latency yet, nor the length of your burst, nor the degree of error correction you need.
Start with that, and figure out (we'll be glad to help) whether processing delays will be a problem or not, before worrying about your decoder.


Modeling ADC in Active Noise Cancellation I'm tinkering in Matlab with a problem that's very similar to active noise cancellation.  In the literature, the secondary path is described as the transfer function from the output of the adaptive filter to the error input sensor.  The algorithm needs to model this path to obtain good results.

However, there is also the transfer function from the noise source to the adaptive filter.  The ADC will always have some transition band below the Nyquist rate.  When I simulate the ADC in Matlab, the adaptive filter tries to invert the ADC response as well as model the system I'm trying to cancel.  In fact, if the system is purely a delay, this is exactly the inverse system identification problem.

In some applications, the ADC can be ignored because the signal may already be digital, but I don't think that is ever the case with ANC.  Why is this not mentioned as being a problem?  Is it because, in practice, the impulse response of the ADC is short compared to the system response?

I am also working on ANC, with my understanding may be I could provide some clarity for you.

I could not see your point exactly on ADC modelling, but as far as I know this ADC transfer function will be covered in secondary path modelling.

  "However, there is also the transfer function from the noise source to
  the adaptive filter."

Do you mean modelling of the reference input microphone, as far as I can see only transfer function involved between noise source to adaptive filter is reference input microphone function.

The way you are using secondary path description is also a bit confusing to me.

But here is what I know from this paper page 945

Here P(z) is the unknown primary acoustic path between reference input sensor and error sensor.

Below text directly quoted from the reference provided

  Summing junction represents the acoustical superposition in the space from cancelling loudspeaker to error microphone, where the primary
  noise is combined with the output of the adaptive filter.

  Therefore, it is necessary to compensate for the secondary-path
  transfer function from y(n) to e(n) which includes the digital-to-analog (D/A)
  converter, reconstruction filter, power amplifier, loudspeaker,
  acoustic path from loudspeaker to error microphone, error microphone,
  pre-amplifier, anti-aliasing filter, and analog-to digital (A/D)
  converter.

The compensation as you might have already known is done by offline or online modelling techniques.

If you want to include your primary path modelling as well i.e reference microphone TF, you can do that or if you are considering feedback effect, then again in the feedback compensation this reference microphone transfer function is included.

Please correct me if I am wrong.
Thanks for bringing this up!


Formula for 2D Image Filter for Scratch Detection Consider Multidirectional Scratch Detection and Restoration in Digitized Old Images Research Paper.

In section 4.1 (Preprocessing), they gave us a formula of the kernel,

$$ H(u, v) = \frac{1}{1+0.414 \times \left(\sqrt{\frac{u^*}{D_h}+\frac{v^*}{D_v}}\right)^{2n}} $$

Where,
\begin{split}
t_x & = \mathrm{center} \times \cos(\theta) \\
t_y & = \mathrm{center} \times \sin(\theta) \\
u^* & = \cos(\theta) \times (u+t_x) + \sin(\theta)\times(u+t_y)\\
v^* & = -\sin(\theta) \times (u+t_x) + \cos(\theta)\times(u+t_y)\\
\end{split}

The formula doesn’t look correct to me.

Why would someone apply $\mathrm{sqrt}()$ and subsequently apply a power of $2n$ as it makes the $\mathrm{sqrt}()$ operation redundant?
$H(u,v)$ transfer function looks like a high-pass filter rather than a band-pass filter.
Shouldn’t $\mathrm{center}$ be $\mathrm{center}_x$, and $\mathrm{center}_y?$
Shouldn’t $(u+t_y)$ be $(v+t_y)$?

I agree with Peter K.  It looks like the authors of the Multidirectional Scratch Detection and Restoration in Digitized Old Images (E. Ardizzone, H. Dindo, and G. Mazzola, EURASIP Journal on Image and Video Processing, 2010) paper tried to reformulate the equation (shown by Peter K.) from the Analyzing Oriented Patterns paper by Kass and Witkin (1987; Computer Vision, Graphics and Image Processing, Vol. 37. No 3, pages 362 to 385).

Comparing the equations posted by the OP with the equations of a Gabor Filter, it feels like the author's are trying to make a novel contribution with their paper by reformulating the work of Kass and Witkin.  Although neither paper includes references to Gabor filters directly, Ardizzone et al. justify the use of their reformulated "steerable" filter by stating,

  The filter order n can control the slope of the subband, so that we chose n=4 to concentrate the filtering in a precise zone of the spectrum.  This explains the use of this bandpass filter instead of a Gabor one.

My experience with image filters is limited, but I believe a Gabor filter can realize the same apparent functionality as the equation given by the OP including bandwidth selection (see this entry, for example).

I think "subband" refers to one of the directionally dependent filters and the "center" is the origin of the Gabor kernel relative to the image on which the kernel is being applied.

I hope this helps.


Impulse and unit step response I'm trying to understand impulse and unit step responses, but I'm not sure if my understanding is correct. In general, are the following statements true:

a) The impulse response for a system $x(t) \mapsto y(t)$ is $y(t)$ when $x(t) = \delta(t)$ where $\delta$ is the Dirac delta, and

b) The unit step response for a system $x(t) \mapsto y(t)$ is $y(t)$ when $x(t) = u(t)$ where $u$ is the unit step (Heaviside step function)

In response to MBaz's comment, if we define the impulse response (edit) for Linear Time Invariant systems as

  "the function $h(t)$ such that the mapping $x(t) \mapsto y(t)$ is computed as $y(t) = x(t) * h(t)$"

then notice that if $x(t) = \delta(t)$ you get $y(t) = \delta(t) * h(t) = h(t)$.

Thus louism's statement a) and MBaz's comment agree.


Psychoacoustic Masking Curve I am trying to implement a Weighted Matching Pursuit algorithm. Essentially I would like to weight the atoms of a Dictionary (redundant matrix) with a frequency masking function, like the one specified for MP3. I have read some of the papers, and theory of frequency masking, but I am still having trouble figuring out how to generate the function in Matlab. My dictionary contains frequencies ranging from -fs/2 to fs/2. fs = sampling frequency.

If anyone can help me out with a quick way to generate this kind of function I would be most grateful.

Psycho acoustic masking is a VERY complicated subject. The masking thresholds (energy required to have one sound be audible in the presence of another) is a very complex function of

Spectral content and overlap
Spatial presentation (there are many books just on binaural masking alone)
Time shapes of the waveforms
Time alignment and envelopes of the wave forms

Even the more advanced audio codecs (such HE AAC) use very simplified models because a more accurate one is very complicated.

The best choice really depends on your application and requirements, monoaural masking of simultaneous sinusoids is probably the easiest case. You may find some helpful stuff here https://engineering.purdue.edu/~malcolm/interval/1998-010/


Window length in matplotlib's specgram In the paper about speaker verification I can read that

Two quantities have
to be set: the length of the window and the shift between two
consecutive windows. For the length of the window, two values are most often used: 20 milliseconds and 30 milliseconds. [...] For the delay, the
value is chosen in order to have an overlap between two consecutive windows; 10 milliseconds is very often used.

and later

Once an FFT algorithm has been chosen, the only parameter to fix for the FFT calculation is the number of points for
the calculation itself. This number N is usually a power of 2
which is greater than the number of points in the window,
classically 512.

However, in the matplotlib's specgram I can only see NFFT and noverlap parameters which, according to the documentation, correspond to the second quotation.

How can I then set the window's length or what is the default value of it?

You'd generally do the FFT without zero padding, i.e. window size == FFT length.

Having an FFT larger than the window leads to an interpolation, but doesn't add any information to the signal, and that's only useful in certain cases – maybe the paper presents one such case.


Dynamic Time Warping of inverted binary signals I am trying to find a match between a reference binary signal and a large set of data. From experience, the reference signal can appear time-shifted(that is why I've chosen dtw) and inverted, meaning a reference [101] may be represented as [010]. Is there a clever metric that i can use, or some alterations of the signals (some kind of norming before dtw), such that, in a single pass, i can check both normal and inverted bit signals? The obvious and labourish way would be to analyze the bit signals, revert them, and make a second dtw pass, but i'm trying to beat that.

I had a similar experience, i had to align side channel traces, I've tried to implement my own DTW that do lot of things not only aligning, but the this method was extremely slower than separating the processing and DTW.

So if you are using python the FastDTW will do the job insanely quicker than any handmade DTW, so the laborious way is the quickest way to get the aligned signals, besides that on very large signals DTW need huge amount of memory, so separating the two process is in my experience is the best solution.


What does "invariant to bias and gain changes" mean? Reading this paper, about the ECC criterion (which is included in opencv3), In section 3, the following is stated, with regard to the following equation:

(4) $\left\|\frac{i_r}{\left\|i_r\right\|} - \frac{i_w(p)}{\left\|i_w(p)\right\|}\right\|^2$

  It is apparent from (4) that our criterion is invariant to bias and gain changes. This also suggests that our measure is going to be invariant to any photometric distortions in brightness and/or in contrast.

What does "bias", "gain" invariant mean?

Thanks

You need to be careful when reproducing formulas! Your $i_r$, $i_w$ are actually $\bar{i_r}$, $\bar{i_w}$ in the paper and two sentences above $(4)$, it says these are

  their zero-mean versions, which are
  obtained by subtracting from each vector its corresponding
  arithmetic mean.

So, the symbols in your formula are always the same, regardless of whether there is a constant value added to them (that's the definition of bias), because that would automatically shift the arithmetic mean by the same amount.

The division by the individual vector's norm also "undos" any multiplication with a scalar (that's the definition of gain).

So, the overall formula is unchanged (that's what invariant means) when you add a bias to the vectors, and unchanged when you apply gain to the vectors.


Is there an example of an infinite signal? In the real world, Is there an example of an infinite signal?

After the proposed heat death of the universe, no active life form will be available to measure the rest of any given signal. That's not a real world to anyone.

One could define a signal not as function of time but as function of location, and wrap around the location. An example is the present ambient temperature as function of distance traveled to the east from Pontianak, Indonesia, on the Equator. We don't know if the universe is infinitely large so measurements on a straight line would be a bit iffy.


Deep space communications BER and FEC? What kind of bit error rate do they get from deep space communications (Pioneer, Voyager, et.al.), and what kind of modulation and FEC allows them to recover messages with that microscopic level of received signal power?

Are there more modern modulation methods and encoding schemes for similar channel conditions?

For many years the state of the art was to use a convolutional "inner code" and a block "outer code".  The "inner" and "outer" terminology come from the following block diagram:

$$\boxed{\scriptstyle \rm Payload}{\longrightarrow}\boxed{\scriptstyle\textrm{Outer Encode}}{\longrightarrow}\boxed{\scriptstyle\textrm{Inner Encode}}{\longrightarrow}\boxed{\scriptstyle\rm Channel}{\longrightarrow}\boxed{\scriptstyle\textrm{Inner Decode}}{\longrightarrow}\boxed{\scriptstyle\textrm{Outer Decode}}{\longrightarrow}\boxed{\scriptstyle\rm Payload}$$

Convolutional codes were used as the inner code because they are very powerful and can correct a large number of bit errors.  They have a weakness, though- when there are lots of errors that are close together they can break down and spit out errors in a burst at that location.  The outer code is used to correct those bursts of errors.  Block codes are not as powerful as convolutional codes (do not use as many parity bits/symbols either), but they are good at dealing with bursts of errors.  Also, there was usually a deinterleaver in between the inner and outer codes that spread the bursts of errors among many blocks, making it even easier for the block code to correct them.

As Wikipedia's Deep Space Telecommunications section says, early on the inner/outer codes were Viterbi (convolutional) and Reed-Muller codes.  Later they were Viterbi and Reed-Solomon codes.

In the early 90's Turbo codes were discovered and took the FEC world by storm.  In the 2000's low-density parity check codes have grown in popularity.  They were discovered in 1960 by Gallagher, but were not feasible to implement until recently because of the computational load that they require.  Both Turbo and LDPC codes are near optimal in the sense that they get very close to the Shannon limit of what it is possible to achieve with FEC.  Presently NASA uses both Turbo and LDPC codes, as far as I'm aware.

Like designing any reliable communications system, designing reliable deep space communications requires more than just adding powerful FEC.  Signal power, free-space path loss, receiver noise, etc., must be taken into consideration.  Deep space communications actually have a lot of advantages and two enormous disadvantages.  The disadvantages are the enormous distance and the limited transmitter power.  The advantages are the really high-gain directional antennas, the low noise that the earth dishes get from looking into empty space, the even lower noise they get by cooling their receivers with liquid nitrogen, etc.  They can also slow down their data rate while keeping the transmitted power constant to give each bit more energy.


How to implement sinc interpolation I'm trying to write my own high quality audio sample rate converter. I barely know anything about signal processing though so I need help.
From what I understand I need to sum together normalized sinc functions that touch every sample to find the value at any arbitrary point. I'm guessing this is hard because that would mean more than a billion calculations for a few second long audio file.
So do I only use the sinc functions for the samples that are closest to my x value? Is this what is meant by "windowed sinc"? How many samples should I go in each direction away from my x value? Additionally, my DAW has something called "32 point sinc" resampling. What is the "32 point" supposed to mean in this case?

"So do I only use the sinc functions for the samples that are closest to my x value?"
Yes, when you are truncating.
"Is this what is meant by "windowed sinc"?"
Yes.  The sinc goes to infinity, calculating that is impractical.
"How many samples should I go in each direction away from my x value?"
The sinc diminishes the further away you are.  At some point additional points become insignificant (no wider than a whisker on a gnat).
"Additionally, my DAW has something called "32 point sinc" resampling. What is the "32 point" supposed to mean in this case?"
I'm pretty sure it means a 32 point window, centered on your current sample.
Sinc interpolation is one form of interpolation.  It is the "Fourier compatible" one as it gives limits on a bandwidth basis.
Here is another link that may be even more useful:  Multi-channel audio upsampling interpolation
On re-reading this reference, Olli has already done this analysis, spectacularly.


Multithresholding- OTSU method I'm trying to segment a specific part in a gray-scale image and I want to achieve this by applying the OSTU algorithm, but after using this function:

im = imread('ready.jpg');
oim =otsu(im,10);
figure,imshow(oim, [1 10]);

I could show the output image, but I do not know how can I save the result to be able to use it later as a mask.

I am sorry that I cannot see where your problem is. What do you mean by "how can I save the result?" If you want to use $oim$ in matlab, then it is already there. If you are interested in extracting one specific class, you can simply get the corresponding $n$th mask by calling $mask=(oim==n)$, which returns a binary mask with only the pixel belonging to the $n$th class marked as 1s.

If you are going to save the classification result outside of matlab, for example a tif image, simply use imwrite function in matlab. Late you can read it via imread function.


Turn acausal Butterworth output into output of a causal Butterworth without input signal I have a dataset which contains a large number of time-series which have been filtered with acausal Butterworth filters. For a real-time application, I can only use causally filtered data.
Is it possible to turn the output of an acausal Butterworth filter into what would have been obtained if the original time-series had been run through a causal Butterworth filter?
(I do not have the original, unfiltered, time-series)

I do know the order of the filter that was applied, the corner freuquencies, and some parameter 'nroll' of which I don't know what it is ... (nroll is usually equal 2.5)

This is possible, as long as you know the filter which was applied to the data in the first place.

Assuming that the original acausal filter is $H_a(\omega)$ and your filtered data is $Y_a(\omega)$, then the original data $X(\omega)$ can be retrieved by inversely applying the filter to the data:
\begin{equation}
X(\omega) = \frac{Y_a(\omega)}{H_a(\omega)}
\end{equation}

This is simply a re-arrangement of the equation which can be used to describe the initial  filtering operation:
\begin{equation}
Y_a(\omega) = X(\omega) H_a(\omega)
\end{equation}

I've illustrated the concept using frequency domain notation, but the concept is the same for the time domain.  Depending upon your acausal filter, you may encounter some problems inverting it, I don't think this will be an issue in this case.


Measuring local "flatness" I have a function which looks like so:

I would like to know if there is a measure for a "weak type run length".
How can I know that a function is approximately constant in some specific neighborhood?
How can I differentiate this function from random noise?

You can run the signal through a low pass filter to get rid of the high frequency content. After that, if you know the expected length of the flatness of the signal, you can perform matched filtering with a template of that length. If the cross-correlation result exceeds some threshold, consider the signal detected.


High pass filter of 0.01 Hz for EEG slow drifts I am trying to make a high pass filter for electroencephalographic (EEG) datasets, in order to get rid of very slow drifts. However frequencies around 0.3 Hz are very important for the research these data serve.

I use Matlab and I tried various things. It seems that a very low dB attenuation, of 3 dB is fine with FIR filters. For example I used that:

h=fdesign.highpass('Fst,Fp,Ast,Ap',0.005,0.01,3,1,250); % i needed to cut
d=design(h,'butter');
fvtool(d)

However this filter makes a dc drift at the start of the data and I would not wish to loose these first datapoints.

I have read through your pages that a clever idea is to design a low pass filter that later I would subtract from the real data. I did it with an equiripple filter and it worked but left a baseline lift. I do not want to run now average removal as this is in a next step in the analysis protocol.

Any suggestions?

You're asking for a time machine: in order to avoid the DC bump at $t=0$, you need to know the state of the filter as if it had been running before the recording ever began. That's not going to happen no matter what kind of filter you use.

Here's a trick which may come in handy for you. Let's say you have an $N$-tap FIR filter (or an IIR filter whose impulse response decays adequately to 0 after $N$ samples). Take the first $N$ samples of your signal, reverse them, and prepend that to the start of your signal. In effect, we're defining a new signal $g(t)=g(-t)=f(t)$ for the purposes of filtering. When you filter that, the filter state will have already been "primed" with the initial signal: delete the first $N$ samples of the output, and the DC bump should have hopefully gone away.

Needless to say, there are numerous gotchas with that scheme — for starters, it's harder to do for realtime processing, and if the signal derivative is high at $t=0$, you are still going to have trust issues with the first few samples. But it should get rid of the bump.


12 bit ADC to 8 bit DAC A bit of background information:
I plan to start experimenting with digital signal processing on microcontrollers. The particular microcontroller I plan to use (TI Stellaris LM4F230H5QR, from the Stellaris Launchpad) has a 12 bit ADC integrated on it, while on the other hand, all DACs I that look like possible candidates for me to use are 8 bit.

Obviously, before getting the output, I'll have to reduce the bit depth somehow. I could just get an 8 bit ADC and avoid the issue altogether, but I'd like to keep as much of the signal processing line at 12 bits as I can.

Questions:
First, am I right in my assumption that by keeping 12 bits for as long as I can, I'll get better results?
Second, is there a smarter way of turning those 12 bits into 8 bits than just shifting everything to the right four times?

The first thing I would enjoin you to do is to measure the actual accuracy of your ADC. The 12-bit performance announced by the manufacturer can often be reached only in ideal conditions, such as:

Using an external precision voltage reference rather than the built-in one.
Using a distinct power supply for the analog section of the MCU.
Using all the necessary EMI isolation/filtering parts as specified in the app notes... that your target dev. board might not have because of cost cutting (it can be as simple as a ferrite between VCC and AVCC).
Using slow sample rates.
Preventing any I/O pin to switch during acquisition (eg, disabling a SPI or UART spitting data in the background, disabling interrupts causing a GPIO to be toggled, disabling PWM peripherals..), or even putting the entire chip in sleep mode during acquisition.

In non-ideal conditions, I have seen ADCs lose up to 3 bit of precision. All this to say that your assumption that doing the processing with 12-bit data all along is correct as long as you truly have a resolution of 12-bit on your input; but if it's actually 8 bits because you are far away form the "sweet spot" of the ADC, there's no point keeping and processing 4 bits of noise all the way - you can just do the reduction at the source.

Converting from 12 bits to 8 bits through rounding or shifting is going to cause quantization noise correlated with your input signal. Dithering algorithms can be used to solve this problem and spread the quantization error. You also have the option of shifting the quantization noise in a frequency band where it is less harmful, through noise shaping methods.


What does it mean that one frequency aliases another? So I was watching a video on sampling and the professor asked which frequencies will alias to 0 and which to $B$ if we are sampling at frequency $f_{s1}$. What does mean that two frequencies alias each other ?

The answer is that $0, f_{s1},2f_{s1},3f_{s1},\ldots, nf_{s1}$ will alias to 0 and $f_{s1}+B, 2f_{s1}+B, \ldots, nf_{s1} + B$ will alias to $B$.

Could anyone explain this please? The image shows the frequency spectrum in question.

I assume that since you're watching lectures on sampling theory, then you have already grasped the core ideas of basic signal operations and Fourier transforms and signal spectrum plots.

Then, consider a signal $x(t)$ whose C.T. Fourier transform is $X(\omega)$. When this signal is multiplied by the (ideal) periodic impulse train $$ \sum_{k=-\infty}^{\infty} \delta(t-kT_s) $$ where $T_s$ is the sampling period (and therefore $F_s = 1/T_s$ is the sampling frequency in Hz.) the resulting signal is called as the sampled signal $$ s(t)=x(t) \sum_{k=-\infty}^{\infty} \delta(t-k T_s) $$

Then it can be shown (by applying properties of CTFT) that the spectrum $S(\omega)$ of the sampled signal $s(t)$ is the periodical extension of the spectrum $X(\omega)$ of the signal $x(t)$, i.e., $S(\omega)$ is obtained by shifting $X(\omega)$ to left and right in frequency by the amount $ k \times \omega_{s} = k \times 2 \pi F_s$.  ( Here I'm using radian frequency $\omega$ which is related to Hertz frequency $f$ by $\omega=2\pi f$ ) :
$$ S(\omega) = \frac{1}{T_s} \sum \limits_{k=-\infty}^{\infty} X(\omega -k 2\pi F_s) $$

For example consider the the shift for $k=-1$, then $X(\omega+2\pi F_s)$ is added on top of original spectrum $X(\omega)$. Now by setting $\omega=0$ in $S(\omega)$, it can be seen that $X(0+2\pi F_s)=X(2\pi Fs)$ has come to the origin in addition to $X(0)$. Setting $k=-2$ would show also that $X(4\pi F_s)$ has come to the origin. This way it can be deduced that those frequencies given by $\omega_k = k 2\pi F_s$ are all appearing at the origin in $S(\omega)$. This appearance of those frequencies $2 \pi k F_s$ at the origin (or more generally any frequency  $\omega_k = 2 \pi k F_s + w_0$ at the frequency $\omega_0$) is called as aliasing in the spectrum of $S(\omega)$


fourier transform to the power of two Could someone explain to me why the computation speed for the fast Fourier transform increases by padding the series with zeros to the point that its length is close to a power of 2? This is common in the matlab environment, for example:

http://www.mathworks.co.uk/help/matlab/ref/fft.html

Firstly the terminology: The thing that you are trying to compute is the Discrete Fourier Series (DFS). This the only flavor for Fourier Transforms that is discrete in both time and frequency and can thus be represented numerically inside a computer. The Fast Fourier Transform is a specific algorithm to computer the DFS (see for example http://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm ). The term FFT is often used to refer to the DFS but that's actually somewhat wrong and often confusing.

You can always calculate the DFS directly using it's definition but this will require $N^{2}$ complex multiplies. The FFT speeds this up by breaking splitting the vector N into smaller sub vectors and then using the symmetry properties of the transform coefficient (often referred to as "twiddle factors"). The vector needs to be split in equal pieces though and that works better the more prime factors N has. A power of two is the best since it has the most prime factors and the factors themselves are the smallest possible. Worst case are N's that are prime themselves.

Let's assume you want to calculate the DFS of sequence of 1021 points. Since it's a prime number you need to use the direct formula which requires 1021*1021 multiples, roughly a million. If you zero pad to 1024, you can use the most efficient FFT version and the number of complex multiplies is 2*log2(N)*N is only about 20000. Hence that's a lot faster.

FFT is efficient if your DFS length can be broken down into lots of small prime numbers. Large prime numbers are slower. In practice you can typically find a spot made our of 2s and 3s that's fairly close to where you want to be so there is no need to always go to a power of 2.


Why is mutual information symmetric but conditional entropy isn't? Conditional entropy $H(X|Y)$ tells us how much the average uncertainty about a channel input $X$ is after observing channel output $Y$, and mutual information $I(X,Y)$ measures how much information about channel input $X$ can be obtained observing channel output $Y$.

Doesn't this mean that $H(X|Y)$ and $I(X,Y)$ complement each other or tell us the same thing? If so, then since mutual information is symmetric $I(X,Y)=I(Y,X)$, shouldn't conditional entropy be symmetric as well?

No,because while mutual information tells you how much uncertainty can a random variable $X$ remove from another random variable $Y$, the conditional entropy of $X$ given $Y$ tells you how much uncertainty remained in X after using the information (given by $I(X;Y)$) that $Y$ gave it.

$H(X|Y)=H(X)-I(X;Y)$

That's why the conditional entropy depends on the value of the entropy before the observation and the mutual information isn't, because it is only the difference ($\delta$) between two entropy states, before and after the observation.


Higher order spline interpolation I noticed that spline interpolation with a degree higher than 3 (everything beyond cubic splines) have a very high interpolation error, hence the prediction is mostly awful. I've come across various lecture notes, slides and Youtube videos that simply indicate that cubic splines (3rd degree) are optimal and that anything beyond that is a bad idea. These sources however never mention why this is the case.

Can anyone explain to me why this is the case and maybe give me a title/link to a journal/conference paper that explains this or maybe even gives a proof.

There is no such proof because it's not always true.  It's a rule of thumb, because I guarantee that you could come up with a situation- an infinite number of situations actually- where higher order splines would do better than cubic splines.  The optimal spline order for any given situation is the exact same order as the system you are trying to model.  If the order is the same and your data points are error free (never the case, of course, except for in theoretical problems), then you should be able to model the system perfectly.

The reason they recommend not to go higher than cubic splines is because overfitting is really, really bad.  Overfitting can greatly magnify errors, while "underfitting" (choosing a spline method with an order lower than the order of the system you are modelling) introduces some low pass filtering that is either not that bad or sometimes even beneficial.


How can I translate what I have learnt in Transfer Functions to differential equations? I have gone through the entirety of K. Ogata's Modern Control Engineering, and I do not understand how can I translate all the transfer function models into differential equations? For example, how would a closed loop control system with PID compensation look like as differntial equations??

You assign variables for each node that branches off or where signals join (eg the feedback adder). Then you write down the transfer functions between those variables. And as a last step you convert the functions back from Laplace domain into the time domain. Et voilà, you have differential equations.

You actually don't need that many variables. Only input, output and one per integrator and per differentiator is required. But starting off with more variables makes the functions simpler and thus the back transformation. At the end you can always get rid of the excess variables.


how to detect short pause in a speech with noise? I would like to count the number of short pauses in a speech and know the duration of the pauses by loading the signal into Matlab but I am very new to Matlab.
Can anyone help and give me some ideas about this?

You can use Voice Activity Detector(VAD) for detecting the pauses in speech and there location(and duration). if your input signal is not very noisy and noise is not varying much you can use a fixed threshold on the energy(calculated for short frames, ex-20 ms), so if energy is above that threshold you declare that frame as speech else pause. if input signal is very noisy(varying noise) you can use adaptive threshold and do the same,i.e. you can have short term(lasting for 10-12 frames) and long term energies(lasting 40-50frames) and long term energy can be used as threshold. once this classification is done for entire input signal you can get duration's of pauses in speech


How to create a swept / chirp signal between a given starting audio signal and a given ending audio signal I can create a swept / chirp signal using matlab / octave with the code below (also see spectrum taken).  But how can I create a swept / chirp signal between a given starting audio signal and a given ending audio signal?  (One way to look at it is morphing one image into another image but using two audio signals)  or a sliding / chirp phase vocoder

clc,clear all
swept_startfreq = 500;
swept_endfreq = 225;
fs = 8000;
swept_dursec = 15;  % duration of signal in seconds
swept_t = (0:swept_dursec*fs)/fs;  % Time vector

if (swept_startfreq == swept_endfreq)
      swept_phase = swept_startfreq * swept_t;
else
    swept_alpha = abs(log(swept_endfreq/swept_startfreq))/swept_dursec;
    swept_b = (swept_startfreq - swept_endfreq)/(1 - exp(-swept_alpha*swept_dursec));
    swept_a = swept_startfreq - swept_b;

    swept_phase = swept_a*swept_t - (swept_b/swept_alpha)*exp(-swept_alpha*swept_t);
endif

swept_sig_out = .8* cos(2*pi*swept_phase); %equation used to create signal
wavwrite([swept_sig_out'] ,fs,32,strcat('/tmp/swept_sig_out.wav'));  % export file

The frequency spectrum of the original code / chirp or swept signal it creates

So I would have an initial audio signal (use as the initial audio signal for testing) and an ending audio signal (use as the ending audio signal for testing).  How can I do a "morphing" between these two signals?

Here's what an example of the final file would sound like example of final audio signal when processed (please not this isn't exact but it shows how the pitch gradually changes from the original signal to the end audio signal over time  If you would like to know how I created the final file I used the initial signal and Audacity with the effect - Sliding Time Scale/Pitch Shift to get as close as possible to what the ending signal sounds like.

Added from comments: Start with the frequency with the largest amplitude in the initial signal and "morph" it towards the frequency with the largest amplitude found in the ending signal

It would sound like "one one one one one" gradually decreasing or increase in pitch depending on the ending audio signals largest amplitude frequency found.   Please note the Initial and ending audio signal may not be the same type of signal, the initial signal could be "one one one one" and the ending signal could be "two two two two"

Ps:  I'm using Octave 4.0 which is like matlab

As noted by endolith in the first comment, the problem is not well defined, in that, there's no one right way of doing this audio morph/crossfade.

One method would be to use some form of interpolation between the source and destination spectra as discussed in this paper and implemented in this toolbox. A similar technique is also discussed in this paper.

This paper proposes another possible method inspired by the problem of heat diffusion. Imagine the left edge of a metal plate is held at a "temperature profile" that resembles the magnitude spectrum of the first signal and the right edge is held at a "temperature profile" resembling the magnitude spectrum of the second signal. The 2D "temperature profile" function between the two edges is reconstructed by solving a heat diffusion problem. An intuitively pleasing property of crossfades obtained with this method is that it produces smooth chirps if two sinusoidal frequencies are close to each other, but produces a fade-in/fade-out effect if they are farther apart. See their audio examples here.

I have not succeeded in locating good Matlab implementations of audio morphing algorithms. But there are a couple of good starting points. For instance, the Java library in the Spectral Toolbox has an audio morphing Java source code which can be converted to Matlab with some effort. Another good starting point might be this Github repo.

W. Sethares, J. Bucklew, Kernel techniques for generalized audio crossfades, Cogent Mathematics 2(1), 2015.

W. Sethares et al., Spectral Tools for Dynamic Tonality and Audio Morphing, Computer Music Journal 33(2), 2009.

M. Slaney, et al, Automatic Audio Morphing, Proc. IEEE ICASSP, 1996.


How to edit audios so that they have the same length I have a dataset of audio files, which differ on their lengths by a couple of milliseconds. Some are little shorter than 1s and some are little longer than 1s. They all have the same sampling rate, namely 20kHz, what differs is the number of samples.

Does anyone know a way to automatically transform them into 1s audio files with 20kHz sampling rate?

EDIT:
I need to transform the audios to having the same length so that I can input them into a standard fully-connected neural network.

Cutting the audio in the beginning and at the end would be the best solution, I think, if the audio is longer than 1 second. For example, if the audio is longer 1.1s, then I would like to cut the audio by 5ms in the beginning and the end to transform it to a 1s audio? And if the data is shorter, somehow putting some "silence" at the end of the data?

Thank you in advance!

ImageMagick is a great tool to edit multiple images in the command line. Thanks to your question, I just discovered SoX, or Sound eXchange:

  SoX is a cross-platform (Windows, Linux, MacOS X, etc.) command line
  utility that can convert various formats of computer audio files in to
  other formats. It can also apply various effects to these sound files,
  and, as an added bonus, SoX can play and record audio files on most
  platforms.

It can trim, crop, cut multiple files (and resample too), see for instance:

Editing Multiple Mono Sound Files in SoX,
Audio format conversion cheat sheet (aka how to)

For shorter files, an option is to create a silence file, happen  it, and crop again, see for instance: Linux command to extend the duration of audio files.


Creating a convolution kernel from filter description I am not quite sure whether the question I am stating is right, or if I am trying to combine two not related topics but here it comes..

I am trying implement and apply a contra harmonic mean filter, but do not know how the kernel should look like, and why it should like that, or even if should do it with a convolutional kernel....

The contra harmonic mean filter is described here.

http://www.blackice.com/Help/Tools/Document%20Imaging%20SDK%20webhelp/WebHelp/Contra-Harmonic_Mean_Filter.htm

Sorry for being a noob within this area, but seems like this information isn't stated anywhere, and I kinda get the feeling that it is some form of knowledge everybody should have, but not me...

I am using filter2D to create my filter..

The key point in the description on the linked-to page is:

  This function will filter the image by the nonlinear contra-harmonic mean method.

that the method in question is nonlinear.

That means convolution kernels are not really applicable, because they are generally only useful for linear, time-invariant systems.

Having said that, you may be able to use a simple linear convolution kernel in the process of calculating this mean.

You need to calculate:
$$
\sum_{(i,j) \in M} A(x+i,y+j)^Q
$$

where $Q$ is either $P$ or $P+1$ and $M$ is an $N\times N$ square mask ($\square$).

One way to do this would be to:

Calculate the image to the power $P+1 = I_1 =   A(x,y)^{P+1}$.
Calculate the image to the power $P =  I_2 =  A(x,y)^{P}$.
Form
$$
I_\mbox{NUM} = \square * I_1\\
I_\mbox{DEN} = \square * I_2\\
\mbox{Contra-harmonic mean} = \frac{I_\mbox{NUM}}{I_\mbox{DEN}}
$$

Here, the convolution kernel, $\square$ is just:
$$
\square(x,y) = \left\{
\begin{array}{cl}
1 & |x|,|y| < N/2\\
0 & \mbox{otherwise}
\end{array}
\right .
$$
i.e. it's just a square $N \times N$ kernel of all ones.


How can spectrum of DTFT of sampled band-limited analog signal be in $[-\pi,\pi]$ Say the analog signal $x(t)$ and its spectrum $X_a(\Omega)$. After sampling with frequency $F_s$ or sampling period $T$ we get

$$x[n] \triangleq x(nT) = x\left(\frac{n}{F_s}\right)$$

and its spectrum $X_d(\omega)$.

The relationship between the two spectrums is

$$X_d(\omega) =\frac1T \sum ^{\infty}_{n=-\infty}X_a\left(\frac{\omega + 2\pi n}{T}\right)$$
$$\omega=T\Omega$$

For band-limited analog signal having spectrum in range $[-B,B]$. $X_a(\Omega)=0$ for $|\Omega| \geq B$ if $F_s$ is above Nyquist rate then
$$X_d(\omega) =\frac1T X_a(\Omega)$$

For digital signal, its frequency is restrict in $[-\pi,\pi]$ but not the case of analog signal. The equation above means that $X_d(\omega)$ is not restricted in $[-\pi,\pi]$

The DTFT $X_d(\omega)$ of a discrete-time signal is periodic with period $2\pi$, so it's sufficient to consider the interval $[-\pi,\pi]$. "Restricted" to this interval means that due to the inherent periodicity there can be no new frequency components outside the interval $[-\pi,\pi]$. Of course, $\pi$ corresponds to half the sampling frequency.

If $X_a(\Omega)$ is the spectrum of a band-limited continuous-time signal, and if the sampling frequency is greater than twice the band-width, then the equality

$$X_d(\omega)=\frac{1}{T}X_a\left(\frac{\omega}{T}\right)$$

holds in the interval $[-\pi,\pi]$. Outside that interval, $X_a(\Omega)$ is zero, and $X_d(\omega)$ is continued periodically.


Need help understanding convolution indexing Ok, this is leaving me increasingly frustrated. This builds on a previous question here but I am still not getting the indexing, although I learned from that question all the same.

Case 1: First convolution equation example:

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

The way I got the right hand side, was:

1) Make the $h$'s have dummy variables, so $h[n]$ and $h^*[-n]$.

2) Flip one of them around, so now I get $h[n]$ and $h^*[n]$

3) Add the delay $k$ into the one that way flipped, so finally I have $\sum_{n=-\infty}^{\infty} h[n] \ h^*[n+k]$.

Great, I have a convolution.

Case 2: Second convolution example: (Unrelated to first)

We have:

$$
\sum_{i=-\infty}^{\infty} h[i-m] \ R_{xx}[i].
$$

Here we want to try and find out what convolution this corresponds to. The way to solve this, (apparently), is to make the assumption that another variable $g[n]=h[-n]$, and call it a day. Then we get:

$$
\sum_{i=-\infty}^{\infty} h[i-m] \ R_{xx}[i] = \sum_{i=-\infty}^{\infty} g[-(i-m)] \ R_{xx}[i] = \sum_{i=-\infty}^{\infty} g[m-i] \ R_{xx}[i] = g[m] * R_{xx}[m] = h[-m] * R_{xx}[m]
$$

This works great.

Case3: The problem:

This is my problem. Case 1 showed the steps for how to do a normal convolution. Case 2 showed how the method of substitution works. However, if I try to solve Case 1 using this substitution method, I get the wrong answer, even though it worked on case 2. Why is that?

To wit - where am I going wrong here?

$$
h[k] * h^*[-k] = h[k] * g^*[k] = \sum_{n=-\infty}^{\infty} h[n] \ g^*[k-n] = \sum_{n=-\infty}^{\infty} h[n] \ h^*[-(k-n)] = \sum_{n=-\infty}^{\infty} h[n] \ h^*[n-k]
$$

However, this is **clearly* the wrong answer! What am I doing wrong here?

This question is plagued by bad notation that the OP is unwilling to abandon because, unfortunately, his book uses it. This greatly
impedes his understanding of the concepts.

  Case 1: First convolution equation example:
  $$
h[k] * h^*[-k] = \sum_{n=-\infty}^{\infty} h[n] \ h^*[n+k]
$$

  The way I got the right hand side, was:

  1) Make the $h$'s have dummy variables, so $h[n]$ and $h^*[-n]$.

  2) Flip one of them around, so now I get $h[n]$ and $h^*[n]$

  3) Add the delay $k$ into the one that way flipped, so finally I have $\sum_{n=-\infty}^{\infty} h[n] \ h^*[n+k]$.

${}{}{}{}{}$

No. 3) is not correct. When convolving, you flip one signal and then move it to
the right by $k$ which will give you $h^*[n-k]$, not $h^*[n+k]$.  Thus you end
up with $\sum_n h[n]h^*[n-k]$ as in the leftmost sum in $(2)$ below.

  Great, I have a convolution.  (emphasis added)

The right side of the last equation is not generally called
a convolution; it is called the autocorrelation function of $h$.
If $h$ and $g$ are two finite-energy sequences (meaning $\sum_n |g[n]|^2$
and
$\sum_n |h[n]|^2$ are finite), then the $k$-th term of
their convolution is
$$\sum_{n=-\infty}^\infty h[n]g[k-n] = \sum_{n=-\infty}^\infty h[k-n]g[n].
\tag{1}$$
It is not a good idea to use $h[k]\star g[k]$ to denote the
convolution in general or its $k$-th term in particular
-- far better to write
$h\star g$ in general and $(h\star g)[k]$ in particular -- as the
OP has already been told repeatedly on math.SE. Be that as it may,

  for the case when $g[n] = h^*[-n]$ for all $n$, we can write $(1)$ as
  $$\sum_{n=-\infty}^\infty h[n]h^*[n-k] = \sum_{n=-\infty}^\infty h[k-n]h^*[-n]
= \sum_{m=-\infty}^\infty h[m+k]h^*[m]
\tag{2}$$
  where the last equality was obtained by replacing $n$ by $-m$.
  Notice the difference from the OP's result. For complex-valued
  sequences, there are two different ways of defining the autocorrelation
  function $R_{hh}[\cdot]$ depending on which sequence we choose
  to conjugate, but both ways give the same conjugate symmetry
  property $R_{hh}[-n] = (R_{hh}[n])^*$.

${}{}{}$

  Case3: The problem:

  This is my problem. Case 1 showed the steps for how to do a normal convolution. (emphasis added)
  Case 2 showed how the method of substitution works. However, if I try to solve Case 1 using this substitution method, I get the wrong answer, even though it worked on case 2. Why is that?

  To wit - where am I going wrong here?

Where you are going wrong is that your analysis
of Case 1 is not quite correct
as explained above, while your Case 2 is not a convolution at all,
and is not an autocorrelation or cross-correlation of complex-valued
sequences either, because these require conjugation of one of the
terms.  So you are really mixing up things and then complaining that
the answers are not the same.


ZOH DAC amplitude formula testing We have an expression shown in red bellow for non ideal DAC ZOH,which links the sampling frequency to the amplitude loss of the input signal.
I took $f_b=1kHz$ and $f_{noise}=4kHz$
$f_{in}=f_b+f_{noise}.$
Our ADC is a sinc function with $x=10kHz$ so I did a convolution for out_put_ADC=conv(fin,sinc).
And basically that is the point i got stuck with the MATLAB code shown bellow.
If everything was correct till this point how to implement the DAC so we will witness the amplitude deterioration formula of ZOH as a function of sampling frequency?
Thanks.

x=0.0001:0.0001:1; %10KHz ADC sampling

ADC_response = sinc(x);

fb=1e3;  %1Khz frequency

f_noise=4e3;  %1Khz frequency

yb=sin(2*pi*fb*x); %1Khz frequency sine input signal

yf_noise=sin(2*pi*4*fb*x); %4Khz frequency sine input noise

y_in=yb+yf_noise;%total signal going into the ADC

ADC_out=conv(y_in,ADC_response);%signal after ADC

Here is the entire simulation code.
First, to verify the spectrum response of ZOH, I replaced the DAC_output with an impulse signal. So I can directly check if the simulation confirms the theory or not. ( As ZOH is a kind of LTI filter, the spectrum of impulse response represents the frequency response of this system. You can try other input signals like the one you originally provided.)

Input test signal:

Output Spectrum: Note that theory spectrum entirely cover the simulated spectrum.

Ps. DAC is simulated via an up sampling and a convolution with hold filter, h = ones(1, sample_hold)

A very high pseudo analog sample frequency should be given because there is no real 'analog' in computer simulation.

You can alter different Tp or Ts (not lower than analog pseudo freq.). These codes should work.

There are some rooms for improvement about the ADC process you provided in question. You can refer the my blog.

clear all;
close all;

end_time = 2e-2;
x=0.0001:0.0001:end_time; %10KHz ADC sampling

ADC_response = sinc(x); %% Please check the frequency response!!

fb=1e3;  %1Khz frequency

f_noise=4e3;  %1Khz frequency

yb=sin(2*pi*fb*x); %1Khz frequency sine input signal

yf_noise=sin(2*pi*4*fb*x); %4Khz frequency sine input noise

y_in=yb+yf_noise;%total signal going into the ADC

ADC_out=conv(y_in,ADC_response);%signal after ADC

%% Started here.

% ADC_out = ADC_out(1:length(y_in)); % trim for same length

ADC_out = zeros(1,length(y_in));
ADC_out(length(y_in)/2) = 1; % impulse

%% pseudo analog frequency & ZOH simulation
% parameters
analog_sampling_rate = 1e8; % 100 MHz.
digital_sampling_rate = 1e4; % according to question.

%% Important Time Parameters
Ts = 1/digital_sampling_rate; % 1e-4 as given.
Tp = 2e-5; % ZOH time in analog unit (sec), tunnable.
hold_sample = round(Tp*analog_sampling_rate);
%  This indicates how many samples need to be hold in this pesudo analog simulation.

%% DAC Processing
t = 0:1/analog_sampling_rate:end_time-1/analog_sampling_rate; % pseudo analog time.

% upsampling
y_up = upsample(ADC_out,analog_sampling_rate/digital_sampling_rate);

% zero order hold filter
h = ones(1,hold_sample); % zero order hold

% filter your signal
y_zoh = filter(h,1,y_up);
y_zoh(1:floor(mean(grpdelay(h))))=[];
y_zoh = [y_zoh zeros(1,floor(mean(grpdelay(h))))];

%% Time domain plot
figure;
plot(x,ADC_out);
hold on;
plot(t,y_zoh);
legend('yb','y_zoh');

%% Freq response of discrete time signal and freq unit conversion
d_N = length(ADC_out);
d_F = fft(ADC_out,d_N);
d_freq_step = digital_sampling_rate/d_N;
d_freq = -digital_sampling_rate/2:d_freq_step:digital_sampling_rate/2-d_freq_step;
d_freq = d_freq(d_N/2+1:end)';

%% Freq response of psudo analog time signal and freq unit conversion
c_N = length(y_zoh);
c_F = fft(y_zoh,c_N);
c_freq_step = analog_sampling_rate/c_N;
c_freq = -analog_sampling_rate/2:c_freq_step:analog_sampling_rate/2-c_freq_step;
c_freq = c_freq(c_N/2+1:end)';

%% Theory transform of DAC Transfer function.

c_thy_F = abs(sinc(c_freq*Tp)*Tp/Ts); %% For verification.

%% Plot Freq domain comparison.

figure;
plot(d_freq,abs(d_F(1:d_N/2))/d_N);
xlabel('Freq. (Hz)');

figure;
plot(c_freq,abs(c_F(1:c_N/2))*d_N/c_N);
hold on;
plot(c_freq,c_thy_F);
xlabel('Freq. (Hz)');
legend('simulated','theory');


reducing MIPS in detecting tones Currently our software is detecting tones using FFT and checking if there is only a single peak in frequency spectrum. Our software needs to know what kind of input it is before taking action. So we differentiate the case when tone is fed as input or if the input is music.

However we are trying to reduce the MIPS and wondering if there is any other better ways to detect tones.

I found correlation with previous data can be used but again don't know the right method to do it. Is there any better method?

Assuming that the tones are sufficiently sinusoidal, here's a possible scheme, with a bandpass filter being controlled by a frequency estimator:

If the tunable bandpass filter's output is at about the same power as the input signal, then the input consists mainly of a sinusoidal tone. It is perfectly fine if the frequency estimator is a simple one that only gives a meaningful output for the tone to be detected.