Plot spectrum of $X(-k)$ Assume FFT of $x[n]$ is $X(k)$ and I would like to find $X(-k)$. I have implemented two techniques:

I am calculating the X(k) = fftshift(fft(x)) and then find X(-k) to be fliplr(X(k)) ;
I write the for loop to correspond each symbol of in $X(k)$ to $X(N-K+1)$ where $N$ is the number of points used for FFT.

I get kind of similar answer but in the second approach I am mapping $X(N)$ to $X(1)$ which is DC, which should not be the case. Therefore I am suspecting I am doing something wrong.

Can anyone tell me what is the correct way to implement $X(-k)$ when you have the $X(k)$ ?

FFT (DFT) is defined as

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

Hence, $X[-k]$ is
$$X[-k]=\sum_{n=0}^{N-1}x[n]e^{jk2\pi n/N}$$

By investigating the formula of inverse DFT (and a quick change of variable $n\leftrightarrow k$ for more clarity),

  $$\mathsf{IDFT}\{X[k]\}=x[n]=\frac1N\sum_{k=0}^{N-1}X[k]e^{jk2\pi n/N}\Rightarrow x[k]=\frac1N\sum_{n=0}^{N-1}X[n]e^{jk2\pi n/N}$$

We can see that
$$X[-k]=N\cdot\mathsf{IDFT}\{x[n]\}$$
So instead of calculating the FFT and then converting the result, you can directly calculate $X[-k]$ by taking the IFFT.


Biquad Filter Non-Linearity I'm using a biquad filter algorithm in my project. I got 3, 16 size variable arrays summing into one and then gets filtered.

Now i want to filter each array seperately and sum after that, which if the filter was linear wouldn't be a problem, but result is non linear after all.

Am I doing something wrong with this process or biquad's non linear in nature?

If latter is true, can you recommend a nice filter to implement in embedded C, and musical sounding as well? (Quality factor is a must)

I'm using this to generate coefficients

http://shepazu.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html

Turns out I was using same state buffer and instance for all three filters.

Whenever an array finishes filtering, new filter uses the old state array for that different timeline and it results in non-linear outputs at the end/start of a new one.

I created different state arrays and instances for each different array and the problem got fixed. Working smoothly as intended now.


How to add additive noise to an image? I have the following image I and I want to add additive noise using the array N below. Can you please give me a calculation example of how to achieve this?

Do I have to move from image space to frequency space?

Thank you in advance!

You don't have to go to the frequency domain to add noise. I'm not completely sure what you mean by adding it, because those number in N are quite big. Are those variances or actually sample values?

Adding them together can simply be done by well... adding them together (value by value). Since they both have the same size, this is possible.

Normally one would give a noise distribution with its paramters. Those random numbers can then be generated and added to each element of the matrix.


Why averaging in Spectral magnitude domain not in complex domain to estimate spectrum of a process Consider we need the magnitude spectrum of the signal. Signal is recorded in $N$ trials to a certain stimuli. Signal varies from each trial because of noise. So, to estimate the actual magnitude spectrum, there are two approaches:

Avg. in temporal domain: Average out trials in temporal domain to obtain signal of less noise and then take its fourier transform to obtain signal's complex frequency content. Visualize magnitude spectrum. The output would be the same, even if we average the fourier transforms(complex), as it is a linear operator.
Avg. in magnitude spectrum domain: Compute fourier transforms of each trial, get magnitude spectrum and average this across trials.

According to me, approach (1) is correct as we are not neglecting any of it. Averaging has to be done in linear domain, which is real and imaginary parts of the complex output (Fourier transform). Whereas the second approach neglects phase spectrum of each trial, i.e. averaging in Polar representation of complex number individually. Many of my friends, even a professor approve the second approach (2). Their argument is that, because the requirement is the magnitude spectrum which is random, because of noise in each trial, so averaging the required output (Fourier transform) over many trials to obtain its estimate is a reasonable step. I am not convinced with it, as spectrum consists of both phase and magnitude and just neglecting one of it to get an estimate individually, by linear approaches, (here Arithmetic mean) is not proper, as the true underlying magnitude spectrum, getting modulated by noise, may be some form of non-linearity (function mapping temporal signal to magnitude spectrum). Thus, averaging magnitude spectrum out won't null out the noise. So, please help by providing opinions of what may be correct and why ?

You seem to understand the issues at hand.  In your situation, if the trial runs would be identical, except for the noise, averaging the signals (or the complex values) will definitely work, and I recommend that it how you do it.

The notion of averaging magnitudes (or magnitudes squared) across DFTs comes from when you are chunking a signal into multiple DFTs and want to aggregate the results.  In that case, if you average the raw DFTs you will get all sorts of cancellation/reinforcement effects depending on the frequency values of the signal compared to your DFT frame length.  When you average the magnitudes, the phase values are ignored and you will get more desirable results.

As is often the case wiht policy, it sometimes gets followed blindly without the understanding of the rationale behind the policy.  That may be what is behind some of the recommendations you have heard.

If you are trying to decompose your signal (assuming it is periodic) into its constituent harmonics, retaining the phase information is essential.

Hope this helps,

Ced


Relation between filter order and window length for FIR filter In case of FIR filter design, is the following relation always applicable or there are other relations/formulas also?
If window length is denoted by $M$ ,then filter order will be $M-1$

Yes. The number of taps or coefficients in a FIR filter (or window length as you call it) tell the order of the filter. A more general way to see it, to avoid confusions, is to think of the order as the largest delay in the filter. For example, consider the following naive reverb:
$$y[n] = a x[n] + (1 - a) x[n-D]$$
Even though there are only two (visible) coefficients, the order of the filter is $D$ because it’s the largest delay. The window length is $D+1$ but all the middle coefficients are $0$, hence not visible.


What is meant by 'combine the luminance and chrominance ' I am new with image processing. I'm trying to implement a research paper in MATLAB. I didn't understand a sentence.
Here it is:

  The luminance is enhanced by applying the proposed algorithm and it is
  combined with the chrominance and converted back to color image.

How can I "combine the luminance and the chrominance"?

Paper reference: Image enhancement with the application of local and global enhancement methods for dark images, 2017 International Conference on Innovations in Electronics, Signal Processing and Communication (IESC)

Intensity of color images are defined through color spaces such as RGB, YUV, and HSV.

RGB is an additive-synthesis primary color decomposition and it's very fundamental.

Whereas YUV is an example of Luminance (brightness or Black & White) and Chroma (color information saturation and hue) decomposition. Usually chroma consists of two components such as U & V in YUV or H & S as in HSV color spaces.

When luminance and chroma are combined together they create the complete definition, YUV or HSV, of the color information which can be converted into RGB. How they are combined depends on the chosen color space and its conversion into RGB.

For the case of YUV to RGB conversion, the combination and conversion is as simple as a matrix multiplication.


Which parameters should I take care of while converting .au to .wav for Genre Classification I am learning to develop Music Genre Classification models. I have a dataset in which files are in ".au"  format, which I need to convert into ".wav". I used SoX (Sound eXchange) to convert it. With a simple:

  sox input.au input.wav

the files are converted smoothly. Which parameters can be lost doing this, which can affect the accuracy of my model?

From a raw "music" point of view, both .au and .wav file formats may introduce a form of lossy compression, as compared to the basic Pulse Code Modulation (PCM). Partly because some parties to drift away for the standard standard, including mu-law companded versions, there are some inconsistencies in the AU and WAV formats.

As the .au is lossy compressed, what is important is that you do not lose even more information while converting to .wav. From your discussion with @Olli Niemitalo, the information seems preserved. However, I am a little worried that the file sizes are the same (.au file) (I tend to think the .au ought to be smaller, with the basic settings):

  The Au file format [...] Originally [it] was headerless, being simply 8-bit µ-law-encoded
  data at an 8000 Hz sample rate. Hardware from other vendors often used
  sample rates as high as 8192 Hz, often integer factors of video clock
  signals. Newer files have a header that consists of six unsigned
  32-bit words, an optional information chunk and then the data (in big
  endian format).

So you can check the header information to be sure of how the .au was compressed, what channels are here, etc. As long as this information remains in the .wav file, the chances of affecting genre classification are very mild.


Response of cascading FIR filters I can't seem to get the right answer for this. Could someone show me the procedure?

As shown in the figure a system is constructed by cascading four FIR filters having the following frequency responses,

If the input signal is $x[n] = 0.9n(u[n] - u[n - 20])$, compute the value of $y[10]$.

[Edit] working thus far

H[z] = (1 - e-jω)x(1 + e-jω)x(1 + e-2jω)x(1 + e-4jω)

 = (1 - e-2jω)x(1 + e-2jω)x(1 + e-4jω)
 = (1 - e-4jω)x(1 + e-4jω)
 = (1 - e-8jω)

y[n] = x[n] - x[n-20]

From this point i'm not sure how to proceed.

In this problem there are a cascade of 4 LTI systems whose multiplication of frequency responses (or convolution of impulse responses equivalently) provides the overall system Frequency Response (or impulse response equivalently) and the output $y[n]$ at n = 10 is being asked for an input of $x[n] = 0.9 \times n \times (u[n] - u[n-20])$

Now I would like to proceed by finding the overall impulse response $h[n]$ from the serial convolution of those four individual impulse responses as it seems simpler to find the output in time domain then in the Frequency domain:
$$ h[n] = (\delta[n] - \delta[n-1]) \star (\delta[n]+\delta[n-1]) \star
 (\delta[n] + \delta[n-2]) \star (\delta[n] + \delta[n-4])$$

Where $\star$ denotes convolution operation and those impulses are directly derived from the Frequency Responses $H_i(e^{j\omega})$ of the individual systems. I hope you can see the simple corresponce with each one of them.

Then convolving from left to right, we can show in 3 steps that:
$$h[n] = (\delta[n] - \delta[n-2]) \star (\delta[n] + \delta[n-2]) \star (\delta[n] + \delta[n-4])$$

$$h[n] = (\delta[n] - \delta[n-4]) \star (\delta[n] + \delta[n-4])$$

$$h[n] = (\delta[n] - \delta[n-8])$$

I hope also that you have no difficulty in obtainig shifted impulses as a result of convolution with a shifted impulse.

Therefore now the question becomes; given LTI impulse response $h[n] = \delta[n] - \delta[n-8]$ , compute $y[n]$ for the given input $x[n]$. It's now easy to show that:
$$y[n] = h[n] \star x[n] = (\delta[n] - \delta[n-8])*x[n]$$
$$y[n] = x[n] -x[n-8]$$

therefore with $x[n] = 9 \times n \times(u[n]-u[n-20]$ , $y[10]$ becomes:
$$y[10] = x[10] - x[2]$$
$$y[10] = 0.9 \times 10 \times (u[10]-u[-10]) - 0.9 \times 2 \times (u[2]-u[-18])$$
$$y[10] = 9 - 1.8 = 7.2$$

Where before the last line I have used the famous property of unit step $u[n]$.


How do you "split apart" the output of a linear estimator? Let's say I have a true signal $x$ which is corrupted by noise $n$ such that:
$y = x + n$
I am attempting to use a Weiner filter to estimate the true signal $x$. This estimator is of the form:
$\hat{x} = Gy$
Now let's say I keep the original equation and variables and add a known input bias $c$
$y + c = x + c + n$
The Weiner filter estimator would now take on the form
$\widehat{x+c} = G(y + c)$
The issue with this form is that I am now estimating $x + c$ instead of just $x$. Assuming I know the input bias value $c$, is it "okay" to simply break up the estimated output like the following?
$\widehat{x} = G(y + c) - c$
My confusion primarily comes about breaking up the output of an estimator. I assume because the Weiner filter is linear, this is okay.

There is nothing wrong about subtracting the bias. For the Weiner filter, the operation is a matrix multiply, there are no non-linear operations taking place which is reason everything works out. There are two possible ways to do it:

Subtract the bias before the filter. In this case you'll have a Weiner filter, $G_1$, which will produce $\hat{x}$ given $x+n$.

\begin{align}
\hat{x}&=G_1(y-c)
\end{align}

Subtract the bias after the filter. In this case you'll have a Weiner filter, $G_2$, which will produce $\widehat{x+c}$ given $x+c+n$ (this is the case you presented).

\begin{align}
\widehat{x+c}&=G_2(y+c)\\
\hat{x}&=\widehat{x+c}-c
\end{align}
Note: Be careful with the notation of $G$. The Weiner filter is constructed using the data so the filter constructed using the $x+n$ data will be different from the filter constructed using the $x+c+n$ data.


What interesting audio can be synthesized with ifft I'm working on a presentation on audio synthesis for some fellow software engineers. I want to include some "fourier synthesis" creating time domain samples from frequency domain samples using ifft.

Does anyone know of interesting sounds that can be easily synthesized with ifft?  I can do (bandlimited) square, saw, and triangle waves already, but i'm thinking more along the lines of voice or stringed instruments or something.  Does anyone know of anything that would be a simple to understand but interesting demo for fourier synthesis?

Thanks!

You can take a look at Photosounder. The idea behind is simple: Take an image as if it were an spectrogram and compute its IFFT, the result is  always interesting at least! It tends to work better with abstract images and geometric shapes, a great example is the video on the upper-right corner, "The sound of fractals".

I would never try to replicate a real instrument using this method, but you can encourage your fellows to try it as a challenge...


Impulse response for an LTI system I'm new to signal processing and working my way through a textbook. There is an exercise where a causal LTI system is given that responds to a rectangular pulse.
I have an exercise where a causal LTI system responds for a rectangular pulse
$$ s(t) = \Pi(t) = \operatorname{rect}(t) $$
with a triangular pulse.
$$ g(t) = \Lambda(2t) = \operatorname{tri}(2t)$$
One task is to find out what the impulse response for this system is. The solution is
$$h(t) = \frac{d}{dt} h_{\epsilon}(t) $$
because $$\delta(t) = \frac{d}{dt} \epsilon(t)$$
where $\epsilon(t)$ is the unit step function, $\delta(t)$ the unit impulse function.
and looks like this

I don't really understand the solution and how it turned to an endlessly oscillating rect pulse. Unfortunately, the textbook doesn't explain its approach. So, why does the output signal looks like this?
EDIT:

Translation:
A causal LTI system responds to a rectangular pulse  $s(t) = \operatorname{rect}(t)$ with $g(t) = \Lambda(2t)$
a) What is the response for $s_1(t) = \operatorname{rect}(\frac{t-1}{2})$?
b) What is the unit step response?
c) What is the impulse response?
And the solutions for a) - c) from the textbook

EDIT : Based on RBJ's comment that the indicated the system cannot possibly be an LTI, due to the fact that the output spectrum contains frequencies that do not exist at the input spectrum..., I've looked at the input-output spectrums, and came to the following conclusion :
Given the input:
$$ s(t) = \begin{cases} { 1 ~~~,~~~-0.5 < t < 0.5 \\ 0~~~,~~~\text{otherwise} }\end{cases} \tag{1}
$$
and the output:
$$ g(t) = \begin{cases} { 2t+1 ~~~,~~~ -0.5 < t < 0 \\ 1-2t ~~~,~~~ ~~~~~~~~0 < t < 0.5 }\end{cases} \tag{2}
$$
The corresponding spectrums (apart from linear scale factors !) are :
$$ S(\omega) = \frac{ \sin(\omega /2) }{ \omega} \tag{3}$$
and
$$ G(\omega) = \left( \frac{ \sin(\omega/4) }{ \omega} \right)^2 \tag{4}$$
Now, the input spectrum contains all frequencies except the zeros of the sine function which occur at $\omega = 2\pi k$, whereas the output spectrum contains all frequencies except those $\omega = 4\pi k$. Therefore, RBJ is correct in his claim that this particular input-output pair violates the LTI principle of output cannot contain frequencies that do not exist at the input.
However, in my opinion, redefining the output triangular-pulse duration as
$$ \hat{g}(t) = \begin{cases} { 2t+1 ~~~,~~~ -1 < t < 0 \\ 1-2t ~~~,~~~ ~~~~~0 < t < 1 }\end{cases} \tag{5}
$$
eliminates the associated problem, as the output spectrum becomes:
$$ \hat{G}(\omega) = \left( \frac{ \sin(\omega/2) }{ \omega} \right)^2 \tag{6}$$
which includes all frequencies except at those $\omega = 2\pi k$. Thus the missing frequencies at the input are also missing at the output.
Hence, I would like to assume that there was a typo at the duration of the triangular-pulse output $g(t)$, and that otherwise the analysis based on LTI remains valid... Thanks to RBJ for the reminder... EDIT-1-ENDS
EDIT-2 :
However, now happens another problem. If the system is to remain LTI (so that subsequent analysis remains valid), then $\hat{g}(t)$ (the redefined $g(t)$) should be used in the remaning progress. But then the observed waveforms in the figures will be different from what's being drawn. Most notably, the step-response obtained in part-b, becomes a flat output (a unit-step with finite length triangular tails), and finally the impulse response found as step-response's derivative will be a single rectangular-pulse at the origin (and one more at positive infinity) and not a periodic waveform.
So I think the problem should be completely redefined to make any sense...
EDIT-2-ENDS
Your system is stated to be causal and LTI. And an example response for a single rectangular-pulse input $s(t)$ is provided to be a single triangular-pulse output $g(t)$:
$$g(t) = \mathcal{T}\{s(t)\} = s(t) \star h(t) $$
where $h(t)$ is the impulse response of the LTI system, and $\star$ denotes convolution.
The followings are easily deduced from the properties of LTI systems:
a-) Observe that $s_1(t)$ is given by a sum of two shifted $s(t)$ pulses side by side as:
$$ s_1(t) = s(t-0.5) + s(t-1.5)$$
then from the LTI property, conclude the output as:
$$g_1(t) = \mathcal{T}\{s_1(t)\} = \mathcal{T}\{s(t-0.5) + s(t-1.5)\} = g(t-0.5) + g(t-1.5) $$
the output for $s_1(t)$ consists of two triangular pulses repeated (and the origin is shifted).
b-) Now, instead of two pulses, consider an infinity of such rectangular pulses shifted to locate side by side. Then the output will also be infinity of triangular pulses located side by side, as given by the figures in the case b:
$$
\epsilon(t) = \sum s(t -0.5-n) \implies h_{\epsilon}(t) = \mathcal{T}\{ \epsilon(t) \}$$
$$
\begin{align}
h_{\epsilon}(t) &= \mathcal{T}\{ \sum s(t -0.5-n) \} \\
&= \sum \mathcal{T}\{s(t -0.5-n) \} \\
&= \sum g(t-0.5-n)
\end{align}
$$
incidentally, the signal $\epsilon(t)$ is a unit-step function, and the corresponding response $h_{\epsilon}(t)$ is known as step-response of an LTI system.
c-) then the impulse response of the LTI system can be found based on the relation of step-reponse to impulse-response. Namely, if $h_{\epsilon}(t)$ is the step-response of an LTI system, then the impulse response $h(t)$ is :
$$ h(t) = \frac{ d h_{\epsilon}(t)} {dt} $$
you can see from the figure in c-) that the impulse response $h(t)$ is indeed the differentiation of the step-response; i.e, the triangular pulses are converted into bipolar rectangular pulses, and they are repeated indefinetely at the rate of triangular pulses in the step-response. Hence, the figure in case-c.
or it could as well be directly verified as the differentitation property of convolution. Actually this property is a consequence of the differentiation property of the convolution operator...


Time-frequency synchronization algorithm for OFDM continuous stream I am studying on the implementation of a continuous stream OFDM baseband modem. Googling and also searching here on stackexchange for synchronization techniques I found the descriptions of several existing algorithms, depending on the specific system features, but no detailed implementation of a time-frequency synchronization algorithm for continuous stream systems (like DVB-T/T2, etc).

Can anyone suggest any link/reference? Any other information will be welcome of course.

Also, many technical papers report synch algorithms for packet oriented systems like 802.11 or WiMax. A continuous stream system should need a less complex synch technique, with the possibility to average the received signal metrics over the time, is my assumption correct?

A good starting point for synchronization in generalized OFDM is: "On Synchronization in OFDM Systems using the Cyclic Prefix", Van De Beek et. al.

AWGN noise will corrupt your symbols, resulting in an error for your synchronization estimates.  You could leverage the fact that you are streaming data, and average the symbol timing offset (STO) and carrier frequency offset (CFO) over multiple symbols in the stream using the cyclic prefix.  Or instead of averaging, you could use a proportional-integral (PI) feedback controller to constantly update your estimates for each new symbol in the stream.  In steady-state, you should have a pretty good estimate of the actual STO and CFO.

EDIT:
Here is a great review paper on the subject.  The paper can be found on Google Scholar.
Furthermore, it describes synchronization using PI control on pages 1407-1409.

[1] Morelli, Michelle, C-CJ Kuo, and Man-On Pun. "Synchronization techniques for orthogonal frequency division multiple access (OFDMA): A tutorial review." Proceedings of the IEEE 95.7 (2007): 1394-1427.


Picking the correct filter for accelerometer data I am fairly new to DSP, and have done some research on possible filters for smoothing accelerometer data in python. An example of the type of data Ill be experiencing can be seen in the following image:

Essentially, I am looking for advice as to smooth this data to eventually convert it into velocity and displacement. I understand that accelerometers from mobile phones are extremely noisy.

I dont think I can use a Kalman filter at the moment because I cant get hold of the device to reference the noise produced by the data (I read that its essential to place the device flat and find the amount of noise from those readings?)

FFT has produced some interesting results. One of my attempts was to FFT the acceleration signal, then render low frequencies to have a absolute FFT value of 0. Then I used omega arithmetic and inverse FFT to gain a plot for velocity. The results were as follows:

Is this a good way to go about things? I am trying to remove the overall noisy nature of the signal but obvious peaks such as at around 80 seconds need to be identified.

I have also tired using a low pass filter on the original accelerometer data, which has done a great job of smoothing it, but I'm not really sure where to go from here.
Any guidance on where to go from here would be really helpful!

EDIT: A little bit of code:

for i in range(len(fz)):
    testing = (abs(Sz[i]))/Nz

    if fz[i] < 0.05:
        Sz[i]=0

Velfreq = []
Velfreqa = array(Velfreq)
Velfreqa = Sz/(2*pi*fz*1j)
Veltimed = ifft(Velfreqa)
real = Veltimed.real

So essentially, ive performed a FFT on my accelerometer data, giving Sz, filtered high frequencies out using a simple brick wall filter (I know its not ideal). Then ive use omega arithmetic on the FFT of the data.
Also thanks very much to datageist for adding my images into my post :)

As pointed out by @JohnRobertson in Bag of Tricks for Denoising Signals While Maintaining Sharp Transitions, Total Variaton (TV) denoising is another good alternative if your signal is piece-wise constant. This may be the case for the accelerometer data, if your signal keeps varying between different plateaux.

Below is a Matlab code that performs TV denoising in such a signal. The code is based on the paper An Augmented Lagrangian Method for Total Variation Video Restoration. The parameters $\mu$ and $\rho$ have to be adjusted according to the noise level and signal characteristics.

If $y$ is the noisy signal and $x$ is the signal to be estimated, the function to be minimized is $\mu\|{x-y}\|^2+\|{Dx}\|_1$, where $D$ is the finite differences operator.

function denoise()

f = [-1*ones(1000,1);3*ones(100,1);1*ones(500,1);-2*ones(800,1);0*ones(900,1)];
plot(f);
axis([1 length(f) -4 4]);
title('Original');
g = f + .25*randn(length(f),1);
figure;
plot(g,'r');
title('Noisy');
axis([1 length(f) -4 4]);
fc = denoisetv(g,.5);
figure;
plot(fc,'g');
title('De-noised');
axis([1 length(f) -4 4]);

function f = denoisetv(g,mu)
I = length(g);
u = zeros(I,1);
y = zeros(I,1);
rho = 10;

eigD = abs(fftn([-1;1],[I 1])).^2;
for k=1:100
    f = real(ifft(fft(mu*g+rho*Dt(u)-Dt(y))./(mu+rho*eigD)));
    v = D(f)+(1/rho)*y;
    u = max(abs(v)-1/rho,0).*sign(v);
    y = y - rho*(u-D(f));
end

function y = D(x)
y = [diff(x);x(1)-x(end)];

function y = Dt(x)
y = [x(end)-x(1);-diff(x)];

Results:


Determining channel frequency response from measured data via IFFT I've been working on a small system that transmits a linear chirp over the 35Hz-20kHz frequency range in a room and simultaneously records the transmitted signal with echoes (from the floor, walls, etc.), using two microcontrollers.
Using the received data $y[n]$ and reconstructing the transmitted chirp $x[n]$ in the receiver, I obtain the channel frequency response by $H[k]=Y[k]/X[k]$, $X$ and $Y$ being the FFTs of $x$ and $y$ respectively. I then attempt to compute the channel impulse response by computing $h[n]={\tt IFFT}(H[k])$. I then try to use time-gating to get rid of the echoes.

The fact that I am not using external memory limits my computational ability to using float arrays of length N=2048 at most, so the signals I'm dealing with are in the single milliseconds regions (this of course limits my resolution, but that's a problem I have to come to terms with).

The problem that's bothering me is, when I look at the calculated impulse response I get something like this:

It seems that the main response is seemingly moved cyclicly a bit to the left, so that its beginning is at the end of the array and there is a large "spike" there (either that or the response has a big delay, which doesn't fit the physical system as the speaker and the microphone are just a few centimeters apart for the tests I'm performing).
What might be the reason for this?
I've read some posts around here which looked similar and there was talk about pre-ringing; could this be the case here?

I want to compute the FFT of the time-gated (and ideally, "fixed") impulse response to get some data out of it, and a part of it is obtaining the delay of the main response- can this be done in this case?

Insight would be greatly appreciated!

Avner

Edit: Attempting to bump this along with a follow-up question I had to Peter's response.
If I should attempt to ask this in a separate question instead- please let me know.

"Thank you very much! So far I've been transmitting a chirp that lasts 25ms, and the reception goes on for 32ms (so there's an underlying assumption that the channel is 7ms long due to short distances between speaker and mic). I've padded the chirp samples with zeroes to match the reception length (in samples) and using this scheme gave me the graph attached in the question.

I've read that using a sine-sweep and performing a deconvolution results in separating nonlinear harmonics to "behind" the channel response in time, so that might be what's moved cyclicly and runs over the channel. Do you think there is a solution to this problem, given the limitations? I'm pretty much filling up the 32kB memory on the microcontroller I'm using with two 2048-float arrays that constitute a general array I'm working on and take up 16kB (doing things in-place) and two more int-arrays that take up 12kB more, so there's not a lot of memory left to use."

Without seeing your code directly, my guess is that you're doing something like:

N = length(x);
X = fft(x,N);
Y = fft(y,N);
H = Y./X;
h = ifft(H);

The problem with doing that is this operation performs circular convolution over a length of N.  The trouble is, the "true" length of the convolution of y with 1/X's impulse response will be of length 2N - 1 (assuming x and y are both of length N).

So you are getting time domain aliasing.

The solution is probably to set

N = length(x) + length(y) - 1

or greater.


Band-limited signal recovery with finite bit depth The sampling theorem tells us that a signal with no frequencies above $f$ can be completely described by sampling it a rate of $2f$. However, the theorem makes no reference to quantization, and so I would think that the result only holds for samples measured with infinite precision.

However, in DSP we never have infinite precision, we instead have a fixed (or at least bounded) bit-depth. So what is the relationship between the bit-depth of each sample and the precision with which we could reconstruct the original signal? Does grossly over-sampling allow us to heavily quantize each sample?

For the purposes of the question, we could even imagine taking countably infinite samples from a signal that is infinite in the time domain, but each sample is itself limited to a resolution of $n$ bits.

Quantization is a loss of information. So if you quantize a signal, you will never be able to reconstruct it perfectly.

That being said usually we can, in certain circumstances, assume that  quantization is a kind of noise. In these circumstances, assuming an N-bit ideal ADC with a full-scale sine wave as an input.

$$SNR (dB) = 1.76 dB + 6.02N$$

So an ideal 16-bit AC would give you a dynamic range of about 98 dB. The hard question, is whether or not that's enough for your application. Unfortunately we cannot answer that for you.

Now, as to oversampling. If the quantization noise is not correlated, then yes oversampling can help you decrease the quantization noise, and thus increase the dynamic range in a certain frequency band. The rule-of-thumb is that oversampling by a factor of 4, then filtering, then downsampling by 4 is equivalent to adding another ADC bit, i.e. adding 6.02 dB to your dynamic range.

However, sigma-delta ADCs also perform quantization noise-shaping which allows them to vastly increase their dynamic range in a certain frequency band. Combined with oversampling, this allows them to implement 24-bit and 32-bit ADCs.

For more information on sigma-delta ADCs https://www.analog.com/en/technical-articles/behind-the-sigma-delta-adc-topology.html


Impact of the WSOLA analysis increment on the time stretching factor WSOLA, basically, is the same as OLA but with a dynamic hop size. That part is clear to me. What is not so clear is the following: The stretching factor is determined by alpha = Hs / Ha and WSOLA uses a small increment, obtained through cross-correlation, to determine where the next input frame will be extracted. If that's the case, doesn't that small increment introduce a deviation between the requested alpha and the real one? After all, that increment determines a new analysis hop size, Ha' = Ha + incr, meaning that, if we maintain the same synthesis hop size, the alpha will change.

Yes, there can be jitter between the requested alpha and the one provided by WSOLA; but that jitter will average out over the long run, and the peak jitter is bounded by on the order of twice the lowest frequency (half the longest period) accepted by the autocorrelation algorithm.


Inverse complex FFTW transform I've been scanning all of the FFTW documentation, trying to figure out how to inverse FFT my FFT spectrum.
The documentation only mentions how to inverse FFT real-to-complex transformations, using the opposite c2r transform.
There are no FFTW_INVERSE flags anywhere but I suspect that the FFTW_FORWARD/BACKWARD transform might be what I need.
i.e. FORWARD for FFT and BACKWARD for IFFT. The documentation speaks about it flipping the sign; I remember something similar happening with the IFFT.
Is the backwards transform just an IFFT (for forward transform)? If not, how do you achieve an IFFT?

Is the backwards transform just an IFFT (for forward transform)?

Yes.
Compare the "what FFTW actually computes", here


how to find the "Natural and Total Response" of the system Hello there I am a Student.

I need to know how to find the Natural and Total Response of the given system:

where: $x(t)=e^{-t}u(t)$ and $y(0^+)=2, \dfrac{d}{dt}y(0^+)=3$. (If the image above is not clear)

This is for my exam preparation, please kindly help or at least specify where
I can find the solution for this question, and  best tutorial for signal and system.

I'll assume you do do what a Laplace transform really is!

All you have to do, is just compute the la place inverse for each part to get the result you want.
If you want the natural response only, just find the inverse Laplace transform of the noted part. If you want the total response just take the inverse Laplace transform for the whole factor.


Computer Vision - Does implementation of Hough Transform for lines require a loop for angle parameter? I am trying to understand how Hough Transform is implemented in real world applications. Consider a line with equation:

$y = mx + b$

In parametric form, this line can be represented using (d, $\theta$):

$d = xcos(\theta) + ysin(\theta)$

This is the Hough Transform algorithm in simple form, using the parametric form (ref: link) :

Initialize parameter space H[d, theta]=0

for each edge point I[x,y] in the image,

   for theta = [theta_min to theta_max ],

      set H[d, theta] += 1

Find the value(s) of (d, theta) where H[d, theta] is maximum

1) In real world implementations, is the line: for theta = [theta_min to theta_max ] required? Since we apply Canny edge detector to an image, prior to Hough Transform, we can obtain angle of gradient($\theta$) for a point (x, y) from Sobel operators. This can be noted down against (x, y) and later used in Hough Transform algorithm, eliminating one for loop. Isn't this the better, real world approach? Or do we actually implement the extra for loop?

In real world implementations, is the line: for theta = [theta_min to theta_max ] required?
. . .
Since we apply Canny edge detector to an image, prior to Hough Transform, we can obtain angle of gradient(θ) for a point (x, y) from Sobel operators. This can be noted down against (x, y) and later used in Hough Transform algorithm, eliminating one for loop. Isn't this the better, real world approach?

The standard Hough Transform algorithm does not assume that direction information is available for the image that is transformed. It basically accepts a set of point values (pixel location, pixel value) ($(x,y,v)$) at its input and converts lines to points. This is because, colinear points 9lines) contribute to large "accumulations" (counts) along the directions of the line they form.
So, a pixel of value $v$ at some location $x,y$ contributes $+v$ to all the Hough domain arcs that pass through it.
Now, if we consider the line formed by two neighbouring pixels $(x_0,y_0,v_0), (x_1, y_1,v_1)$ that both happen to have values $v_0=v_1=v$ (for simplicity), then this line segment's contribution along the directions that are parallel to the line is $2 \cdot v$ but along the vertical direction is still $v$ for two arcs.
As the length ($l$) of the line segment that is formed by any two $(x_m,y_m), (x_n,y_n), m \neq n$ increases, even if the pixels are not neighbours now, then they contribute $l \cdot v$ to that one direction of convergence but still $v$ to $l$ arcs in the vertical direction. But also, as that length increases, it starts contributing between $v .. l \cdot v$ to a wider range of $\theta$s around the direction that is parallel to the line.
Therefore, the answer to the first question is "yes". The answer to the second question is "under certain conditions".
This is because, you can use some of the information of the Canny edge detector to inform the accumulation process of the Hough Transform, but you have to do it carefully.
The Canny edge detector implies a scale parameter because it uses a Gaussian filter. This scale parameter, along with the gradient information can be used to further subdivide the Canny edge detector's output into a set of composing elementary line segments at direction $\theta$ and length $l$. Then, the Hough Transform would be approximated by adding together the individual Hough Transforms of those elementary lines. But of course, notice here that you still have to "scan" the elementary line segments at a range of $\theta$s to obtain their contribution during this pre-computation stage.
Hope this helps.


Has this Principal Component Analysis (PCA) been done correctly? I have a set of 3D data points, indicated by the blue color in the picture below.

I then project them onto the x-y plane, i.e. setting z values of all the points to 0, shown by the yellow color below.

Then, I perform PCA on the 2D data points (the yellow ones).

I end up with the 1st Principal Component being the green line, and the 2nd being the yellow line.

As I am a newbie in PCA, I am not sure whether I have done something wrong or not.

Because in my intuitive, the 1st Principal Component should be like in the picture below:

I mean unlike the above picture given in wiki, my 1st component is NOT obvious and seems to make no sense. Did I do something wrong?

Edit:

To avoid the confusion, the 2D view from the above is as follows:

If the PCA was correct, than green and yellow line would be perpendicular. So something is wrong with your result.


Compute real signal from a discrete analytic signal I have a 128MHz-wide down converted signal that gets processed in a FPGA via a polyphase filter bank to give 8x16 MHz baseband analytic signals. How do I convert this analytic signal to a real-valued signal? The Wikipedia article on analytic signal mentions that interpolation/up sampling will be required if the original signal was down converted. Here is how I plan to proceed:

Extract the real part of the analytic signal
interpolate (up-sample by a factor of 2, followed by a low pass filter)

Would the above give me a real-valued signal?

EDIT 1: I believe, I have to up-sample (i.e add zeros to every alternate point) and then use a low pass filter, e.g. a FIR to interpolate.

@user85858302: If designing a complex-valued bandpass filter is the solution to your problem, you can design such a filter without using the Hilbert transform. Start by using your favorite method for designing a real-valued N-tap FIR lowpass filter whose one-sided bandwidth is B Hz. That will give you N real-valued coefficients. Next, multiply those coefficients by a complex exponential whose frequency is Fc Hz. (Fc must be sufficiently less than half your data’s Fs sample rate.) That will produce N complex-valued coefficients that can be used in a tapped-delay line FIR filter block diagram. The complex-valued bandpass filter will have frequency response whose center frequency is Fc Hz and a bandwidth of 2B Hz.  Model this whole process in software to demonstrate what I'm sayin' here.

For those readers interested in both real-valued and complex-valued FIR filters, I suggest you have look at the following blog. You may well learn something new about FIR filters. http://www.dsprelated.com/showarticle/808.php


Hamming Window, python 2.7 Hi I have a FFT which is quite noisy. How to apply to my code Hamming window to make it less noisy. Look at my code:

# Loop for FFT data
for dataset in [fft1]:
    dataset = np.asarray(dataset)
    psd = np.abs(np.fft.fft(dataset))
    freq = np.fft.fftfreq(dataset.size, float(300)/dataset.size)
    plt.semilogy(freq[freq>0], psd[freq>0]/dataset.size**2, color='r')

for dataset2 in [fft2]:
    dataset2 = np.asarray(dataset2)
    psd2 = np.abs(np.fft.fft(dataset2))
    freq2 = np.fft.fftfreq(dataset2.size, float(300)/dataset2.size)
    plt.semilogy(freq2[freq2>0], psd2[freq2>0]/dataset2.size**2, color='b')

I put an answer in the questions.

In fact Welch was a good idea. End of post. Problem solved.

# Loop for FFT data
for dataset in [fft1]:
    dataset = np.asarray(dataset)
    freqs, psd = welch(dataset, fs=266336/300, window='hamming', nperseg=8192)
    plt.semilogy(freqs, psd/dataset.size**2, color='r')

for dataset2 in [fft2]:
    dataset2 = np.asarray(dataset2)
    freqs2, psd2 = welch(dataset2, fs=266336/300, window='hamming', nperseg=8192)
    plt.semilogy(freqs2, psd2/dataset2.size**2, color='b')


Recognizing bodily noises (ie, not spoken words) I was pointed here by a chain of stackoverflow posts. I'm a software developer looking to make a module that recognizes noises like snorts and coughs. Here was my original question on stackoverflow and here is where I found this SE board.

In short, I'm trying to use a Raspberry Pi and a microphone to pick up noises like snorts, coughs, sneezes, etc as a proof of concept to tell if someone is sick. I'm not quite sure where to start, as those posts I linked above have a lot of information. I've also come across this post, but I'm not sure if that helps.

I'm very new here and to this particular concept, so any and all help would be appreciated.

Most likely you are interested in a very simple approach that will run on your Pi. Some of methods I am going to mention are possible to run in real-time with not-badly written C code. Probably you will start with HMM's and MFCC's - it's not a bad idea for the beginning, but as you will see - people developed more advanced methods for such application. Choice of features is rather obvious, pre-processing and Machine Learning is up to you - this will have biggest impact onto your results.
Here is a list of some papers you might find interesting. All of them are targeting cough recognition problem and are applicable to other types of sounds.
HMM-like approach:

Takahashi S., et al. - Cough Detection in Spoken Dialogue System for Home Health Care
Shasha Le, Weiping Hu - Cough sound recognition based on Hilbert marginal spectrum
Matos S., et al. - Detection of cough signals in continuous audio recordings using HMM's
Chunmei Z., et al. - Recognition of cough using features improved by sub-band energy transformation
Hu W., et al. - Cough sound detection bases on EMD analysis and HMM recognition

Karhunen–Loève-like approach:

Liu Y., et al. - A cough sound recognition algorithm based on time-frequency energy distribution

Decision Trees-like approach:

Larson. E, et al. - Accurate and Privacy Preserving Cough Sensing using a Low-Cost Microphone
Martinek J., et al. - Distinction of cough from other sounds produced by daily activities in the upper airways

SVM-like approach:

Jia-Ming L., et al. - Cough signal recognition with Gammatone Cepstral Coefficients

And this is absolutely lovely ;)


Chop out frequencies outside human hearing range I have a bunch of audio files all sampled at 44100 Hz sample frequency. I am trying to remove all the frequencies which are outside the human hearing range (I use the following as reference: Frequency Range of Human Hearing, as well as robert.b's answer for determining the frequency value of a given bin: How to get Frequency from FFT result), and my basic approach is:

Perform an FFT on the files
Convert all frequency bin indices to actual frequencies
Remove those which are outside the range (i.e. below 20 and above 20,000 Hz)

Would that be a valid approach, and would it be possible to "put back" my signals to a state (using the inverse FFT) where I can play them and check if they still sound "normal"?

FFT -> zeroing coefficients -> IFFT is not a good approach for implementing filters. This will lead to a filter with poor performances and/or complicated implementation and/or inefficient implementation. This has been frequently discussed on this site, for example in this question.

Look into digital filters instead.


Symbol Synchronization for GFSK/FSK signals in GNU Radio I am trying to demodulate GFSK signal, and extract the packets using GNURadio.

TX information:

It is basically BLE transmission but customized. BW: 2MHz, data rate: 1MSymbols/sec, 1 bit per symbol. Each packet is 160 bits (which equals to 181us duration). And the packets are repeated every 100ms, hence approx. 10 packets per sec

RX information:

I am using a BladeRF with 4Msps and storing it to a .complex and read back in GNU Radio.
The samples are given through an LPF performing decimation by 2 which reduces the sample rate to 2Msps. Then I use a quadrature demodulation block following which I use a Symbol Sync block. After this I use a custom OOT module, to count the packets.

So, everything is perfect until the demodulation, after demodulation I can clearly see the data stream. Only thing is that we can't make a hard decision here  (we have 2 samples per symbols). In order to fix that we pass it through the Symbol Sync block. Then we do packet detection

Problem:

My main issue is not every packed is synced/recovered by this block. As most of the symbol timing recovery is a kind of estimation/control loop algorithm. Setting the correct parameters is very important. I played a lot with the parameters. And whatever I have set now works best for me.
Even in this best-case scenario, I am losing some packet. The capture file I have is for 32.03 seconds, which approx. contains 319 packets, but I can only recovery 311 packets.

I have been experimenting with this for nearly a month, with different packet rates, like 10, 25, 50, 75, 100 packets per second. All the test capture is done in a very low noise environment, and RX close to TX. So, any channel effects can be ignored. I even tried inserting a matched filter before demodulation, use a polyphase clock recovery, etc. but still, packets are lost.

And every time packet loss rate is always 5% to 7%

Questions:

1) why exactly this happens? is it because that the limit of symbol recovery can do or Am I doing something wrong?

2) is there any other alternative for this problem, to recovery 100% or at least with less than 4% loss? I have seen a technique which reverse engineer like M.ossmann's clock recovery use (https://github.com/mossmann/clock-recovery). But, I want to do it in real-time.

Any other suggestions to modify the algorithm or parameter or even approach are highly welcomed and appreciated....!!! I stuck with this problem for soo long and it's getting on my nerves.

Please, I am calling out experts in the field for help and looking for discussions.

Thanks in advance

P.S: I am attaching my screenshot of the GNU Radio flow graph. And also a plot showing how some packets are distorted after passing through the symbol sync block. The red is the bits of a packet after the symbol sync block and blue is before that. You can see that it's been lost many bits. Hence, the packet detection block is unable to detect it.

UPDATE:

Based on the suggestion from the comments, apparently the loop bandwidth and damping factor value are not very optimal. As my preamble is 8 symbols the loop should settle to linear operation within this time. This corresponds to the loop bandwidth of ~40KHz (normalized loop BW 20m). After these changes also I was still getting a 5% packet loss. In order to analyze a bit deeper, I created a bit more detailed graph with the debug signals from the Symbol Sync block. I have superimposed the recovered packet with error, instantaneous period, average period debug signals. I am attaching the graph for both cases, where the packet is recovered/detected good and the failure case. The first image corresponds to a good packet, the second image is for a bad packet.

The Symbol Synchronizer block is a PLL-based synchronizer that is trying to estimate the symbol clock period and symbol clock phase (aka timing offset) based on the samples coming in that represent the data symbols.

Being a PLL configured with static parameters, there is a fundamental trade off between acquisition speed and tracking stability of the symbol clock estimate.  The more reactive the loop configuration, the faster it will acquire, but the more likely it may loose track later and cause an error.  The more sluggish the loop configuration, the longer it will take to acquire, but it will track much better when locked.

You can speed acquisition by using a correlator out front (e.g. the corr_est block) that detects your preamble and marks the center of a symbol with a "time_est" tag.  The symbol sync block will notice the "time_est" tag (or a "clock_est" tag) and reset it's tracking estimates right there to rapidly acquire lock.

But that is not your immediate problem.

For the symbol sync block to work well, you must understand the requirements, assumptions, and behaviors of the Timing Error Detector you have selected.

The M&M TED, with the default "bpsk" decision constellation, is expecting:

symbols that are peaked in the center
symbols that ideally have peak(trough) levels of $\pm 1$
you to set an expected TED gain that comes from modeling and simulation of the TED in this particular situation (expected SNR, samples symbol in use, matched filter in use, etc.) in MatLab, Octave, Python or other M&S tool of your choice.

When you model and simulate the M&M TED, you may notice something peculiar about its behavior:

If you have a good estimate of the symbol clock period, then the M&M TED S-curve, and hence TED gain, collapses to $0$ on sequences of alternating symbols: $..., -1, 1, -1, 1, -1, 1, ...$

This is likely why you are failing to acquire lock at the beginning of some packets when starting with some random symbol clock phase (aka timing offset) estimate.

Either use a correlator that marks a preamble symbol to jump-start the timing synchronization, or use a different TED.  The M&M TED is old technology.

When using a different TED, make sure you understand its derivation and the assumptions that went in to its derivation.

Also design and modeling and simulation in MatLab, Octave, R, Python, or whatever, is absolutely essential if you want something other than suboptimal results.  GNURadio is a fine prototype tool, but is not a very good up front design tool.

Please read this brief, if you haven't already:

https://www.gnuradio.org/grcon/grcon17/presentations/symbol_clock_recovery_and_improved_symbol_synchronization_blocks/Andy-Walls-Samples-to-Digital-Symbols.pdf


Recovering signal for psk I am currently working on a product(cant disclose too much details of the actual project) in which some data is being sent via audible(but on the higher freq) sound that is encoded with dqpsk.

However, it gets very inaccurate at times since the receiver and the broadcaster does not have the same sampling windows. I am not trained academically for DSP and therefore having a lot of problem finding out what are the correct terms to search for.

Basically after a lot of debugging and simulation, i found that the most major issue is that the phase information gets more and more inaccurate as the output and the input gets more out of sync... Are there some techniques to find the beginning of the actual signal in the samples? The data isnt contained in just one sample window(they can span through quite a lot of windows)

Thanks! (I would appreciate any sort of help, even if it is just some term i should use to search on google... )

And on a seperate note, when the timing is relatively accurate, I am already able to get the data properly decoded, so I am not really having much problem with the psk side

To add to the excellent information given by Cassman in his response, here is a block diagram of a carrier recover loop for QPSK and QAM modems using a decision directed approach. I have detailed the decision directed phase detector in this post Phase synchronization in BPSK and this one How to correct the phase offset for QPSK I-Q data, while the block diagram below shows how it fits in a complete recovery loop (and demodulator).

The NCO puts out the complex frequency (Sine and Cosine for moving the carrier frequency offset in one direction) necessary to steer the input signal to 0 frequency (removes the carrier offset). I have more details on NCO implementations specifically here: Numerically Controlled Oscillator (NCO) for phasor implementation? and If the idea of complex frequency and "one direction" frequency shift is at all confusing, see this post: Frequency shifting of a quadrature mixed signal

The frequency setting is provided by the Loop Filter which accumulates the phase error up or down as needed to maintain the error at 0. As shown, this implementation is a Type 2 2nd Order Phase Lock Loop.

Where this block diagram below shows its use in a QAM demodulator. IF the Phase Detector just has one decision threshold instead of the 4 shown (for 16 QAM) then this will track and demodulate QPSK coherently.  As Dilip has said in the comment, you do not need coherent demodulation with DQPSK but this will work in either case.

Further this works really well with Gardner Loop for timing recovery (which requires 2 samples per symbol), as The Gardner Loop is not very sensitive to carrier offsets. (While the M&M synchronizer, another common timing recovery approach is quite sensitive).  The timing recovery would be ahead of the carrier recovery in the block diagram at 2 samples per symbol (or more), which is then downsampled to one sample per symbol at the sampling position in the center of the symbol for the carrier recovery implementation shown. Further details of the Gardner Loop are in this post: Gardner Timing Recovery for Repeated Sybmols. It's position in a receiver can be as done in the block diagram below, where depending on implementation it may be best to use the samples after a matched filter or before the matched filter.

This is a control loop and has a loop bandwidth that you set by the gain constants in the Loop Filter. For guidance on what Loop BW to use see this post: Loop bandwidth for symbol timing recovery. I have experimented with this specific to the Gardner and M&M recovery approaches and found that the Gardner performs best before the matched filter while the M&M performs best after, but that was specific to waveforms with raised cosine pulse shaping and different pulse shaping will impact timing recovery. Details on that comparison are given in this post: Location of Matched Filter


pitch detection: which is the most robust way to distinguish voiced from unvoiced audio I am trying to write a voice/whistle driven synth but I am struggling with pitch detection. As a premise I must say I am quite skilled and experienced at dsp coding plus I have pretty read almost everything I could find in literature about this very thorny topic and I understood that a single robust solution does not exist. I tried and invented every possible and fancy approach in frequency domain based on analysis of the harmonic series or HPS and variations thereof, but I am still plagued by octave errors,  plus the lack of a robust criterion to distinguish pitched from unpitched frames (i.e containing noise). At the end I decided to give up spectral approaches and use time domain autocorrelation (by FFT, squared magnitude plus IFFT) which turned out the most reliable approach of all even if still not 100% error prone. But again, I cannot find any reliable and robust rationale to decide when a frame is voiced/pitched and when it contains noise or transients and therefore the detected highest peak in autocorrelation has to be deemed invalid and discarded.

I found the sourcecode of a plugin called Autotalent which I judged interesting,  the author there relies on a "confidence" value obtained by multiplying thr height of the highest peak in autocorrelation by the value at the corresponding position of the autocorrelation of the window used to window the processed frames but I really failed to understand the logics behind that. I tried for curiosity to use the same trick but what I get is a value which seems completely unrelated to the fact a frame is pitched or unpitched.

I kindly ask where should I start looking for a solution which is not,  possibly,  too complicated (everything is meant for realtime operation). Any hints ?  Do you know any robust method to distinguish pitched and unpitched frames (either in time or frequency domain ) ? Thanks

SWIPE is a pitch estimator that worked well for me in the past. You can find dissertation & code in MATLAB here. It might give you some ideas.


Matlab: is a highpass filter what I want, to obtain my required spectrogram? How do I design/apply this? (pic included) The axes are wrong in the spectrogram. I have an audio wav file, and I know the sampling rate. this is read into the audioData variable. The audio is 1 channel. selected window length I believe (can't remember) is 512. I believe overlap is 50% ( the default?), 1024 fourier points. sampling rate is 250,000.

Using Matlab 2014a.

spectrogram made using this:
    s = spectrogram(audioData, 512, [], 1024, 250e3);

In the spectrogram, the bottom bit is near the 20kHz range and I want to get rid of this and concentrate on the thing at the centre (which is around 25khz). To my knowledge, I need a high pass filter to allow all of the frequencies above 20kHz.
It seems like there are many ways to design filters which involve a lot of parameters and I'm a bit overwhelmed with how to start.

I essentially just want to cut out the bottom bit, so that I can image process the upper bit (ideally kept as close to intact, if not exactly, as possible)

If you don't have any specification for the filter besides the cut-off frequency, I believe the simplest way is to use the fir1 function. Note that the frequencies are normalised to half the sampling rate. For example:

b = fir1(n,20e3/(250e3/2),'high'); %Wn = Cutoff_req/(Fs/2)
figure; freqz(b,1,[],250e3); %visualize the filter if you wish
FilteredData = filter(b,1,audioData);

If you wish to have more control over the specifications (attenuation, ripple, etc.), I suggest using fdatool.


Receiver function, frequency domain deconvolution not giving logic results I'm working on some code for receiver function method in seismology. For anyone one not into the topic, it's just a deconvolution of two time series (seismograms). This can be done in the time domain as well as in the frequency domain, however, the method I'm working on is doing it in the latter. The function I'm applying is this one, as described in Langston, 1979:

The phi term is what we call water level, and what it does is that it "fills" the low frequencies in the DFT so it can ignore the noise and smooth the result. The Gaussian filter is also intended for smoothing of the result.
Both R and V seismograms looks like this, an usual seismogram:

However, after applying the deconvolution I get something like this:

I'm very confused by those two high spikes at the beginning and the end of the result, it should be smooth and not spiky.
This is my code:

verfrec = sp.fft(stv.data, numpoints)[:numpoints//2]
radfrec = sp.fft(stradial.data, numpoints)[:numpoints//2]
frec = np.linspace(1, numpoints//2, numpoints//2)

gauss = np.exp((-frec**2)/((4*60)**2))

waterlevel = np.abs(verfrec * np.conjugate(verfrec))
waterlevel = np.maximum(waterlevel, 0.1*np.max(waterlevel))

erefrec = radfrec * gauss * np.conjugate(verfrec) / waterlevel

inverse = sp.ifft(erefrec)

fig = plt.figure(figsize=(12, 4))
ax = fig.add_subplot(1, 1, 1)
ax.plot(frec, inverse, 'sandybrown')
plt.title('Sismograma')
fig.show()

I'm splitting the transforms range in the middle since the real input will produce symmetric outputs. The 60 in the gauss definition is the gauss parameter and it should be between 0.1 and 1.5 or so, with those values I get wrong results so I'm guessing its a magnitude thing.
Any help is very appreciated.

The "spikes" that you see are coming from a "half-done" inverse Discrete Fourier Transform (DFT).
The main thing to note is that when you specify a filter in the frequency domain (as it is hapenning here), that filter has to "fit" the DFT spectrum ordering too.
You conclude (correctly) "I'm splitting the transforms range in the middle since the real input will produce symmetric outputs. " but that is not what is expressed in the code.
If you use numpy.fft.fft() then you have to specify a symmetric filter. Which is clearly not what gauss will be, because it is defined over freq whose linspace goes from 1 (why one?) to numpoints//2. If you plot this, you will see that it is only half the Gaussian and it has to be extended with its mirror image (since you are using numpy.fft.fft()) to include the negative frequencies part.
What you describe (that you only specify half the spectrum because the rest is symmetric) is what the pair numpy.fft.rfft(), numpy.fft.irfft() implement. But if you choose this way of doing it, you have to get rid of numpoints//2 because it is already taken care for you.
Having done this, we then turn to the "cut off" (or equivalent cut off) of the filter. The formulas in Langston79 are in angular frequency and the "standard deviation" (the width of the Gaussian filter) is specified at some $4 \alpha$, which is something like 6Hz (for an $a=1.5$). This has to be defined appropriately on the discrete frequency axis. If you have (for example) an $Fs=1000Hz$ and you are looking for 6 Hz in a 500 point signal, then after the DFT, each bin is $\frac{Fs}{N_{FFT}}$ Hz wide (or 2 Hz in this case). So your discrete $\alpha$ (let's call it $\alpha_d$) is now $3$. If your signal is very long (which I think it is, here, just looking at the duration of those signals), the $\alpha_d$ will work out much higher of course, because the resolution of the DFT is now higher (to still achieve the desired $4 \alpha$).
Hope this helps.


Cepstral Mean Normalization Can anyone please explain about Cepstral Mean Normalization, how the equivalence property of convolution affect this? Is it must to do CMN in MFCC Based Speaker Recognition? Why the property of convolution is the fundamental need for MFCC?

I am very new to this signal processing. Please help

Just to make things clear - this property is not fundamental but important. It is the fundamental difference when it comes to using DCT instead of DFT for spectrum calculation.

Why do we do Cepstral Mean Normalisation

In speaker recognition we want to remove any channel effects (impulse response of vocal tract, audio path, room, etc.). Providing that input signal is $x[n]$ and channel impulse response is given by $h[n]$, the recorded signal is linear convolution of both:

$$y[n] = x[n] \star h[n]$$

By taking the Fourier Transform we get:

$$Y[f] = X[f]\cdot H[f] $$

due to convolution-multiplication equivalence property of FT - that is why it's so important property of FFT at this step.

Next step in calculation of cepstrum is taking the logarithm of spectrum:

$$Y[q] = \log Y[f] = \log \left( X[f] \cdot H[f]\right) = X[q] + H[q]$$

because: $\log(ab) = \log a +\log b $. Obviously, $q$ is the quefrency. As one might notice, by taking the cepstrum of convolution in time domain we end up with the addition in cepstral (quefrency) domain.

What is the Cepstral Mean Normalisation?

Now we know that in cepstral domain any convolutional distortions are represented by addition. Let's assume that all of them are stationary (which is a strong assumption as a vocal tract and channel response are not changing) and the stationary part of speech is negligible. We can observe that for every i-th frame true is:

$$Y_i[q] = H[q] + X_i[q] $$

By taking the average over all frames we get

$$\dfrac{1}{N}\sum_{i} Y_i[q] = H[q] + \dfrac{1}{N}\sum_{i} X_i[q]$$

Defining the difference:

$$\begin{array}
&R_i[q] &= Y_i[q] - \dfrac{1}{N}\sum_{j} Y_j[q]\\
& = H[q] + X_i[q] - \left(H[q] +  \dfrac{1}{N}\sum_{j} X_j[q]\right) \\
& = X_i[q] - \dfrac{1}{N}\sum_{j} X_j[q]\\
\end{array}$$

We ending up with our signal with channel distortions removed. Putting all above equations into simple English:

Calculate cepstrum
Subtract the average from each coefficient
Optionally divide by variance to perform Cepstral Mean Normalisation as opposed to Subtraction.

Is Cepstral Mean Normalisation necessary?

It's not mandatory, especially when you are trying to recognise one speaker in a single environment. In fact, it can even deteriorate your results, as it's prone to errors due to additive noise:

$$y[n] = x[n] \star h[n] + w[n] $$

$$Y[f] = X[f]\cdot H[f] + W[f] $$

$$\log Y[f] = \log \left[X[f]\left(H[f]+\dfrac{W[f]}{X[f]} \right) \right] = \log X[f] +\log \left(H[f]+\color{red}{\dfrac{W[f]}{X[f]}} \right)$$

In poor SNR conditions marked term can overtake the estimation.

Although when CMS is performed, you can usually gain few extra percent. If you add to that performance gain from derivatives of coefficients then you get a real boost of your recognition rate. The final decision is up to you, especially that there are plenty of other methods used for the improvement of speech recognition systems.


How could one can say that frequency spectrum having Nyquist frequency 'x' is the most appropriate representation of a given signal? I have a sinusoidal signal which has  a period of 2 and I want to find the frequencies that are present.I assumed sampling interval as 0.01 and based on this I sampled the signal.The Nyquist frequency i got from this is suppose 'x'.Now if I change my sampling interval to 0.001, Nyquist frequency in this case is '10*x'
So my question is,which frequency spectrum will represent the signal correctly, the one with Nyquist frequency 'x' or the one one with Nyquist frequency '10*x'

The most important thing is to sample at more than twice your signal's bandwidth; the period is largely irrelevant. Say your signal with period 2 is a square wave -- its bandwidth is infinite and you'll never be able to sample it correctly (at least in theory). If your signal is a sine wave, then a period of 2 means its bandwidth is 0.5 and you'll easily sample it adequately.

Personally, whenever possible, I'd sample:

a bit above Nyquist (say, at $2.2B$, where $B$ is the signal bandwidth), to maximize the number of useful frequency bins, and
for as long as possible, to increase the spectral resolution.


Real-time resampling pitch-shifting (without corrective time-stretching)? Is real-time resampling pitch-shifting impossible (without corrective time-stretching, such as WSOLA)? I've been told that it is.

Real-time resampling pitch-shifting would mean applying resampling in real-time, but it's claimed that it's impossible to overlap differently resampled segments in real-time (since the total length of the "window" varies in every window and in connection to every frequency of resampling).

okay, the answer to the question is "no", from my understanding of what "real-time" in DSP means.  a real-time DSP process is one that can run forever without either starvation of input samples or overflow of input samples.  for every second of sound input, there is a second of sound output.  just speeding up the sound or slowing it down, to change the pitch, is an operation you can do on a segment of sound or a sound file.  but it is not a continuous operation that runs real time.  this is not because the CPU or DSP cannot handle the number of operations, it's because, either the input or the output will have no time when the other still has time remaining.

now, there are DAW functions that play back sound files faster or slower that recorded without changing pitch.  what they are really doing is playing back the sound file faster or slower (with an increased or decreased number of input samples) and the DSP algorithm is doing real-time and complementary resampling and time-scaling to get to the output.


How good or bad my downsampled data is? Suppose I have time series data at a one-minute resolution. Now I downsample data by taking mean of every 10-minute window, i.e., after downsampling, 60 readings will reduce to 6 readings. How should I show how good or bad my downsampling technique is? In a good downsampling technique, the downsampled sample should closely represent the original data.

First look at the spectrum of the original data $x[n]$. Downsampling a discrete-time signal $x[n]$ by $M=10$ requires that (to prevent aliasing) the signal should be bandlimited to $|\omega|<\pi/M$.

If your signal is bandlimited enough, then averaging 10 samples (as a crude lowpass filtering) will suffice. Otherwise, a better filter can be designed by different filter design techniques such as a simple windowed linear phase lowpass filter.

the following matlab excerpt gives you a better result.

K = 15;
h = fir1(2*K, 1/10); % impulse response of a 20th order FIR LPF.
yc = conv(x,h);      % convolve h[n] with x[n] of length L
y = yc(K+1:K+L);     % the filtered signal


Determine whether the system is linear? $T(x[n]) = ax[n] + b$

$ T(\alpha_{1}x_{1}[n] + \alpha_{2}x_{2}[n]) = \alpha_{1}ax_{1}[n]+ \alpha_{2}ax_{2} [n] +b $

$ \alpha_{1}T(x_{1}[n] ) +  \alpha_{2}T(x_{2}[n]) = \alpha_{1}(ax[n]+b) + \alpha_{2}(ax_{2} [n] +b) $

therefore they are not equal so the system is non-linear?

or would it be something like this:

$ T(\alpha_{1}x_{1}[n] + \alpha_{2}x_{2}[n]) = \alpha_{1}ax_{1}[n]+\alpha_{1}b + \alpha_{2}ax_{2} [n] +  \alpha_{2}b $

$ \alpha_{1}T(x_{1}[n] ) +  \alpha_{2}T(x_{2}[n]) = \alpha_{1}(ax[n]+b) + \alpha_{2}(ax_{2} [n] +b) $

and therefore they are equal so the system is linear?

In the following, I suggest that, before using the generic $T(\alpha_1 x_1+\alpha_2 x_2)$ versus  $\alpha_1 T( x_1)+\alpha_2T( x_2)$, it can be more informative to try with simpler partial tests, or try counterexamples, based on your intuition.

I don't really understand the motivation behind the second group of equations, or why $b$ gets multiplied by $\alpha_i$.

Here, you are using the generic version of the linearity test, which is good, but can raise some doubts, as apparent from your question. You can, in complement, use other simpler tests, that can show that the system is not linear, as a double check. Those are counter-examples. For instance:

is the output of the zero signal zero? For a linear system ,  $T(\vec{0})=0$.
is a single  output  with a zero coefficient  zero? This tests if $T(0.\vec{x})=0.T(\vec{x}) = 0$
is a single  output linear? This tests if $T(\alpha.\vec{x})=\alpha.T(\vec{x})$
is a simple addition   linear? This tests if $T(\vec{x}+\vec{y})=T(\vec{x})+T(\vec{y})$

Those, if not passed,  prove that the system is non-linear. And instead of using the generic version, they can show to your (clever) professor that you have some intuition about what is going on, and the risk of errors is reduced. I personally  appreciate a lot when students use minimal arguments: they go straight to the point, spend less time on such questions, to focus on more involved ones.

Of course, if the system is linear, more is required.

Here, your system is non-linear... unless $b=0$.


What are bits in the context of channel capacity units? 'bits' are used with 2 different meanings.  One can use 'bits' to mean the binary digits i.e. 1's and 0's.  Bits are also the units of information in a event of discrete source.  I think when channel capacity is measured in bits/sec, we use the second meaning of bits.  Not how many 1's and 0's we can send down the channel but how much information we can send per second through that channel.  Is my interpretation right?

in Information Theory, the amount of information in a message is proportional to the logarithm of the probability of occurrence of the message:

$$ I(m) = -C \log\left( P(m) \right) $$

the only difference between different bases of logarithms is a scaling constant.  if $ C = \frac{1}{\log(2)} $, that is

$$ \begin{align}
I(m) & = -\frac{1}{\log(2)} \log\left( P(m) \right) \\
     &= - \log_2\left( P(m) \right) \\
 \end{align} $$

then we say that the $I(m)$ is in units of bits.  the reason why is if the message $m$ had only two possibilities and both had equal likelihood, then $P(m) = \frac12 $ and $I(m) = 1$.  there are two messages, two possibilities of equal likelihood, like heads and tails of a coin.  and you need only one bit, 0 or 1, to represent the two possibilities.


Optimal DC-free line encoding I have a block of n bits that I'd like to encode into m bits that have an equal number of zeros and ones.

I looked at protocols like 8b/10b but it wastes 25% of the data to achieve balance.

As an example, it turns out that 128 bits have ${128 \choose 64}$ combinations where there is an equal number of zeros and ones.  This number is greater than $2^{124}$.  This means that I should be able to encode 124 bits of data into 128 bits that have an equal number of zeros and ones.

I have no idea how to do this mapping though.  Obviously I can't make a look-up table of size $2^{124}$.  Anyone have any ideas?

Your problem can be restated as the following:

"Find an easily computable bijection from the integers $0 .. {n \choose n / 2} - 1$ to the set of binary strings of length $n$ with exactly $n / 2$ zeros and $n / 2$ ones."

Formulated this way, this looks like a combinatorial numbering problem.

Decoding is quite simple and involves accumulating binomial coefficients for each bit set in the codeword. The encoding looks more complicated. Having enough ROM/RAM to store a big binomial coefficients table helps a lot in any case. You might want to repost on cs / math stackexchange, mentioning only the combinatorics aspect (not the digital communications application) if you are in search of more memory-efficient algorithms.

See python code for encoding and decoding here.

High-level explanation.

Let us consider the case $n = 6$ (3 ones and 3 zeros). The goal is to find a mapping between the strings 000111, 001011, 001101 up to 111000... and the numbers 0, 1, 2, 3 up to 19. This is equivalent to coding $\log_2 20 = 4.32$ bits into a 6-bit DC-free code.

There are ${6 \choose 3} = 20$ strings made of 3 zeros and 3 ones. ${5 \choose 3} = 10$ of them start with a 0 and ${5 \choose 2} = 10$ of them start with 1. So we can assign the numbers 0 to 9 to the first group, and 10 to 19 to the second group.

Within the first group, there are ${4 \choose 3} = 4$ which start with 00 and ${4 \choose 2} = 6$ which start with 01. So we can assign the numbers 0 to 3 to the 00 group and 4 to 9 to the second group.

Within the 00* group, there are ${3 \choose 3} = 1$ which start with 000 and ${3 \choose 2} = 3$ which start with 001. So we can assign the number 0 to the first one (000111), and the numbers 1 to 3 to the second group (001*). And so on...

Observe that in all the ${n \choose k}$ we have written, $n$ is always the number of remaining positions, and $k$ the number of unallocated ones. So the partitions of the $0 .. 19$ set we have built are created by adding ${\text{position at which we allocate a one} \choose \text{number of remaining ones}}$ terms.

Let us formalize this procedure. If the bitstring to decode is $b_5b_4b_3b_2b_1b_0$, the number we have associated to it is:

$b_5 {5 \choose b_4 + b_3 + b_2 + b_1 + b_0} + b_4 {4 \choose b_3 + b_2 + b_1 + b_0} + b_3 {3 \choose b_2 + b_1 + b_0} + b_2 {2 \choose b_1 + b_0} + b_1 {1 \choose b_0} + b_0$.

That's all for the decoding (DC-free bitstring to number)!

The encoding (number to DC-free bitstring) works by trying to "fit" the terms of this sum to the integer to encode. For example, given the number 15 to encode:

This number is above ${5 \choose 2}$ so the first digit is 1, and the first term of the sum is $1 {5 \choose 2} = 10$. 2 ones remains to be allocated. 5 remains in the sum.
The remainder is less than ${4 \choose 2} = 6$ so the next digit can't be a 1, so the second term of the sum is $0 {4 \choose 2}$. 2 ones remain to be allocated. 5 remains in the sum.
The remainder is above ${3 \choose 2} = 3$, so the next digit is 1, and the third term of the sum is $1 {3 \choose 2}$. 1 one remains to be allocated. 2 remains in the sum.
The remainder is equal to ${2 \choose 1} = 2$, so the next digit is 1, and the fourth term of the sum is $1 {2 \choose 1}$. All ones are allocated. All terms of the sum are accounted for.
The remaining terms are $0 {1 \choose 0} + 0$.

Thus, the bitstring associated to 15 is 101100.


Limiter, compressor, AGC, dynamics processor? What are the primary differences between the usage or definitions of the terms “AGC”, “compressor”, “limiter”, and dynamics processor? They all seem to be forms of playing with gain knobs.

one category has a protection property and has a very fast response to large amplitude transients.

another category attempts to maximize sensitivity.

another group tries to put the receiver in a predictable “sweet spot” or operating point.

some can be “slow” and others are “fast”.

none of these are mutually exclusive.

advances in converter bit resolution (dynamic range) moves these techniques along.

One might consider schemes that adapt to be related.


Playing sound in left and right speaker alternatively and storing output in array in matlab I have a play function that is called twice, each plays a channel (a column of the two dimensional array) in both speakers. What I am trying to do  is to play both channel at the same time. How do I silence a channel and set the  the elements in the corresponding column to 0. I also need an output array that im trying to make look like this when the data is output into it. The code is in the link below.
1     0
2     0
1     0
0     4
0     1
0     1

Fs = 44100;
N = 400
t = (1:N)*(1/Fs);
freq = 600;

sound_array = zeros(N,2);
%starts N x 2 matrix
sound_array(:,1) = sin(2*pi*freq*t);
sound_array(:,2) = sin(2*pi*freq*t);

%Sound 1
sound(sound_array(:,1), Fs);
%Sound 2
sound(sound_array(:,2), Fs);
%Play Consecutive
sound(sound_array, Fs);

If you want to play only one channel then you must set zeros at one of your matrix columns, i.e.:

sound_L = sound_array;
sound_L(:,2) = 0;

sound_R = sound_array;
sound_R(:,1) = 0;

So your code becomes:

Fs = 44100;
N = 400
t = (1:N)*(1/Fs);
freq = 600;

sound_array = zeros(N,2);
%starts N x 2 matrix
sound_array(:,1) = sin(2*pi*freq*t);
sound_array(:,2) = sin(2*pi*freq*t);

%Left and Right channels

sound_L = sound_array;
sound_L(:,2) = 0;
sound_R = sound_array;
sound_R(:,1) = 0;

%Sound 1
sound(sound_L, Fs);
%Sound 2
sound(sound_R, Fs);
%Play Consecutive
sound(sound_array, Fs);


How to get started analyzing audio signals for detecting cracks in a material? I have two processes, one is normal and the other with some cracks of material. And I have the acoustic signals got from the materials by a sensor. I want to compare these signals to see if this signal can be used to detect cracks. I think two different processes (with and without cracks) will create two different signals. Based on that, I want to compare these signals, but I'm new to DSP.

Could you please give me some suggestions what I should do, step by step to process these signals?

(I will of course read books to get knowledge about DSP, but it's really huge. I therefore hope that someone can tell me what I should do, step by step to process those signals. By doing so, I can save some time)

Not really a programming question!

In any case if you need to read books to understand DSP and are trying to compare two non stationary processes I suggest you start reading right away because most steps I tell you won't make sense and you need to answer a host of questions even before anyone can answer this question.

Signal processing deals with lot of statistical phenomena and you have to work with it emperically many a times. Not to say there are no fixed rules but the step-by-step to process this signals is hard enough when data is stationary. Given that you have non stationary data, get reading and once you have a hang of it ask some pointed questions on dsp.stackexchange.com


Maximum cross-correlation coefficient value for time delay estimation I am using cross-correlation for time delay estimation of two synchronized recordings ($x_1$ and $x_2$) of a fixed sound source from two different locations.

I understand that the delay is associated with the maximum cross-correlation coefficient.

My question: should I use the maximum of the cross-correlation coefficient vector or absolute cross-correlation coefficient vector?

This

$\Delta t \approx \tau_{peak} =  \underset{\tau}{\arg\max} (\rho_{x_1 x_2}(\tau))$

or this

$\Delta t \approx \tau_{peak} =  \underset{\tau}{\arg\max} (|\rho_{x_1 x_2}(\tau)|)$

As your plot shows, the second form allows for the correlation peak to be negative. Now, what does a strong negative cross correlation mean? It means the signals are very similar, except one has a negative sign in front of it, i.e., $x_1 \approx -x_2$. Whether or not this makes sense depends a lot on the actual application.

In the application you describe, it would totally make sense to me. See, signals that are transmitted from a source and travel to a destination through some medium are subject to changes in their phase, either due to reflections or the transmission itself. A minus sign corresponds to the opposite phase, which can very likely represent the signal you are looking for, only having gone through an appropriate phase change.

So if I were you, I would allow for the correlation to be negative and look for peaks in the magnitude of the correlation function. Eventually, it still depends on the exact requirements of your application.

If you have the chance, I'd suggest to try some calibration measurements with known source and receiver positions to see what's going on. Then you could even compare both approaches.


If a time-series has odd number of samples does it have no energy at Nyquist frequency? Suppose I have real time series A with n samples and time-spacing dt and I want to analyze its frequency content.
Af = fft(A)

If dt=1 Nyquist frequency is 0.5.
According to both Numpy and Matlab manuals, the frequency vector is defined as
f = [0, 1, ...,   n/2-1,     -n/2, ..., -1] / (dt*n)   if n is even

f = [0, 1, ..., (n-1)/2, -(n-1)/2, ..., -1] / (dt*n)   if n is odd

If n=9 (odd), f= [-0.44444 -0.33333 -0.22222 -0.11111  0.  0.11111 0.22222  0.33333  0.44444]
If n=8 (even), f= [-0.5   -0.375 -0.25  -0.125  0.     0.125  0.25   0.375]

Note the maximum frequency for odd number of samples is 0.44444 Hz not 0.5 Hz. Does this mean a time series with odd number of samples has no energy at Nyquist frequency ?

No, that's not what it means. You just don't get a frequency bin at Nyquist if the DFT length is odd. If you want you can always append one or more zeros to the time domain sequence such that the total length becomes even.
If you're just interested in the energy at Nyquist, you could also simply compute
$$E_{Nyq}=\left|\sum_{n=0}^{N-1}x[n](-1)^n\right|^2\tag{1}$$
where $N$ is the length of the sequence, no matter if it's even or odd.


Oversample or undersample correlated noise? I wish to measure nontrivial (non-Gaussian) broadband current noise for its distribution and power spectral density.
The noise is amplified (with transimpedance amplifier) with some characteristic risetime (tau ~30 usec), that I assume makes the signal correlated with the same timescale. The noise itself may be correlated with different timescale (~1 usec to 1 msec).

Is it best to oversample (sample frequency > 1/tau) or undersample (< 1/tau) the signal?
Does it really matter and what is the difference/effect?

I thought that slight undersampling draws directly from the distribution and the distribution is not biased to "double counting" (measure distortion) of the signal, however a colleague advised to overasample, thus no information gets lost.

To make it clear, the signal is random noise, and the purpose is to properly characterize it (so I'm not sure how Nyquist theorem fits in or relevant to it. Maybe other sampling theorems may be useful).

I would strongly suggest oversampling at a frequency much higher than the bandwidth of concern, and specifically setting up the test with proper filtering to ensure that there is insignificant energy beyond half the sampling rate or $f_s/2$ (assuming real sampling). This would imply using a low pass filter prior to the sampler and also reviewing the resulting power spectral density that the energy is sufficiently lower as you approach the Nyquist frequency of $f_s/2$.  The analog input frequency that extends to $f_s/2$ is referred to as the first Nyquist band.

You may be "bandpass" sampling in which case "undersampling" would apply with the same guidance above to "oversample" at more than 2x the bandwidth of interest with sufficient margin for filting. In this case the filter is a bandpass filter centered on the higher Nyquist zone (the first Nyquist zone for real signals is $f=0$ to $f_s/2$ as noted above, the second Nyquist zone is $f_s/2$ to $f_s$, the third Nyquist zone is $f_s$ to $3f_s/2$ etc all of which may be undersampled.

If this is still confusing, please see my post here with further details:

How do you simultaneously undersample and oversample?

However in all cases it is imperative to be filtering the signal of interest and to oversample it with a sampling rate > 2x the resulting bandwidth of the filter.
The reason is when you undersample below that bandwidth, while you will still be able to measure the overall total integrated power density, you will have lost all information as to what frequency your noise issues are at which would typically be of great interest in dealing with the noise. When sampling, everything that is in a band over a multiple of the sampling rate will fold into the same frequency band in your digital spectrum that extends from $-f_s/2$ to $+f_s/2$ (which is equivalently $f=0$ to $f_s/2$ for real signals).  For example, every frequency $Nf_s + \Delta$ will alias to $\Delta$ in your digital spectrum.

The easy way to see that the total integrated power density is conserved is to consider sampling white Gaussian noise and take the histogram of the results to visually see the bell curve and related standard deviation, the square of which is the variance proportional to the total power in the signal. From this visual it is easy to convince yourself that you can sample at ANY rate and with enough samples the histogram would have the same shape, standard deviation and variance. (This is true as long as the sampling clock is uncorrelated to the signal which would be the case for white noise.) But the content of the signal itself is all the energy at the input to the sampler at all frequencies where energy exists, so the actual analog frequency is not recoverable if not all isolated to one Nyquist band.


How to use 'randi' in MATLAB to generate equally probable samples which are either 0 or 1? How can I use ‘randi’ function in MATLAB to generate a signal sequence s(t) which contains 1e6 samples which are either 0 or 1 and are equally probable. Any suggestions? Thanks!

Generate random integer uniformly distributed between 1 or 2. Then you can shift it to your liking.
X=randi(2,size,1)-1


Finding the fundamental frequency from Autocorrelated data I'm writing an app in which I need to find the fundamental frequency of a note produced by a trombone. To do this I'm taking the FFT of audio data from a microphone and then using autocorrelation code taken from this question.

Once I've done that I'm then attempting to find the fundamental frequency by finding the highest and second highest peak locations and how many samples they are apart. I then find the amount of time one sample lasts with a sampling rate of 44100 and multiply that by the number of samples between the two peaks to give me the period of the fundamental frequency. Finally I use 1/period to find the fundamental frequency itself. Here's the code I've written for that (in Java):

public double findFundamentalFrequency(double[] autoCorrelatedData){

    double max = Integer.MIN_VALUE;
    double secondMax = Integer.MIN_VALUE;
    int maxLoc = Integer.MIN_VALUE;
    int secondMaxLoc = Integer.MIN_VALUE;

    for(int i = 0; i < autoCorrelatedData.length; i++){
        if(autoCorrelatedData[i] > max){
            secondMaxLoc = maxLoc;
            secondMax = max;

            max = autoCorrelatedData[i];
            maxLoc = i;
        }
        else if(autoCorrelatedData[i] > secondMax){
            secondMax = autoCorrelatedData[i];
            secondMaxLoc = i;
        }
    }

    double samplingPeriod = 1/44100.0;
    //Log.i("a", String.valueOf(maxLoc));
    //Log.i("b", String.valueOf(secondMaxLoc));
    double period = samplingPeriod * Math.abs(maxLoc - secondMaxLoc);
    double fundamentalFreq = 1.0/period;
    return fundamentalFreq;
}

Is this the correct way to find the fundamental frequency using autocorrelation? The answers this gives seem to be incorrect so I'm unsure if I've made a mistake somewhere or if I'm going about this the wrong way.

Three suggestions:

Try searching only over your expected pitch range: from a lag or period a bit shorter than that of your highest expected pitch for the instrument in question, to a lag a bit longer than twice the period of the lowest pitch expected.
Make sure there is a sufficient dip between your highest and your 2nd highest peak, or else you might find the shoulder of the 1st peak, and thus completely miss finding any 2nd peak.
You might also want to try using parabolic interpolation to better estimate each pitch period if any are between autocorrelation lags.


Frequency extracted from MATLAB's $\tt fft$ is not very accurate As a test I made up a sine wave in MATLAB of this form

y = 5*sin((2 * pi * freq).*x + 1.4) - 6;

where freq is 10 and x varies from $0$ to $1.5$ with a resolution of 1/1000 as shown below

fs = 1000;
x = 0:1/fs: 1.5 - (1/fs);

So I already know the frequency to be able to verify it with fft. After computing the amplitude FFT abs(fft(yy)), I find that the frequency bin with the highest magnitude is $16$. Since I have $1500$ samples which correspond to a sampling frequency of $1000$ then the 16$^\rm{th}$ bin corresponds to

$$\mathrm{\frac{Frequency \ Bin \times Sampling \ Frequency}{Number \ of\ Samples} = \frac{16 \times 1000}{1500} = 10.6667\ Hz}$$

However I know that my frequency I hardcoded is actually $10\ \rm Hz$. This can be repeated with different values and the same inaccurate result keeps occurring. And the smaller the hardcoded frequency the  larger the error in the result. Why is this happening?

MATLAB uses one based indexing, DFTs are zero based, so you have to shift your index value by one.

The bin index corresponds to the the frequency in units of cycles per frame.  Your freq is in units of cycles per second.  The sample count is in units of samples per frame and the sampling rate is units of samples per second. The bin index calculation should be:

\begin{align}
\mathrm{bin} &= \mathrm{freq} \cdot \frac{N}{f_s}\\
15 &= 10 \cdot \frac{1500}{1000}
\end{align}

Where the units are:

$$\rm \frac{cycles}{frame} = \frac{cycles}{second} \cdot \frac{\displaystyle\frac{samples}{frame}}{\displaystyle\frac{samples}{second}} $$

Bin 15 is in index 16 of MATLAB.  Since this is a whole number of cycles per frame, all your other bin values should be zero.

Hope this helps.

Ced

Follow up

Your calculation should have been:

$$\rm \frac{Frequency\ Bin \times Sampling\ Frequency}{Number\ of\ Samples} = \frac{15 \times 1000}{1500} = 10$$

Bin $-15$, aka bin $1485$, will also be non-zero and the complex conjugate of bin $15$.  Each should have a magnitude of $2.5$.  Your DC bin, aka bin $0$, should have a value of $-6$.  The corresponding MATLAB indexes will be $16$, $1486$, and $1$.


Normalized frequency in matlab I have a progam in matlab that generates a linear and quadratic chirp, and now I want to study the frequency components in wvtool. However, I am not able to decode the frquency domain scale "the normalized frequency". The usual process of deviding the normalized frequency with the sampling frequency doesnot provide the desired value. Please refer to the program and image for details.

Fs = 500;
tf = 1;
t = 0:1/Fs:1;
f1 = 10;
f2 = 90;
slope = (f2-f1)./t(end);

F = f1 + slope* t;
F1 = f1 + slope* t.* t;
% quadratic chirp
y = 1.33*cos(2*pi*F1.*t);
% linear chirp
y1 = 1.33*cos(2*pi*F.*t);

In order to get the actual frequency of each bin, not normalized frequency, you must compute the following:

fres = fs/fft_size; % Bin spacing
f = fres * (0:fft_size/2); % Bin centers and frequency axis


Are there other filter models besides IIR and FIR for processing digital signals? I am very familiar with IIR and FIR filters for processing digital signals (audio is my main interest). I am keen to know if there are alternative filters based on different mathematical models to IIR and FIR?

I have a degree in maths but learning about DSP is just a hobby (I didn't study it at university) so I haven't a clue whether this is a stupid question or not :-)

FIR and IIR filters are discrete-time linear time-invariant (LTI) filters. There are also (linear) time-varying systems, which could be implemented just like FIR or IIR filters, but with time-varying coefficients. Audio effects like chorus, flanger, or phaser are time-varying filters.

Of course there are also non-linear filters, the most well-known of which is probably the median filter, which is very useful for removing impulsive noise. Many different types of non-linear filters are used in image processing.

Note that none of the above mentioned filters can be characterized by their impulse response or, equivalently, by their frequency response. This is only possible for LTI systems.


how to detect signal noise I have data from a power transformer that has a frequency of 60 hz.  Sometimes there are faults in the system and causes a large spike in the reading of the transformer.  What I want to do is run the data into a signal processor to detect when these large spikes occur and either remove the spike value or present the time of occurrence.  how would I be able to detect this spike in Octave?

I'll assume that your data is in the form of discrete samples equispaced in time (as given by a standard analog-to-digital converter).

If the nominal signal is a clean 60Hz sine wave, I would suggest filtering out the 60Hz with a digital notch filter. This will essentially subtract the correct signal, leaving you with the noise. Then, compare the magnitude of the filter output to a threshold (which you will choose empirically) in order to detect the spikes. Note that the filter may or may not introduce a delay between the original signal and the filter output (esp, if you are doing this in realtime), which you will need to compensate in your time-of-occurrence estimate. This time compensation may also apply to the methods below.

If the nominal signal is not a clean sine wave, you can probably track the nominal frequency characteristics of the signal, e.g. using short time fourier transform, and watch for undesirable deviations from this.

You may also be able to perform a linear prediction on the signal, and detect spike noise events based on the magnitude of the prediction error signal.

These are all just suggestions, many other approaches likely exist.


Digital Lowpass butterworth filter with cut off 500Hz and sampling rate 1.25MSPS I am trying to simulate a distributed sensing system and I need to filter only frequencies lesser than 500 Hz(low pass filter) from acquired signal with sample rate 1250000 Samples/sec using below mentioned program:

 %Time specifications:
   Fs = 1250000;                   % samples per second
   dt = 1/Fs;                   % seconds per sample
   StopTime = 0.25;             % seconds
   t = (dt:dt:StopTime-dt);     % seconds
   %Input Sine wave:
   Fc = 300;                     % hertz
   dataIn = cos(2*pi*Fc*t);

   % Plot the signal versus time:
   figure;
   plot(t,dataIn);
   xlabel('time (in seconds)');
   title('Input signal');
   zoom xon;

   wc=500; % cut off frequency
   Wn=wc/(Fs/2); % normalized frequency
   [b,a] = butter(6,Wn,'low');
   figure;
   freqz(b,a)

%Plot filtered output
    dataOut = filter(b,a,dataIn);
    figure;
    plot(t,dataOut)
    title('dataout')
    zoom xon;

However I am unable to filter the desired frequencies below 500 Hz, as the cut off frequency is dependant on the sampling rate. Kindly suggest what to be done as I am new in MATLAB.

Thanks in advance.

Solution: Take a lower order butterworth filter or sampling frequency

Reason: As discussed here, a high order butterworth filter with a low (relative) cutoff frequency may be numerically unstable due to quantisation noise.

Explanation with figures:

I changed your frequency plot to include the region of interest:

freqz(b,a, logspace(1, 5, 1000), Fs)
ax = findall(gcf, 'Type', 'axes');
set(ax, 'XScale', 'log');

The frequency response is:

which seems very noisy at low frequencies, indicating an unstable filter.

This results in an extremely high output (blue = input, red = output):

If I change the butterworth order to 4, I get:

The output (red) is a shifted version of the input (blue) as can been expected from the frequency response.

Keeping a 6th order butterworth filter, but lowering the sampling frequency also solves your problem.