Zero padding - High amplitude I have a set of 512 samples that I have collected with a sensor. I'd like to apply zero padding to it, for better frequency bin resolution. However, all the tutorials I see to make this occur center around the idea that your wave is localized around zero. Here is an example:

Adding 10,000 zeros to the end of this signal seems logical to me.

However, say that the Y values on the above graph are shifted 5,000 points.

In that case, would you still add 10,000 zeros to the end of the singal?

Or, would you average out your signal ( find the median ), and use that as the padding? Given that the Y coordinates shifted 5,000 points, then we might be adding 10,000 pads of 5,005.

I'd like to apply zero padding to it, for better frequency bin resolution.

First of all, let's state it one more time that zero padding does not improve frequency resolution of DFT. It'll only interpolate the existing spectrum on a finer frequency grid, but will not add any new information to it, otherwise.

In order to improve the true frequency resolution of a DFT computation, you have to take more samples and therefore increase the length of the signal.

Zero padding is applied by adding zeros to the end of the signal to make it length N, where you then compute $N$-point DFT of it afterwards.

Lets consider the following concrete example of the sum of two sinusoidals: $$x[n] = 0.73 \cos(0.25\pi n +\pi/12) + 0.87 \cos(0.28 \pi n + \pi/17)$$

And we shall represent $x[n]$ by a finite length $L$ segment and take $N$-point DFT of it to see a spectral vision of it.

In the first case, let signal length $L=16$ and let DFT-length be $N = 64$. The following is the resulting DFT magnitude plot:

Now, apply zero padding to $x[n]$ and make it length $N=1024$ samples (by adding zeros) and take 1024-point DFT over the still $L=16$ length signal, with the result:

As you can see although increasing DFT length (equivalently zero padding) increases the smoothness of the spectrum, it does not resolve the two sinusoidals embedded in $x[n]$ at the window length of $L=16$.

Now lets increase the signal sample count by setting $L=512$ (taking more samples of $x[n]$ instead of adding zeros to it) and again taking $N=1024$ point DFT of the segment, whose result is:

As can be seen, now the two sinusoidals are finely resolved while still computing the same $N=1024$ point DFT but over a longer signal observation length ($L=512$), hence, demonstrating the fact that the spectral resolution is improved by longer signal observation and not by longer DFT (by zero padding).


What is the maximal frequency resolution for Matlab's STFT implementation spectrogram()? Matlab's spectrogram() function calculates the STFT of a signal. It describes its NFFT argument as follows:

  S = SPECTROGRAM(X,WINDOW,NOVERLAP,NFFT) specifies the number of
     frequency points used to calculate the discrete Fourier transforms.
     If NFFT is not specified, the default NFFT is used.

Am I correct in that NFFT is a trade-off only between frequency resolution and number of computations? For my offline work, there's no need to save cycles. Is there any maximum limit for NFFT, imposed e.g. by spectral leakage, or any other problem that I should know about, or can I set that argument to as high as possible?

The length of the FFT is a tradeoff between frequency and time resolution. Spectrograms are often generated by calculating overlapping FFTs on the signal of interest. If you make the FFT longer, then the effective bandwidth of each output bin becomes smaller, so the resolution along the frequency axis improves. The only limiting factor to the frequency resolution that you can obtain is the total observation time that you have from the signal.

At the same time, however, your ability to resolve a feature that is localized in time diminishes. An intuitive way to think about this is to view the FFT as a complex downconversion followed by an integrate-and-dump operation:

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

Viewing it in this way makes the loss of time resolution more apparent. The product in parentheses shifts $x[n]$ down in frequency by $\frac{2 \pi n k}{N}$, and the resulting signal is integrated over a window of $N$ samples. If there is a feature in $x[n]$ that is only located in a limited time span, then as $N$ becomes larger, more of the overlapped FFTs will contain that time period inside their integration time windows. Therefore, the feature will appear in more rows of the spectrogram image (assuming time is along the Y axis). If you then make a cut down the columns (i.e. the frequency bins) of the spectrogram that the feature is located in, you would notice a broader, smeared out peak. You therefore have a lesser ability to resolve the actual time location of the feature's onset.

You are also right that increasing the FFT length does require more computations, which might be relevant for real-time applications.


which of the analog sinusoidal frequency can not pass through the filter? First time I am encountering this type of question so i just tried but not getting whether my logic correct or not.

First let the sinusoidal signal be $X(t)=\cos(2\pi ft)$.

After sampling this signal I will get $X(nT_s)=\cos(2\pi fnT_s)$; now at sampling rate $f_s=18\, \text{kHz}$ it becomes $X(n)=\cos(\frac{n\pi f}{9})$.

Now I am checking optionS $(C)$ (given correct answer is also C) to see if the output of filter will zero or not at $f=12\,\text{kHz}$ and at  12kHz  X(n) becomes
$X(nT_s)=\cos\left(\frac{4\pi n}{3}\right)$  .

$H(z)=1+z^{-1}+z^{-2} $ can be written as
$H(e^{jw})=1+e^{-jw}+e^{-2jw}$, so from eigenvalue concept the value of $H(e^{jw})$ at $w=\frac{4\pi}{3}$ is $H(e^{jw})=1+e^{\frac{-j4\pi }{3}}+e^{\frac{-j8\pi }{3}}=0$, so it's proven that output at 12 kHz will be zero.

Therefore, 12 kHz frequency will not pass.

I am doing by checking options, option B also satisfying this condition.

Is there any alternate method to solve this?

The impulse response of the given transfer function $H(z)= 1 + z^{-1} + z^{-2}$ is $h[n] = \delta[n] + \delta[n-1] + \delta[n-2]$.

This is a rectangular window with three taps and its Frequency response magnitude is given by:

$$ H(w) = \frac{\sin( 3 w/2 )}{\sin(w/2)} $$

It'll have zeros (nulls) at the frequencies given by:

$$ \sin(3w/2) = 0 \longrightarrow 3 w_k / 2 = \pi k \longrightarrow w_k = 2\pi k / 3 ~~~~, \text{for} ~~~k = 0,1,2 $$

Note that the zero for $k=0$ is cancelled by the denominator and therefore there exists only two zeros $w_1 = 2\pi/3$ and $w_2 = 4\pi/3$

Now note that those two zeros are at mirror locations about $w=\pi$ which means that they will be the same, indistinguishable, discrete time frequencies.

Also you shall consider only the first zero , $w_1 = 2\pi/3$ to be the frequency of the analog sinusoidal that will be nulled by the filter.

The given sampling frequency of $f_s = 18$ kHz maps $w_1= 2\pi/3$ to the analog frequency of

$$ f_1 = \frac{f_s w_1}{2\pi}  = \frac{18k 2\pi/3}{2\pi} = 6 kHz$$

The answer therefore is $f_1 = 6$ kHz which cannot pass through your filter when sampled properly at the sampling rate of $f_s = 18 $ kHz. This assumes that there is no aliasing during the sampling process which is prevented by an analogue anti-alasing filter before sampling. If this filter is omitted and sampling is performed then a sinusoid at the frequency of $f_2 = 12$ kHz will be aliased into a frequency of $f_1 = 6$ kHz and therefore it won't also pass through the filter. Note that there will be infinetely many new frequencies which will map into $f_1 = 6$ kHz and won't pass through the filter in such a case.


discrete digital filtering in python I am not from DSP background and also not sure if I am posting this in right community. I have basic understanding of resistors, capacitors and filters.
My requirement is as follows-

The supply voltage is fixed and I would like to simulate the output of the filter using Python. With a fixed input voltage, I would like to know what is the output voltage at R3. The current is sampled at a frequency of 100 milliseconds. The input is a DC voltage, which can change.

I am confused where should I start from. Should I use scipy library?
Since I want a time domain output, what is the best approach?
Is there some specific reading material for such problem?
Which kind of filter would be this?

Please excuse for a basic question.

Obtaining the system equations

Using Kirchhoff's voltage law:
$$
\begin{align}
  V(t) &= V_{R_1}(t) + V_{R_2}(t) + V_{R_3}(t) + V_I(t) \tag{1}\\
  V_{R_1}(t) &= V_{C_1}(t) \tag{2} \\
  V_{R_2}(t) &= V_{C_2}(t) \tag{3}
\end{align}
$$
And similarly using Kirchhof's current law:
\begin{align}
  I_{R_3}(t) &= I_{R_1}(t) + I_{C_1}(t) \tag{4}\\
             &= I_{R_2}(t) + I_{C_2}(t) \tag{5}
\end{align}

Now using the Voltage-Current relationships of the resistances and capacitors (assuming ideal components):

\begin{align}
  V_{R_1}(t) &= R_1 I_{R_1}(t) \tag{6}\\
  V_{R_2}(t) &= R_2 I_{R_2}(t) \tag{7}\\
  V_{R_3}(t) &= R_3 I_{R_3}(t) \tag{8}\\
  I_{C_1}(t) &= C_1 \frac{d\,V_{C_1}(t)}{dt} \tag{9}\\
  I_{C_2}(t) &= C_2 \frac{d\,V_{C_2}(t)}{dt} \tag{10}
\end{align}

Based on your comments, a current sink with a resistive internal resistance can be assumed. For analysis, it is common to replace such a load with its Norton equivalent circuit (see Norton theorem) which consists of putting the internal resistance (say $R_4$) accross an ideal current sink $I(t)$.
\begin{align}
  V_I(t) &= R_4 I_{R_4}(t) \tag{11} \\
         &= R_4 (I_{R_3}(t) - I_I(t)) \tag{12}
\end{align}

Numerical solution

To discretize the problem to discrete time step $t_n = n \Delta t = n/f_s$, we can make approximations for $(9)$ and $(10)$:
\begin{align*}
  I_{C_1}[n]
    &\approx C_1 \frac{V_{C_1}[n+1] - V_{C_1}[n]}{\Delta t} \\
    &\approx C_1 f_s \left(V_{R_1}[n+1] - V_{R_1}[n]\right) \\
  I_{C_2}[n] &\approx C_2 \frac{V_{C_2}[n+1] - V_{C_2}[n]}{\Delta t} \\
    &\approx C_1 f_s \left(V_{R_1}[n+1] - V_{R_1}[n]\right)
\end{align*}
This gives us the recurrence equations:
\begin{align*}
  V_{R_1}[n+1]
    &\approx V_{R_1}[n] + \frac{1}{C_1 f_s} I_{C_1}[n] \\
    &\approx V_{R_1}[n] + \frac{1}{C_1 f_s} \left(I_{R_3}[n] - I_{R_1}\right) \\
    &\approx V_{R_1}[n] + \frac{1}{C_1 f_s} \left(\frac{V_{R_3}[n]}{R_3} - \frac{V_{R_1}[n]}{R_1}\right) \\
  V_{R_2}[n+1]
    &\approx V_{R_2}[n] + \frac{1}{C_2 f_s} I_{C_2}[n] \\
    &\approx V_{R_1}[n] + \frac{1}{C_2 f_s} \left(I_{R_3}[n] - I_{R_2}\right) \\
    &\approx V_{R_1}[n] + \frac{1}{C_2 f_s} \left(\frac{V_{R_3}[n]}{R_3} - \frac{V_{R_2}[n]}{R_2}\right)
\end{align*}
And $V_{R_3}[n]$ begin obtained by substituting $(8)$ and $(12)$ into $(1)$:
\begin{align*}
  V_{R_3}[n]
    &= \frac{V[n] - V_{R_1}[n] - V_{R_2}[n] + R_4 I[n]}{1+R_4/R_3}
\end{align*}

which can readily be implemented in python (with $V(t)$ and $I_I(t)$ being constant for illustration, but the extension to arbitrary function is trivial) with:

# numerical solution
T    = 1.0/fs
Tmax = 20
t = np.arange(0,Tmax,T)

vr3 = np.zeros(len(t))
vr1 = 0
vr2 = 0
for i in np.arange(len(t)-1):
    vr3[i] = (V-vr1-vr2+R4*I)/(1+R4/R3)
    i3  = vr3[i]/R3
    ic1 = i3 - vr1/R1
    ic2 = i3 - vr2/R2
    vr1 = vr1 + ic1/(C1*fs)
    vr2 = vr2 + ic2/(C2*fs)
vr3[-1] = (V-vr1-vr2+R4*I)/(1+R4/R3)

Steady-state solution

As a validation of the above numerical solution, we can check whether it converges to a steady state which can be computed analytically fairly easily (assuming a steady state can be reached, which typically requires $V(t)$ and $I_I(t)$ to be near constant). To do this we note that in the steady state derivative terms go to zero. As a result, $I_{C_1}(t) = I_{C_2}(t) = 0$ and consequently, $I_{R_1}(t)=I_{R_2}(t)=I_{R_3}(t)$. Thus $(1)$ can be written as:
\begin{align*}
  V(t) &= R_1 I_{R_1}(t) + R_2 I_{R_1}(t) + R_3 I_{R_1}(t) + R_4\left(I_{R_1}(t) - I_I(t)\right) \\
     &= \left(R_1+R_2+R_3+R_4\right) I_{R_1}(t) - R_4 I_I(t)
\end{align*}
From which we can solve for $V_{R_3}(t)$:
\begin{align*}
  V_{R_3}(t) &= \frac{R_3}{R_1+R_2+R_3+R_4} \left(V(t) + R_4 I_I(t)\right)
\end{align*}

For illustration, the graph below shows an example for the following parameters:

# parameters
R1 = 2200.0
C1 = 0.001
R2 = 1000.0
C2 = 0.0047
R3 = 5600.0
R4 = 100000.0
V  = 5
I  = 1e-5
fs = 10

As a final note, for voltage source $V(t)$ and current sink $I_I(t)$ which are constant or are from a set of fairly simple functions, the solution could also be obtained analytically with the use of Laplace transforms.


harmonic waves as integer multiple in spectrum i have a motor that is rotating with a certain frequency. If i check the frequency spectrum it contains a peak on 150 hz. Also i have peaks at 300,450,600 ... i guess that those peaks are harmonics.
But i dont really understand why harmonics appear.

I think the reason is:
a periodic signal can be described by sine, cosine waves (except the if periodic signal itself is a sine or cosine wave) So because a rotating motor is periodic, the frequency contains multiple cosine or sine waves that together build the signal. And the harmonics are the parts of the signal?

but I still don't understand why they occur if my machine is rotating at 150 hz. How are the other frequencies build?

Every periodic signal can be described as sum of "basis" functions. These basis functions are the sin and cosines at frequencies of integer multiples of a fundamental freqeuency. This fundamental freqeuency is a consequence of the periodicity of the signal. So if a signal is periodic with time $T$, then the fundamental frqeuency is $F_o = \frac{1}{T}$ Hz. That's the basic Fourier series theorem for periodic signals.

It is the waveform of the signal which will define which of these sines and cosines are present.

Consider the figure below:

It is not straightforward to look at the periodic signal above and tell that it will be comprised of so many frequencies. But it is, that's just a consequence that periodic signals are made of a sum of weighted sinusoids at integer multiples of a fundamental frequency.

Just as you require basis vectors in a vector space to define "any signal" in the vector space, the contributions of these basis vectors can vary depending on the vector being defined.Similarly periodic signals live in a vector space whose basis vectors are these sinusoids (that are defined by the periodicity) of the signal and hence there are contributions (big or small) from all these frequencies to define the periodic signal.

So your observations regarding the frquencies appearing in the spectrum are justified. Infact from the fundamental freuency of the motor (150Hz) we even know how fast the motor is rotating :)


Why are MFCCs of two equivalent signals completely different? I have an audio file, of which I calculate 16 (actually 15 because I omit the first one) MFCCs in R. When I stream this file via VLC-Player and Icecast2-Server, receive it in Java (with the Player of Javazoom-Lib) and then pass it on to R again for MFCC calculation I get completely different values? Has anyone an idea why this could be?

Additional info:

When I write the received data to a file again and view it next to the original file they look almost the same (Waveform + Spectrogram).
The file has a length of 3 seconds and contains the sound of a passing vehicle (car).
MFCCs are calculated for centered 44100 samples.
I need valid data in order to automatically classify cars and trucks with a SVM.
This classifier is trained with features calculated from 150 audio files (each
containing one vehicle).
Since R had problems with reading and processing mp3-files directly I first converted wav to mp3 and back again (in order to "simulate" the loss of information due to compression).

Furthermore, I automatically detect the vehicles in a continuous audio stream and therefore (for the particular example data above) the waveforms for calculating the MFCCs have an offset of about 4400 Samples. Does this matter with an analysis window of 44100 samples?

You say:

  I stream this file via VLC-Player and Icecast2-Server, receive it in Java (with the Player of Javazoom-Lib)

When you receive the stream in Java, what is the stream format?

You say that it is the "same" because the waveform and/or spectrogram look similar, but the MFCCs will come out different if the stream format (sample rate, bit depth, etc) are different.

Can you confirm by printing out the stream formats

Before sending to Java
Upon receiving from Java?

I recommend using Soxi to print out the stream format (if you cannot otherwise dump it from R, or Java)

If you can eliminate stream format issues, mp3<->wav issues (deal only in wav), if you analysis windows are identical (as jojek says), then your MFCCs should come out identical.


How to get rid of discontinuities at the start of a signal caused by the filter I'm processing an audio signal using an overlap add method.Everything at the block boundaries is fine. However I have a discontinuity at the very start of the first block and I'm not sure how to (or if I can) get rid of this.Zero padding the start seems to make no difference as the discontinuity is obviously "triggered" when the filter first encounters the non-zero signal.

If you are post processing you can instead use a "zero-phase" filter which passes the signal through the filter twice in the forward and reverse direction, resulting in an output that is completely aligned with the input as would be expected with a filter with no delay, and with a significant reduction in start-up deviations (this is a non-causal filter hence the restriction on postprocessing).
Such a filter is available in MATLAB, Octave and Python scipy.signal as filtfilt which has further documentation on its use.
I demonstrate this as shown with the plot below using a traditional 91 tap FIR low pass filter for an in-band cosine signal and again with the same coefficients with filfilt (in red). Notice there are no transients or delay.

As far as any residual transient as @MBaz astutely points out, below is the plot zoomed in on the first 4 samples showing the filtered output as (+) and the original signal as (o). The line plot simply connects the output samples. This is due to the fractional cycle remaining at the end of the sampled data where I did worst case with a 1/2 cycle residual in the capture; if done with an integer number of cycles there would be no such mismatch.

And the same portion of the waveform later in the sequence showing close alignment.

However in comparison to the blue line in the first plot even in this worst case condition, we do see the significant reduction in the start-up transient due to the use of the zero-phase filter (limited to post-processing applications). The alignment of the input and output is also a huge plus for analysis applications where we many otherwise have to go through additional efforts to offset our waveforms for the delay when comparing input to output.


Algorithm to effectively find perspective transform in motion estimation I've shot a series of photographs in a sort of stop-motion fashion, and although the camera was mounted on a tripod, there are tiny motions between each two successive photos.

After some experimentation, I figured that using a perspective transform would be appropriate to apply a correction to a photo so that the difference frame is minimised. This transform basically just moves the four corners of a picture by a given offset, allowing for a trapezoid distortion.

To do this process by hand is impossible in this case, as it takes me like half an hour for each frame. So I am looking for an effective algorithm to estimate the transform parameters (the four corner offsets). I would give the algorithm the maximum magnitude of the offset, and it should then find, as quickly as possible, the parameters (eight values).

Is such an efficient algorithm known? The photos are in colour but could be converted to grey scale, although the lighting changes also over time, so taking a difference image and summing the pixels is not enough, but probably one needs to auto-adjust the contrast and brightness as well.

The suggested technique seems to be to gather points of interest, calculate a feature vector for them, and then find matching pairs of points in both images, followed by a calculation of the transform.

Here is a sample project: http://www.codeproject.com/Articles/95453/Automatic-Image-Stitching-with-Accord-NET

Here is a paper that claims improvement by preparing through a phase correlation stage: http://www3.nd.edu/~kwb/ThomasKareemBowyerIGARSS_2012.pdf


Tradeoff between Filter Order and Sections I have designed a bandpass IIR butterworth filter with desired response using the MATLAB FDAtool with the configuration shown in the figure below

Since this is multistage implementation, when I export the filter object to workspace I don't get numerator and denominator of transfer function as filter parameters and instead I get sos matrix and scale values as shown below

FilterStructure: 'Direct-Form II, Second-Order Sections'
              Arithmetic: 'double'
               sosMatrix: [12x6 double]
             ScaleValues: [13x1 double]
     OptimizeScaleValues: true
        PersistentMemory: false

The problem is I have to port this design to C code for implementation in firmware and the multi-order specification looks complicated to define an input-output expression. So I tried converting this design into a single stage filter(using convert to single section option under edit), but the following is magnitude response that what I get.

Not only is the response way-off, but the system is unstable.

My questions are the following,

1) Why is FDAtool not able to produce a good single stage equivalent filter?

2) How to calculate groupdelay of the multi-stage filter?

3) Is the groupdelay of the multistage filter simply FilterOrder*Sample time as in a direct form FIR or is it more complicated for IIR and multistage filters?

4) Will the groupdelay of multi-stage filter be the same as a possible acceptable(in terms of magnitude response) single stage equivalent?

5) How important is stability of the filter?

6) Does FDAtool have an option to visualize the filter structure?(Edit>Show Filter Structure doesn't seem to work)

7) What is the tradeoff between Filter Order and Number of sections in the filter other than structural complexity?

Large order single stage filters are inherently problematic. Calculating the numerators and denominator from the poles and zeros is an ill defined problem and numerically unstable.
Furthermore Direct Form II is the worst possible filter topology. I strongly recommend Direct Form I or transposed Form II.

To the questions:

Numerically ill defined problem
Just add the group delay of the indivdiual stages
Group delay is always defined as the negative first derivative of the phase of the transfer function vs. frequency. That's the same for FIR and IIR. It's frequency dependent.
Yes. If you had infinite precision math, the results would be the same.
Very important. You can't implement and unstable filter, you'd only get noise NaN or infinity as output
No idea
Pretty much everyone implements filter ins biquad sections. So the number of stages and the filter order is directly proportional


Difference between applying a window to a whole signal and frames. I'm going to do feature extraction using a speech signal.

The first step is the windowing of the signal using Hamming window.

I'd like to know if I should apply one single window to the whole signal or if I should apply several windows to different frames of the signal?

What are the differences between them?

Thank you.

Each frame to be FFT’d needs to be individually windowed to avoid rectangular window artifacts (sometimes called “leakage”) due to the finite length of the FFT.

The frame length you choose depends on the time-frequency resolution trade-off needed to extract the features in which you are interested.


Accelerometer Inference from power spectrum I have 3d accelerometer data sampled as 60hz. When I am looking at the power spectrum in all three accelerations (x,y,z). I am noticing a trend with a negative slope. I trying to figure a few things out:

What could the cause of this trend be? System bias or something else?
Best way of handling this - Matlabs detrend tool does not have much of an effect. What would be the best way to correct this?

What could the cause of this trend be? System bias or something else?

Most mechanical systems are low-pass, i.e. they let through less energy at higher frequencies. That would be perfectly plausible, but I don't know what you're observing.
Of course, your accelerometer is inherently frequency-selective, too, and you might need to calibrate it, together with your recording setup.
Lastly, we don't know exactly how you post-process.

Best way of handling this - Matlabs detrend tool does not have much of an effect. What would be the best way to correct this?

"Correcting  this" means "I saw something that is clearly visible in my measurement, but I don't know where it comes from and whether it's actually there in the observed phenomenon. Let me get rid of the evidence".
That's called a lie.
So, before changing anything about the symptom (plot), make sure you're not masking a cause. This is less about the signal processing than about physically understanding your measurement setup. It can, and probably will, involve you measuring known acceleration amplitudes at different frequencies, and looking at whether the results are as expected, or whether you need multiply with a different factor for different frequencies.


Modulation and DeModulation I have a question regarding Modulation and Demodulation.
In question 5.21 (b) from the "Discrete-time Signal Processing" book by Oppennheim, he asks to find out the type of filter used in the attached figure:

The question is what is the effect of $(-1)^n$ on the original filter and y[n].
I saw the answer to be high-pass filter, but I am not sure why.
Please help

Thanks

As hinted out, you firstly shift signal by π by multiplying it by cos(πn). Now, your low frequency components are around π. In the next step of diagram, you filter the components around frequency 0 since the filter in diagram is lowpass. But since your signal was shifted by π, you filter the high frequency components of the original signal. In the next step, you shift your original signal back to 0 frequency.

So, overall you filter the high frequency components.


Shannon theorem and zero padding Sorry if I ask a basic question but I am a bit lost.

I have a data set in frequential domain, I need to go to time domain to perform time-gating and then going back to frequencies. In order to do this, I compute an inverse Discrete Fourier Transform.

Shannon theorem tells us that two times the difference between the greatest and the least frequency, respectively noted as $f_N$ and $f_0 ( \ne 0)$, contained in the signal must be less then the sample-rate in time domain ($F_t$). Do I need to do some zero padding to avoid aliasing ?

Without zero padding on frequential domain, what is the relation between $F_t$ and $f_N, f_0$ ? My guess is that $F_t=f_N-f_0$. Is that correct ?

If I must do some zero padding, should I do it this way :

$$(0,\dots,0,f_0,\dots,f_N,0 \dots,0)$$

With enough zeros such that the first frequency contained in the padded signal is $0$ and the last is $2 f_N$ ?
If my guess above is correct, then the new sample rate is high enough to avoid aliasing, because it will be $F_t=2f_N>2(f_N-f_0)$.

Thanks for your help.

You'll not experience aliasing – in the end, the IDFT result per definition only contains frequency components that can be represented by the original frequency domain.

If you will, just stop worrying about aliasing: the DFT and inverse DFT are just base change matrices in $\mathbb C^N$, and are invertible.

Aliasing can only happen when something is non-invertible. This is not the case here; so, no problem.

What you might see is leakage, but you're already taking care of that by windowing.


What is the Effect of Multiplying a Function by the Unit Impulse Function in the Frequency Domain? I know about the the shifting property of the impulse function in the time domain as can be seen in equation $(1)$.

$$
\int_{-\infty}^{\infty} f(x)\delta(x - a)dx = f(a)\tag{1}
$$

But what is the effect of multiplication of a function by the impulse function in the Frequency domain? I.e

$$
X(\omega) = \delta(\omega - \omega_0)\cdot H(\omega)
$$

I think there is a slight typo in robert bristow-johnson's answer. Should be

\begin{align}
x(t) &= \frac{1}{2 \pi} \int\limits_{-\infty}^{\infty} X(\omega) e^{j \omega t} \, \mathrm{d}\omega\\
\\
 &= \frac{1}{2 \pi} \int\limits_{-\infty}^{\infty} H(\omega_0) \delta(\omega - \omega_0) e^{j \omega t} \, \mathrm{d}\omega\\
\\
 &= \tfrac{1}{2 \pi} H(\omega_0) e^{j \omega_0 t}
\end{align}

What this actually means is that during the inverse transform, $\delta(\omega - \omega_0)$ selects point frequency $\omega_0$ from the frequency response $H(\omega_0)$, and the value of $H(\omega_0)$ at that point determines the magnitude and phase of the time domain sinusoid $e^{j \omega_0 t}$.


Practical advices for applying phase correlation in image registration task I'm using OpenCV to detect shift between 2 images, here is sample code (based on cv::phaseCorrelate function):

cv::Mat ShiftTransform(const cv::Mat &source, const cv::Mat &target, bool bUseHanningWindow)
{
    CV_Assert(source.type() == CV_32FC1);
    CV_Assert(target.type() == CV_32FC1);

    Point2d shift;
    if(bUseHanningWindow)
    {
       Mat hann;
       createHanningWindow(hann, source.size(), CV_32F);
       shift= phaseCorrelate(source, target, hann);
    }
    else
    {
        shift= phaseCorrelate(source, target);
    }

    cout << "Detected shift: " << shift << endl;

    Mat H = (Mat_<float>(2, 3) << 1.0, 0.0, shift.x, 0.0, 1.0, shift.y);

    Mat res;
    warpAffine(source, res, H, target.size());

    CV_Assert(res.size() == target.size());
    CV_Assert(res.type() == CV_32FC1);

    return res;
}

For example images from wikipedia seems it works.

source

target

result

response map

But for some real life images it gives wrong shift, maybe because of this problems described in wikipedia:

  In practice, it is more likely that   $g_{b}$ will
  be a simple linear shift of   $g_{a}$ , rather
  than a circular shift as required by the explanation above. In such
  cases,  $r$ will not be a simple delta function,
  which will reduce the performance of the method. In such cases, a
  window function (such as a Gaussian or Tukey window) should be
  employed during the Fourier transform to reduce edge effects, or the
  images should be zero padded so that the edge effects can be ignored.
  If the images consist of a flat background, with all detail situated
  away from the edges, then a linear shift will be equivalent to a
  circular shift, and the above derivation will hold exactly. The peak
  can be sharpened by using edge or vector correlation.

  For periodic images (such as a chessboard), phase correlation may
  yield ambiguous results with several peaks in the resulting output.

But I'm not sure how to debug such cases if something goes wrong.

Here is suggested not to blur images.

What are other practical advices for applying phase correlation in image registration task?

The effect of circularity could be decreased by bigger border filled up with zeros (up to width of image). But the method is not perfect inherently, it works perfectly just for circular shift of image. Otherwise the method can fail, mainly when one part of image is very light (high values of pixels). The rest of image should have as uniform distribution of light intensity as it is possible.


Butterworth High Pass Filter Amplitude Response I have learned the Butterworth filter, normally it is used for low pass design. And the amplitude response is:

$$|H(jw)|=\frac{1}{\sqrt{1+(\frac{w}{w_{c}})^{2n}}}$$

mentioned in this Butterworth Introduction.
But as I know, the Butterworth can be designed as a high pass as well. Then what is the amplitude response for high pass? Is it simply like this?

$$|H_{hp}| = \frac{1}{|H(jw)|}$$

If $H_{LP}(j\omega)$ is the frequency response of a Butterworth low-pass filter, then its inverse would go to infinity for $\omega\rightarrow\infty$ (because $H_{LP}(j\omega)$ goes to zero for high frequencies). So that's not what you want.

The correct high-pass frequency response is obtained by replacing $\omega/\omega_c$ by $-\omega_c/\omega$, resulting in

$$|H_{HP}(j\omega)|=\frac{\left|\frac{\omega}{\omega_c}\right|^n}{\sqrt{1+\left(\frac{\omega}{\omega_c}\right)^{2n}}}$$

which approaches unity with increasing frequency, as should be the case for a high-pass filter.


generating log mel spectrogram using librosa I have come across 2 different ways of generating log-mel spectrograms for audio files using librosa and I don't know why they differ in the final output, which one is "correct" or how different is one from the other.
#1
path = "path/to/my/file"
scale, sr = librosa.load(path)
mel_spectrogram = librosa.feature.melspectrogram(scale, sr, n_fft=2048, hop_length=512, n_mels=10, fmax=8000)
log_mel_spectrogram = librosa.power_to_db(mel_spectrogram)
librosa.display.specshow(log_mel_spectrogram, x_axis="time", y_axis="mel", sr=sr)

#2
path = "path/to/my/file"
scale, sr = librosa.load(path)
X = librosa.stft(scale)
Xdb = librosa.amplitude_to_db(abs(X))
librosa.display.specshow(Xdb, sr=sr, x_axis='time', y_axis='hz')

The respective images are:

** EDIT ** Now that I specify the number of mel bins to be = 64, I obtain the spectrogram as below:

If I want to process many such spectrograms, should I trim off the bold blue portion above as it is common to all? What does the bold, dark region represent? Is it advisable to use fmax parameter to trim it?

The mel spectrogram additionally includes a step of projecting (power of) STFT bins onto Mel-frequency bins via a Mel filterbank; I don't have access to path so I made demo on exponential chirp:

You can visualize the kind of projection taking place by plotting the mel basis:

Note in general the two won't look alike unless filters.mel are carefully selected (nor do they have to).

Code
import numpy as np
import librosa
import librosa.display
import matplotlib.pyplot as plt

def echirp(N, fmin=.1, fmax=None, tmin=0, tmax=1):
    """https://overlordgolddragon.github.io/test-signals/ (bottom)"""
    fmax = fmax or N // 2
    t = np.linspace(tmin, tmax, N)

    a = (fmin**tmax / fmax**tmin) ** (1/(tmax - tmin))
    b = fmax**(1/tmax) * (1/a)**(1/tmax)

    phi = 2*np.pi * (a/np.log(b)) * (b**t - b**tmin)
    return np.cos(phi)

#%%
N = 2**16
n_fft, hop_length, n_mels, sr = 2048, 512, 128, 16384
x = echirp(N)

#%%
mel_spectrogram = librosa.feature.melspectrogram(
    x, sr, n_fft=n_fft, hop_length=hop_length, n_mels=n_mels)
log_mel_spectrogram = librosa.power_to_db(mel_spectrogram)
librosa.display.specshow(log_mel_spectrogram, x_axis="time", y_axis="mel", sr=sr)
plt.show()

#%%
X = librosa.stft(x, hop_length=hop_length)
mel_basis = librosa.filters.mel(sr, n_fft, n_mels=n_mels)
Xmel = np.dot(mel_basis, abs(X)**2)
Xdb = librosa.power_to_db(Xmel)
librosa.display.specshow(Xdb, x_axis='time', y_axis='hz', sr=sr)
plt.show()

#%%
mel_basis = librosa.filters.mel(sr, n_fft, n_mels=10)
plt.plot(mel_basis.T, color='tab:blue')
kw = dict(weight='bold', fontsize=16, loc='left')
plt.title("n_mels=10", **kw)
plt.show()

mel_basis = librosa.filters.mel(sr, n_fft, n_mels=n_mels)
plt.plot(mel_basis.T, color='tab:blue')
kw = dict(weight='bold', fontsize=16, loc='left')
plt.title("n_mels=%s" % n_mels, **kw)
plt.show()


What is the frequency response function for the Bode Plot pictured here? I've forgotten how to get the frequency response function out of a Bode Plot. The phase component of $H(j\omega) = e^{j\pi/2}$, right? But I am having trouble finding $|H(j\omega)|$. How would I get it from decibels?

Your Bode plot is a straight line. As we are on a semilog-scale, the equation is given by

$$y = m \log(x) + b$$

In this case, $y$ would be the magnitude of the frequency response $|H(j\omega)|$ in $\mathrm{dB}$, and $x$ would be the frequency $\omega$. Taking two points of the plot and solving the $2\times2$ equation system should give you the function you're looking for.

As you didn't put the values of the $\omega$-axis, I'm going to assume that the plot starts at $\omega=1$. In this case, when $\omega = 1$ the magnitude is $0 \ \mathrm{dB}$, so:

$$0 = m \log(1) + b\implies b=0$$

Also, when $\omega = 10$, the magnitude is $20 \ \mathrm{dB}$:

$$20 = m \log(10) +0\implies m=20$$

So the function is:

$$|H(j\omega)|_{\mathrm{dB}}=20\log(\omega)$$

Note that if your plot doesn't begin at $\omega=1$ (what I just assumed to be true), the reasoning would be the same but $b$ would not be zero.


Extracting audio from picture of time-domain signal I saw this nifty art thing where they dress up the time-domain audio signal in some pretty way.
https://www.etsy.com/ArtbloxShop/listing/621561709/sound-wave-art-custom-soundwave-print

My question is whether I can take a photograph of that object and somehow reconstruct the actual audio?  I'm assuming that the printing and my photograph won't have the resolution to directly capture the wave data points.  But I was wondering if using some basic assumptions about frequency if we could reconstruct something crude.  I suspect that too much information is lost, but I was hoping there might be some clever tricks.

There is clearly some resolution of printing that would preserve enough information, but I'm wondering if we're talking about a poster-sized version of a short clip or something.  That wouldn't be very practical.  It would be fun to have time-domain audio hanging as art if I could take my smart phone and take a picture of it, and after some processing, play back the audio.

As an alternative, could this be done using a spectrogram?  Like if that block of art displayed that instead?  Or maybe some combination of that plus time-domain envelope?  If I did want to represent audio as an image as part of art (maybe presume something around 8.5x11 size), is there some way I could have an app photograph it and turn it back to audio?

In general no, there is way too much information lost

If you can constrain the audio to a reasonably small corpus of songs, phrase, keywords, etc. you may be able to train a neural network to recognize the basic envelope shape. There used to be an app that was trying to do something like this, but I can't find it anymore.


signal uncorrelated with delta in the origin I have a stochastic process completely uncorrelated. Why the autocorrelation function has a delta of dirac in the origin?
Which is the reason of that?

$R_{XX}({\tau})=A{\delta(\tau)}$

where $R_{XX}$ is the autocorrelation function of $X$ process and $A$ is a constant dependent by process spectrum

It has a delta at the origin because you are multiplying your time series with itself (i.e. there is no time delay).  Any non-zero time series correlated with itself will give a positive value.


Redefine Median Filter in Matlab How can I redefine median filter function instead of using the available function in Matlab? I want to smoothen my sound chart

You can shadow a built-in function as long as the newly defined function is higher on the function precedence order.

Note, however, that this is not recommended because:

It makes your code harder for someone else to understand
It makes your code harder to debug if something goes wrong

and plenty of other reasons.

As the commenter noted, it is usually better to just define a new function that does what you want with a different name.

EDIT:

If you are looking for specific ways to do something similar to a median filter, one place to look would be pseudomedian filters.


which are saturation and intensity images in this figure? The attached article talks about hue,saturation and intensity images. Where is each of them in fig6.25 ? I know that fig(b) shows Hue image but where are other two(saturation and intensity images)

When the text refers to Hue component image, it is to represent what hue value these colors hold, not an actual picture in the text. When the text says H is shown in figure b, it means the numerical representation of these three colors (or hue) are illustrated more clearly when we take away saturation.

That being said, the saturation for all four figures (a to d) is being held constant as well. Figures a and d having high saturation, while figures b and c have no saturation. So I wouldn't say any of these figures are a Saturation image, by the definition you might be looking for. To make that picture you would need a square of a single color (Hue) and then change the saturation of each of the four corner squares.

Likewise, all of these figures have constant intensity in each corner. So I wouldn't say any of these are Intensity images either.

All of these figures have "component images" in them: H, S, and I.

Figure a and b illustrate the effects of changing Hue on an image with static saturation and intensity, while figure c and d show what happens when you "blur" (taking color averages over small portions of the image) the first two figures.


Is there a function in python which returns the Amplitude/Sound Pressure Level of a sound file(like .wav) file? I found a function called scipy.io.wavefile.read() which returns the sampling rate and an array. But in the documentation, it is not specified that whether the array returned by this function is of Amplitude or Sound Pressure Level. Which one is it?

No, you need to know the full Transfer function from the acoustic response at the sensors face to your A/D converter.  Transducer responses vary with frequency. You need a calibrated reference at a minimum.

A reference sensor will have a response curve from the manufacturer.  An average would be useful as well.


Power Spectral Density vs Energy Spectral Density for Discrete Signal I'm sure this is a silly question but reading a bit about the PSD vs ESD for continuous signals, I'm a bit confused how that applies to discrete signals.

Usually people seem to be using the PSD when analysing real world signals, why not the ESD, I would assume a discrete real world signal is never infinite and thus never encounters the problem of infinite energy?

Energy spectral density is mostly applied to time-limited signals. In such cases, energy is finite, and average power (throughout all time) is zero. These signals are not stationary.

Normally you work with (or model your signals as) stationary signals. These have constant power (and infinite energy). Power spectral density applies to stationary signals, and is defined as the Fourier transform of the autocorrelation function.


Origin of term spectrogram I'm writing up a manuscript and I use the term spectrogram quite a bit.  I've defined a spectrogram as the magnitude squared of the short time Fourier transform of my signal.  It's an estimator for a time-varying spectrum of a non-stationary process.  I'd like to cite a paper for the origin of the term but I'm not sure what to cite.  Could anyone give me some advice?

IIRC, etymology of "Spectro" is "of radiant energy". "Gram" is derived from the Greek word gramma, γράμμα which is a written learning.


Fluctuation of autocorrelation of a signal due to signal's noise I have a question about the fluctuation of autocorrelation of a signal due to signal's noise.

I have a signal defined in $-1\leq t \leq 1$ as the following:
$V(t)=kt+R(t)$, where $R(t)$ is the random noise, and $k$ is non-negative.

Then, when I calculate the autocorrelation of $V(t)$, $G(\tau)=\int_{-1}^{1-\tau} V(t)V(t+\tau)dt$ (using the commend CorrelationFunction of Mathematica), I could observe that the fluctuation of $G(\tau)$ reduces as $k$ increases from zero. Could you explain the origin of this observation? Is it because as $k$ becomes larger, there is more autocorrelation and $R(t)$'s contribution to $G(\tau)$ is diminished?

If you view $V(t)$ as desired signal $kt$  affected by noise $R(t)$, the clearly as you keep increasing $k$, this has the effect of increasing signal energy. Hence if the noise power is kept constant then $V(t)$ tends towards being more deterministic (in a loose sense) with improving SNR, and autocorrelation hence is more deterministic with reduced fluctuations.


Demonstrate that $x(t)$ and $\hat{x}(t)$ are orthogonal I have to demonstrate that a function $x(t)$ and his Hilbert's transform are orthogonal, it is said:

$$\int^{\infty}_{-\infty} x(t) \cdot \hat{x}(t) dt = 0$$

I have tried the exercise using Parseval's theorem:

$$\int^{\infty}_{-\infty} x^2(t) \cdot \hat{x}^2(t) dt = \int^{\infty}_{-\infty} |X(f) \ast \hat{X}(f)|^2 df = \int^{\infty}_{-\infty} |X(f) \ast (-jsgn(f)X(f))|^2 df =$$

$$= \int^{0}_{-\infty} |X(f) \ast jX(f)|^2 df + \int^{\infty}_{0} |X(f) \ast (-jX(f))|^2 df$$

Using properties of convoution ($ax(t) \ast y(t)) = a(x(t) \ast y(t))$:
$$\int^{\infty}_{-\infty} x^2(t) \cdot \hat{x}^2(t) dt = \int^{0}_{-\infty} | j(X(f) \ast X(f))|^2 df + \int^{\infty}_{0} |-j(X(f) \ast X(f))|^2 df =$$ $$=\int^{0}_{-\infty} |X(f) \ast X(f)|^2 df + \int^{\infty}_{0} |X(f) \ast X(f)|^2 df = $$
$$= -\int^{\infty}_{0} |X(f) \ast X(f)|^2 df + \int^{\infty}_{0} |X(f) \ast X(f)|^2 df = 0$$

But I don't know how to relationate both integrals. Maybe that isn't the way to demonstrate the exercise.

Thanks for your help.

You have to use the following form of Parseval's identity:

$$\int_{-\infty}^{\infty}x(t)y^*(t)dt=\int_{-\infty}^{\infty}X(f)Y^*(f)df\tag{1}$$

where $^*$ denotes complex conjugation, and where $X(f)$ and $Y(f)$ are the Fourier transforms of $x(t)$ and $y(t)$, respectively. With $y(t)=\hat{x}(t)$ and $Y(f)=-j\,\textrm{sgn}(f)X(f)$, we get from $(1)$

$$\int_{-\infty}^{\infty}x(t)\hat{x}^*(t)dt=j\int_{-\infty}^{\infty}|X(f)|^2\textrm{sgn}(f)df\tag{2}$$

For real-valued $x(t)$, $|X(f)|^2$ is an even function, and, consequently, the integrand on the right-hand side of $(2)$ is odd. This implies that the integral on the right-hand side of $(2)$ equals zero if it is interpreted as a Cauchy principal value.

Note that we've used the fact that $x(t)$ is real-valued. For complex-valued $x(t)$, the integrals in $(2)$ are generally not zero.


Time stretching allowing for pitch change I have implemented the WSOLA time stretching algorithm for audio data, but there are still many artifacts from the stretching.

I don't actually care about maintaining pitch and am totally ok with the "chipmunking" effect. What options are available for time stretching in this case? Are there any algorithms I can use that will produce fewer artifacts?

I've only been able to find questions that DO care about maintaining pitch. I've read about resampling, but don't understand exactly how that would work when not trying to double or halve the speed of the sample.

Any information or links are appreciated!

Resampling does appear to be the answer for the OP. The point is that you will not change the speed of the sample as you will not change the sampling rate by resampling just the number of samples, so therefore the time represented by each sample must change after resampling. An example will best clarify this :

We have 10 seconds of audio data sampled at 44.1 KHz.

This is a total of 441,000 samples over 10 seconds.

We resample to twice the number of samples by inserting zeros between each sample and then low pass filtering out the images created (standard interpolation techniques that will minimize distortion).

Now we have 882,000 samples.

The thing is we don't actually change the sampling rate! So therefore time must have stretched.

882,000 samples / 44,100 samples/sec =  20 seconds duration!

This will result in a time stretching changing all frequencies accordingly.

Most resampling techniques are described as changing the sampling rate but would be equally applicable to this. So if the resampling technique describes changing the sampling rate by $M$, we can equally run the new samples at the same sampling rate and time would have changed by the same $M$.


How to remove or filter the drift problem in measured Strain signal? I have the strain signal of a lateral beam of a car measured at sampling rate 1,200Hz from data acquiring system. Even after using temperature compensation in strain gage, we are getting drift. So I wanted to remove drift in post processing of the strain signal. I have tried using High pass filter (HPF) of order 10th with cutoff freq = 200Hz, but I am losing peaks of the strain signal. So HPF is not useful.

Which filter is suitable to remove drift for strain signals? considering the fact that I want to retain peaks  and what other methods I can use to remove the drifts?

Also, what are the exact factors to consider in choosing the cutoff frequency for HPF? and how to check it is correct?

I have attached the strain signal and HPF filtered signal of 10th order.

The usual first approach to this is to use a DC Blocker:
$$
y[n] = \alpha y[n-1] + x[n] - x[n-1]
$$
where $0 < \alpha < 1$, $y$ is the DC blocked signal and $x$ is your original signal.

If I simulate your signal and apply the DC Blocker to it, the results are in the figure below.

It should preserve the peaks well.

R code to implement it below.

How to choose $\alpha$ ?

This really depends on what you are going to use the data for later.  If it's a control application, then phase is probably important to you.

More generally, you want to JUST remove the DC component and not much else.  In that case, just look at the magnitude frequency response of the filter.

As you can see from the plots below, low values of $\alpha$ (less than $0.9$ tend to have a bigger impact on the lower frequencies.  Hence, I tend to use values of $\alpha$ between $0.9$ and $1.0$.

R Code Below

#27468

T <- 1000
Npeaks <- 10

idx_peaks <- runif(Npeaks,1,T)

strain <- rnorm(T) + c(seq(1,2*T/3,1) ,seq(2*T/3+1,0,-2))/100

strain[idx_peaks] <- strain[idx_peaks] + 10

detrended <- 0*strain
alpha <- 0.9

for (k in 1:length(strain))
{
  if (k>1)
  {
    detrended[k] <- alpha*detrended[k-1] + strain[k] - strain[k-1]
  }
  else
  {
    detrended[k] <- 0
  }
}

par(mfrow=c(2,1))
plot(strain,col="blue", type="l")
plot(detrended,col="blue", type="l")


Loop Bandwidth for Carrier Demodulation loop I am trying to implement a simple MATLAB model for the demodulators for phase modulated signals (PSK) used in communication systems.

For the detection of the phase modulated output (PM out in the above figure) the bandwidth of the demodulator should be sufficiently lower than the symbol rate (modulation bandwidth) so as to ensure that the PLL is not able to track the modulation. But clearly, low frequency phase noise will be tracked which is the benefit of this demodulator.
Since the demodulator is working as a high pass filter, I assume this can be implemented with just a high pass filter in MATLAB.
I generate the input data as follows:

data = randi([0 7], 1e6,1);
modData = step(pskModulator, data);
scatterplot(modData);

I assume that my symbol rate is 1GS/s and since we want the loop bandwidth to be much lower, I keep it at 100kHz. I generate a high-pass filter as follows:

[b,a] = butter(1,100e3/1e9,'high');
H = tf(b, a, 1/1e9);
freqz(b,a)

Clearly close to 1GHz the filter reaches its pass-band.
MY QUESTION
My issue is that if I pass my input data through this filter, it still gets filtered out. I check it as follows:

angle_out = arrayfun(@(x) angle(x), modData);
angle_out_loop = filter(b,a,angle_out);
out = exp(1i*angle_out_loop);
scatterplot(out)

What am I doing wrong here? How can I make a filter which filters the low frequency noise but not the signal?
Although I have assumed 1GS/s as my symbol rate but how does MATLAB know that the random input data is indeed coming at 1GHz so that the data doesn't get filtered?
I have read similar qustions:Loop Bandwidth for symbol timing recovery but they were less implementation specific.
Thanks!

Your results are not showing that the data is being "filtered out" necessarily but that you are using a non-linear phase filter resulting in group delay variation over your passband, as evidenced by the spread in phase.

Don't use a Butterworth filter to model this. Instead use filters derived from either firls or firpm (remez in Octave) which are linear phase and have constant group delay.


Do Audio Interface manufacturers use A-weighted Measurements to make specs look better? Recently I came across a reputed music-tech teacher who said audio interface manufacturers use A-weighted measurements to show dynamic range in their specs so that they look better on spec sheets and to hide noise in some frequencies. He also said the A-weighted measurements are not used to measure audio interfaces and advised against buying any audio interface that has A-weighted measurements in specifications. Then I was looking at some really high end audio interfaces and noticed they're also displaying their specs in A-weighted measurements.

So what's the real deal about this? Is it fair to say an audio interface is bad just because A-weighted measurements are used in specs?

Most audio interface devices are built for the commercial electronics industry which produces devices to be used by people to listen to music, radio, TV etc; for multimedia reproduction purposes.

Therefore the commercial electronics standardisation organisations, suggest or enforce the use of a number of frequency weighting filters (A,B,C etc.) for getting more relevant results for the technical specs involved, considering that human beings will be the final consumer of the audio end result. It's not applied for getting an illusive improvement of the numbers (mostly SNR), but for providing an acceptable performance at the lowest cost.

If your purpose of audio (acoustical) interaction is some engineering or scientific measurement and computations, then look for the unfiltered (bare) specs if available. You may ask manufacturers to provide you those numbers, if possible. But most would restrain from disclosing that information due to confidential reasons...


FIR filter not eliminating out of band noise? I have a simple FIR filter that should eliminate all lower frequencies:

H = zeros(60,1);
H(30:end) = 1+1i;
H = [H;conj(flipud(H(2:end)))]; % double sided
h= ifft(H);

Then I assign a sample time to generate a time vector:

dT = 0.002;
Tspan = dT*length(h);
dF = 1/Tspan;

Here my nyquist is 250Hz.  So the filter should eliminate everything below 125Hz.

I generate a signal at 15Hz, and convolve it with the filter coefficients.

t = 0:dT:5;
x = cos(2*pi*15*t);

But, I get an output signal that clearly does not have the 15Hz eliminated:

Input signal above, response signal below

Below is the FIR filter (plotted single sided):

Does it look like my convolution algorithm is incorrect? Or am I misunderstanding something?

Any help is appreciated!

Thanks

You can solve your problem by an fftshift of the impulse response, as this line shows:

h = fftshift( real( ifft(H)) );

now it should work.

NOTE that there will be a transient of length about $60$ samples (half the linear phase FIR filter length) at the beginning of the convolution output, after which the steady state result emerges...


Determine maximum frequency of input signal to make system LTI This is question 4.26 from the third edition of Alan Oppenheim's textbook "Discrete-time Signal Processing". It is stated as follows:

Firstly, the answer is the input signal should be bandlimited to 12 kHz. I'm having trouble seeing why it is the case.

My reasoning:

The cutoff frequency of the low pass filter is obviously 8 kHz, so any frequency component of the discretized input signal beyond 8 kHz will be removed. But since we do not want aliasing to occur, we want the input signal to be bandlimited such that after passing the sampled signal through the low pass filter we don't get aliasing below 8 kHz. But in the case of input signal bandlimited to 12 kHz, we get the following (a quick drawing):

which causes aliasing, making the overall system no longer LTI.

So where did I go wrong? Any help is greatly appreciated!

This question, and a number of similar others, from DTSP book can be a little tricky to recognise the fact that it's not actually asking alias free operation (which would require 8 khz cutoff signal as you expected), rather it is actually asking how much aliasing (due to initial sampling block) is tolerable: you can allow aliasing in those regions of the DT filter $H(e^{j\omega})$ whose stop band will help you get rid of the aliased spectrum, hence enable an LTI processing possible, but not with the samples of the complete original fullband signal $x_c(t)$ but with that of the reduced bandwidth input signal.

To see this you shall better consult to the spectral figure below:
The mathematical condition on the maximum allowable cutoff frequency $\Omega_c$ of the input signal, whose aliased spectra will be deleted by the ideal lowpass filter with cutoff frequency $\pi/2$, is given as:

$$\Omega_s - \Omega_c > \frac{\pi/2}{T_s}$$
$$\frac{2\pi}{T_s} - \Omega_c > \frac{\pi/2}{T_s}$$
$$\Omega_c < \frac{2\pi}{T_s} - \frac{\pi/2}{T_s}$$
$$\Omega_c < \frac{3\pi/2}{T_s}$$

Thats stated in hertz frequency as:
$$f_c < \frac{3}{4}f_s$$

which gives a maximum allowable input signal bandwidth as $f_c=12,000$ Hz when the sampling rate is $f_s=16,000$ Hz.


How can I better extract the hidden images in FEZ's soundtrack? I believe I have found a (new?) secret in FEZ's soundtrack, but my visualizer isn't clear enough to capture it 100%.
In FEZ's soundtrack (by Disasterpeace), there are "secrets"- stuff like encoded images. These exist in almost all of the songs, even some of the more atmospheric ones. The game is about a shift in dimensional perspective, like from 2D to 3D and back, and the soundtrack represents that same shift beautifully when you pull up a spectrogram.
Every note is deliberate, in either creating an emotion or image, or even just setting up the rest of the song. White/blurry noise is used to create a sense of space, while mirrored sounds (high-low frequency, like a guitar playing the opposite notes, ascending instead of descending) give a sense of either spookiness or The Ancient Ones (the space squids that mastered the fourth dimension or whatever). Every. Single. Track... has either a secret, is atmospheric for the hell of it, every note deliberate. Disasterpeace did a good job on not making any "throwaway" tracks.
So while messing with setting up a spectrogram for streaming, I remembered the FEZ soundtrack existed, and had a listen to my favorites... and saw something regular. Very regular. HORIZONTALLY, not just the usual vertical of drums or shifting horizontal lines/planes of synths. Tweaking the settings, I got the best screencap I possibly could, but as you can plainly see without the context of the game, it's too blurry to properly see anything besides the telltale grid and tiny lines within, juuuust clear enough to recognize a few, but not enough to read at all.
The song is Sync by Disasterpeace.
This is the original capture:

This is the same image but cropped with the context on top:

Original Upload - Imgur
This is the FEZ alphabet (for reference), in case you don't know what you're (possibly) looking at:

A few pointers to help you produce more detailed results:
First, you may want to read about a phenomenon called Pareidolia. By tweaking and tuning parameters, you may be simply creating results through coincidence.
There is no such thing as the spectogram.
There is such a thing as a spectogram.
I've looked at some of the results already found that have been posted online, but none actually talk about the specifics of the spectogram used. This makes those results hard to reproduce (outside of running the exact same program with default settings). Why?
Spectograms have many, many parameters. Some explanation below:

Windowing function
First, there is the windowing function, which is technically any function taken from a function space (this stuff gets very complicated). Typically it's taken to be continuous or at least semi-continuous (a finite number of discontinuities).
Theoretically, this kind of space is of the order aleph-3, although practically we're limited by the resolution of the computer, or, the total resolution is N floating-point values big. Typically these have a 53-bit mantissa and are between zero and one, so we're already looking at 53KB of parameters for a 1024 FFT size spectogram.
Window size
This is typically taken to be a multiple of the FFT size, usually W = N. Ideally, at least one full FFT size is inside the window size. Complicated mathematical expressions can identify the resolution of a spectogram, which depends on this parameter.
Window overlap
In order to improve time resolution somewhat, when using non-rectangular windows windows are usually chosen to overlap to reduce spectral leakage losses. High window overlap factors tend to smear out/deform the image. Typically chosen to be 0.5 with most standard windows as that satisfies COLA.
FFT size
Changing the FFT resolution can change the shape and image of the spectogram through rounding effects and also has an effect on the time-resolution capability of the spectogram. You could compare it to the Heisenberg principle: increasing the sample count improves frequency resolution, especially at low frequencies, but reduces time resolution.
Must be a power of two for the traditional FFT algorithm to run quickly enough.
Has an impact on the speed of the visualization. Older programs tend to run with lower FFT sizes to enable real-time visualizations on slower computers. FFTs scale with O(N log N). See here.
Sampling frequency
May not match the recording frequency, in which case the audio is resampled before analyzing. Iff so, then the resampling method used can also impact the result.
Base frequency
Typically not directly configurable (typically the user is asked to select min/max frequency shown in the resulting image), the base frequency of the FFT is equal to fs/N, or the sampling frequency
Filters, transformations, etc.
Post- and pre-processing is typically used to improve the output quality for visualizers. Typically done to remove noise, compression artifacts and the like, they can corrupt hidden images, which are often of very low volume and may be caught by a noise filter.
Image width in px
Finally, the chosen image width and, iff different from the size of the FFT output, the resize filter can make an image fuzzier in quality. Ideally the image width is equal to the FFT output, though, if this is too big, the FFT output is an integer multiple of the image width. Take care to convert the image if needed with an integer ratio as well to reduce artifacting, and to use a lossless format.


Downconversion / Demodulation of RF signal to DC My question is around down conversion of a signal (e.g. RF) to DC and the affect of negative frequencies.

If I have a signal whos spectral content is say between 2-4kHz and I mix this with an intermediate frequency of 3kHz, ignoring the upper mirror, I now have a signal whos spectral content is centered around DC +- 1kHz with a bandwidth of 2kHz.

What does the signal at -1Khz to DC represent- I believe this is mirrored onto the DC-1kHz signal but does this create artifacts if the spectral content is asymmetric?

"If the signal spectrum is asymmetric, the original signal spectrum will be corrupted"

If your mixer is a simple multiplication with a single oscillator, then yes, that will be the case.

As Dilip explained, real (RF) signals are always symmetrical in spectrum to the 0Hz axis. That means your original doesn't only exist around 3kHz, it also exists around -3kHz, only mirrored!

I think your image captures that pretty fine:

Now, if you shift your 3kHz signal (black) to 0Hz, the higher spectral peak ends up on the positive frequency side, and the lower on the negative frequency side.

Problem is that for real-valued signals, there's no physical difference between positive and negative frequencies – and thus, your mirrored signal also gets overlayed with this. So you have the small and the big spectral peak overlay - and since they were originally two different parts of your signal, this overlay is not actually just a "shifted copy" of your original signal, but a "corrupted" version.

The only way this wouldn't happen is if the "upper" half of your signal (3-4 kHz) and the "lower" half (2-3 kHz) were already symmetrical in spectrum to its center frequency – because then, the your mirrored spectrum would be identical to the original, and nothing was lost.

But: simply read on. The fact that your text uses negative frequencies definitely works towards showing you how a complex baseband mixer works. In two or three pages, you'll understand what needs to be done to still be


Fourier transforms and time shift There is probably something trivial behind this, but I am missing something. I need to create a stationary random time series data v(t) which is the the sum of another time series u(t) and u(t) with a time shift $\tau$, i.e

$v(t) = u(t) + u(t - \tau)$

The series u(t) is itself a stationary colored Gaussian time-series with some PSD $S_u(f)$. No problem! I wrote a small function (using Python) to make u(t) for a given PSD. I calculate the PSD of the the generated time series using Welch's method and it checks out.

If we take the Fourier transform of u(t - \tau) it should just gain an extra phase compared to the Fourier transform of u(t)

$F[u(t - \tau)] = \int dt \, u(t -\tau) \, e^{2 \pi j f t } =  \int dt \, u(t) \, e^{2 \pi j f t }  e^{- 2 \pi j f \tau } = \tilde{u}(f) e^{- 2 \pi j f \tau }$

Since $v(t) = u(t) + u(t - \tau)$,

$\tilde{v}(f) = \tilde{u} (f) \, (1 + e^{- 2 \pi j f \tau })  $

And the PSD of V should then be (using Euler identity)

$S_v (f) = 4 S_u (f) \cos^2 (\pi f \tau )$

In particular $S_v$ should vanish when $f = \frac{2n + 1}{2 \tau}$ with the smallest zero being at $\frac{1}{2 \tau}$. However when I test my this with the u(t) series I generated above the PSD looks different. There are zeros low frequencies which shouldn't exist. For example  here is the PSD of the v(t) series and the spectra I would expect for $\tau = 2$. (The first zero should be at 0.25 Hz)

Moreover the zeros in the data seems to change with the sampling frequency I am using. The v(t) series for the above plot has a sampling frequency of 0.25 Hz. Here is one with 0.5 Hz, keeping everything else the same.

I am kind of baffled by this. Any help would be greatly appreciated. Sorry about the long post, and thanks in advance!!

EDIT: After using a higher sample rate as suggested by Hilmar below, and fixing the error I described in the comment everything looks good. Spectrum for the fixed case below.

Your math is correct, for a delay of 2 seconds, the first null in the PSD should be at 0.25 Hz.

I'm guessing the problem is in the code you use to generate the signals and/or the PSD. The sample rates seem awfully low. At a sample rate of 0.25 Hz, your 2 second delay is just half a sample, that means you would have to implement a fractional delay to generate the signal. T


MFCC process confusion I have read through many papers related to the Mel-frequency cepstrum coefficient (MFCC): I have some confusion about the MFCC process:

Many of the papers say that I have to do a FT (Fourier Transform) on the log-spectrum to obtain cepstrum (as mentioned in Wikipedia : definition of power cepstrum)
Many say I have to do an IFT (Inverse Fourier Transform) on log-spectrum to obtain cepstrum.
Many refer to the DCT-II (Discrete cosine Transform-2) on the Mel-log-spectrum to obtain MFCC.

My question is, which one is right:

FT or IFT?
If the IFT is correct than why are they using DCT-II instead of IDCT-II, since the DCT-II is similar to DFT or FT?

(hence, the IDFT or IFT should be similar to IDCT-II)

First of all, FT and IFT of a discrete sequence are the same up to complex conjugation/scaling - these two operations will thus yield qualitatively the same result (in the sense that subsequence processing/statistical modeling will make abstraction of these differences anyway). The textbook definition of cepstrum requires a Fourier transform. People sometimes use the inverse Fourier transform to emphasizes that the X axis on the cepstrum plot can indeed be interpreted as a time (rather than as a quefrency).

The MFCC extraction process found in speech and audio analysis systems uses the DCT since it has the advantage of being a real -> real transform, and because the resulting coefficients are in practice decorrelated from each other - which is not the case if you take the sequence of real/imaginary part of a FT or IFT. This makes the resulting data very good for subsequent modeling with gaussian models with diagonal covariance matrix.

Another view on this: take thousands of hours of audio, extract log mel frequency coefficients (through a bank of filters with triangular responses, with center frequencies distributed on the Mel scale) for each analysis frame. You'll get a big dataset... Do a PCA on it to extract the dimensionality reduction transform that captures most of the variance of the mel coefficients. The transform matrix the PCA will produce will be uncannily close to the DCT matrix.


Frequency response of marginally stable LTI systems The frequency response of a system is defined as:
$$\int_0^\infty{h(t)e^{-j\omega t}dt}$$
where $h(t)$ is the impulse response. But in marginally stable systems, $h(t)$ does not decay so the integral doesn't converge. Does the frequency response "exist" in this case?

As long as the system is linear and time-invariant, it can be described using the frequency response / transfer function $H(j\omega)$. For marginally stable systems, $H(j\omega)$ does exist, but cannot be expressed in terms of the Fourier-transformed impulse response $h(t)$.

The Form $H(j\omega) = \frac{Y(j\omega)}{X(j\omega)}$ however is still valid for systems not stable in the more narrow sense simply because they are both linear and time-invariant.


About length of the signal In my research work i have the equation y(n) = hx(n)+w(n). My understanding about this equation from MATLAB point of view are:

h has to be a single number. (for e.g., rayleigh generated).
x(n) , w(n) will be vector because they have n associated with them. So is my understanding correct or is their any other concept. Thanks.

I beleive that h can be a vector or scalar .. depends on channel taps i.e. h

y[n] = hx[n] + w[n]  -- indicates  h convolved with x[n]  and then added to noise w[n]


Can convolution of one signal with different signals give the same answer? Let us consider $x_1(t)$, $x_2(t)$, $x_3(t)$, all the same within some some duration 0 to $T$ but all different outside this interval.  Now let us multiply each of these signals with $w(t)$, a window function - nonzero from 0 to $T$ but 0 outside this interval.  So the multiplication of this $w(t)$ with each of $x_i(t)$ will give the same signal.  This should be same as Inverse Fourier Transform of $W(\omega)$ convolved with $X_i(\omega)$.  Does this mean convolution of the same signal with different signals can give the same result?  Any comments?  Is there a pitfall in my interpretation?

Also I see mathematically equations neatly entered in questions and answers on this site?  Where should I start to learn on how to write the equations?

Convolution is filtering.

Consider creating a filter with a rectangular response in the frequency domain.  Then any signal with the same content within the filter's passband, no matter what the signal content outside the passband of the filter, should result in only the content within the passband after convolution with the impulse response of the filter.  Thus, multiple inputs (signals with different stuff outside of the passband) to a convolution with this filter should result in the same output.

(Your opening statement merely swaps the frequency and time domains for the rectangular function, making this relationship harder to recognize.)

This should work even for a non-rectangular and finite length FIR filter with at least one zero in it's transform (maybe at a nice integer submultiple of the sample rate).  Add any magnitude of sinusoidal input at the frequency of the zero, and the convolution result or filter output remains unchanged (except perhaps for numerical issues/clipping/quantization/etc.).


Does MATLAB code retain its accuracy when is executed in dsp kit? Ladies, Gantlemen,
Please let me ask your help on my problem. MATLAB code (at https://metrw-pitch.blogspot.com/ you can see code in .mlx (for run it you want delete comments) and .m format) executing in MATLAB online R2020b, outputs good accuracy even for low fundamental frequencies (fundFreq) but code (originally in my answer to a question of mine) of similar algorithm in C++ (acctually I have translated this code into Matlab's language), executing in Visual Studio 2019, outputs bad accuracy in low (lower than 150 Hz) frequencies.
Question is whether MATLAB's code executing in dsp kit retains its accuracy. I'm homeless (in France) and get internet access in Social Services and Post Offices like this I'm just now. So it's impossible for me to run code in dsp kit. You are kindly requested test MATLAB's code in your dsp kits and say me its output for fundFreq 50, 60, 70 Hz. Be aware when change FundFreq, change, if necessary, FFTfundFreq and/or GridDemiSpan so that grid's span covers FundFreq.
I have asked similar question on 8 Feb 2021 at 13:25, at MATLAB but nobody answered yet.
Copy of Matlab's code without comments is following.
format long;
SampFreq = 16000;
Segm = 1:1600;
FundFreq = 50;
FFTfundFreq = 41;
GridDemiSpan = 10;

FirstHarmAngles = FundFreq*2*pi/SampFreq*Segm+1.9*pi;
SinFirstHarmAngles = sin(FirstHarmAngles);
SecondHarmAngles = FundFreq*2*2*pi/SampFreq*Segm+0.9*pi;
SinSecondHarmAngles = sin(SecondHarmAngles);
ThirdHarmAngles = FundFreq*3*2*pi/SampFreq*Segm+0.3*pi;
SinThirdHarmAngles = sin(ThirdHarmAngles);
Xn = 170000*SinFirstHarmAngles+220000*SinSecondHarmAngles+...
150000*SinThirdHarmAngles;

Freqs = FFTfundFreq-GridDemiSpan:0.1:FFTfundFreq+GridDemiSpan;
MagnSqrd = ones(1,201);

for f = 1:201
   Angles = Freqs(f)*2*pi/SampFreq*Segm;
   XnCos = sum(Xn.*cos(Angles));
   XnSin = sum(Xn.*-sin(Angles));
   MagnSqrd(f) = XnCos.^2+XnSin.^2;
end

[maxMagnSqrd, maxMagnSqrdIndex] = max(MagnSqrd);
GRIDfundFreq = Freqs(maxMagnSqrdIndex);
disp(GRIDfundFreq);

With regards and friendship, Georges Theodosiou

Yes, C\C++ code generated from Matlab/Simulink code should behave the same as the original code.
However there are a few caveats, here are 3 problems I've experienced in real-life.
1 - Simulink\Matlab performs an implicit reset just before starting the simulation/script. Make sure you perform this reset in your C\C++ code, otherwise your results won't be the same and won't be comparable. Simulink coder generates a reset function for you. I've seen a corner case where the reset function would not reset all global variables, so be careful.
2 - Variable types. Matlab/Simulink typically use double-precision floating point variables. If your code uses doubles then great. The behaviour should be the same. However, for performance reasons, you might choose to use single-precision floating point variable aka floats. Then the behaviour could change slightly. There are ways to mitigates the effect of quantization, some algorithms are really sensitive to precision, others less so. Even if you use doubles there are some difference between platforms. So if you use doubles on an ARM µProcessor, the behaviour could be slightly different in corner cases.
3 - Integer saturation. I made the mistake in the past of using a 32-bit integer in Matlab as a free-running counter. It worked well, but after like 10 days the code would freeze. The counter would be stuck at 2^32-1 instead of wrapping back to 0. Obviously, you won't see that behaviour in simulation. So be careful, test your code thoroughly.


Autocorrelation and Power Spectral Density (Discrete) The Autocorrelation, $\phi_{aa}[\kappa]$, of a discrete time random process, $a[k]$, is defined as:

$$ \phi_{aa}[\kappa] = \mathrm{E}\left\{ a[k+\kappa]a^*[k] \right\} $$

Taking its fourier transform would give the Power Spectral Density, which is given in my book as:

$$ \mathscr{F}\big\{ \phi_{aa}[\kappa] \big\} = \Phi_{aa}\left(e^{j2\pi fT} \right) = \sum\limits_{\kappa = -\infty}^{\infty} \phi_{aa}[\kappa] e^{-j2\pi\kappa fT}  $$

Now, why don't they denote simply as: $\Phi_{aa}[\kappa]$?

Is the form above with the exponential signify something related to sampling interval? I'm quite confused with this form and any explanation to better understand this form would be helpful.

First of all lets state more correctly that a discrete-time 1D auto-correlation sequence (ACS), $\phi_{XX}[\kappa]$, of a single parameter $\kappa$, of a discrete-time random process $X[n,s]$ will exist as long as the process is stationary or at least WSS (wide sense stationary).

And therefore the power spectral density (PSD), $\Phi_{XX}(e^{j2\pi fT})$, of the process will also exist iff the random process is stationary (or WSS at least) and its ACS is stable, in the sense that its Fourier transform will exist. (In case ACS is not stable, generalized definition of the Fourier transform is used to represent nonconvergent integrals with impulses or with discontinuties at those particular frequencies...)

Note that the notation changes slightly and I prefer to use $\phi_{XX}[m]$ to denote the ACS of the WSS random process $X[n,s]$ and use $\Phi_{XX}(e^{j\omega})$ to indicate its PSD, which is indeed given by the discrete-time Fourier transform (DTFT) of its ACS :
$$ \Phi_{XX} (e^{j\omega}) = \mathcal{F} \{ \phi_{XX}[m] \} = \sum\limits_{m=-\infty}^{\infty} \phi_{XX}[m] e^{-j\omega m} $$

Since the PSD of a WSS random process $X[n,s]$ is a function returned by the DTFT of a sequence, it must be a periodic and continuous function of frequency $\omega$ , as all other Fourier transforms. We are only interested in the first period from $\omega=-\pi$ to $\omega=\pi$ (or from $\omega=0$ to $\omega=2\pi$)

Coming to your confusion, in most practical applications, a discrete-time random process $X[n,s]$ is represented by a discrete-time random sequence $x[n]$ which is created by the process of sampling a particular realization, $x(t)$, of a continuous-time random process $X(t,s)$ with a sampling period $T_s$. In such a case $\omega = 2\pi f T_s$ could be used as indication of the sampling relationship between the continuous-time and the discrete-time random processes involved. Also in communications terminology, frequency in Hertz $f$ is preferred over frequency in radians per second $\Omega$.

Note also that I've deliberately used a capital letter $X$ , instead of a lowercase $x$, to denote the ownership of those ACS and PSD functions; They belong to the random process $X[n,s]$ and not to a single particular realisation sequence denoted by $x[n]$. The sequence $x[n]$ belongs to the family of all functions that constitude the ensamble of the random process $X[n,s]$ ($n$ denotes the (time) index and $s$ denotes the outcome of the experiement for those indexed random variables $X_n(s)$ associated with the process $X$)

The random process $X[n,s]$ exists as a theoretical entity, an abstract being with its own characterisation that generates all those possible realizations $x[n]$ being observed. In practice we can only observe one or more $x[n]$ and therefore we have to estimate the true (but unobservable) quantities of ACS and PSD that belong to random process $X[n,s]$ from those realisations by means of statistical estimation methods. This process is not exact, hence in practice we can only compute an approximation to those true theoretical entities of ACS and PSD from particular obsrevations of instances $x[n]$ of the random process $X[n,s]$.


Energy Spectrum of a signal after FFT in Matlab I have a vector of signal x(t) with its time vector. I want to obtain a frequency representation of the signal, in particular the energy spectrum of x(t).

Can someone please show some light into obtaining the energy spectrum ?
Is energy spectrum and energy spectral density ,ESD produce the same result ?

Any help is much appreciated !

 t = 0:.001:.25;
x = 5*sin(2*pi*50*t) + 1.2*sin(2*pi*120*t);
y = x + 2*randn(size(t));
y=y';
Fs = 1000; % Sampling frequency

% Fast Fourier Transform
    L = length (y); % Window Length of FFT
    nfft = 2^nextpow2(L); % Transform length , Next power of 2 from length of signal
                          %  Exactly same as nfft = pow2(nextpow2(L))

% Take fft, padding with zeros so that length(input) is equal to nfft
    Ydft = fft(y,nfft)/L;

    y_HamWnd = y.*hamming(L);           %hamming window fft
    Ydft_HamWnd = fft(y_HamWnd,nfft)/L;

%FFT is symmetric, throw away second half

    % Take the magnitude of fft of x
%Amplitude Spectrum 01 (Single Sided)
mY = 2*abs(Ydft(1:nfft/2+1)) ;

%Amplitude Spectrum 02 (Single Sided HamWmd)
mY_HamWnd = 2*abs(Ydft_HamWnd(1:nfft/2+1));

%Amplitude Spectrum 03 (Single Sided; Exclude DC-- FoldingFreq)
mY2 = abs(Ydft);
mY2= mY2(1:nfft/2+1);
mY2 (2:end-1) = 2* mY2 (2:end-1);

    %To get the Power
% P ; single sided; then coefficient square**********Multiply by 2 ?
YdftCoeffi_Sqr = Ydft(1:nfft/2+1).^2;

% P = y*(conj y) == y^2  ;     %********DO WE NEED TO MUL 2 ?
Pyy_NormNfft = Ydft.*conj(Ydft);

% P ; single sided; First square, then Multiply by 2
m1YL = abs(Ydft(1:nfft/2+1));
mSqrYL = m1YL.^2;
m2YL = mSqrYL*2;

% Frequency vector
f1 = Fs /2 * linspace(0,1, nfft/2 + 1);
f2 = (0:nfft-1)*(Fs/nfft);

% Generate the plot, title and labels.

figure;
plot(t,y); title('Time Domain'); xlabel('Time(s)');

figure;
subplot(3,1,1)
plot(f1,mY)
title('Amp Spec 01 --mY = 2*abs(Ydft(1:nfft/2+1)) /L; '); xlabel('Frequency (Hz)'); ylabel('|mY|')
        %axis ([0 500 0 5]); %Zoom in
subplot(3,1,2)
plot(f1,mY_HamWnd)
title('Amp Spec 02 -- mYHamWnd =2*abs(YdftHamWnd(1:nfft/2+1))/L;'); xlabel('Frequency (Hz)'); ylabel('|mY_HamWnd|')
        %axis ([0 500 0 5]); %Zoom in
subplot(3,1,3)
plot(f1,mY2)
title('Amp Spec 03 DC-- FoldingFreq -- mY2 ;First, single sided; then square, lastly Multiply by 2'); xlabel('Frequency (Hz)'); ylabel('|mY2|')
        %axis ([0 500 0 5]); %Zoom in

figure;
subplot(3,1,1)
plot(f1,YdftCoeffi_Sqr);
title('Ydft Complex Coefficient Square; Single Sided-- YdftCoeffi_Sqr = Ydft(1:nfft/2+1).^2');  xlabel('Frequency (Hz)'); ylabel('YdftCoeffi_Sqr');
subplot(3,1,2)
plot(f1,m2YL);
title('Ydft Complex Coefficient Magnitude Square; Single sided; then Multiply by 2');
xlabel('Frequency (Hz)');
ylabel('Power m2YL');
        %axis ([0 500 0 5]); %Zoom in

The cited passage gives a method for obtaining the energy in the harmonics of a specific frequency.  If your frequency of interest is $f_0$, then take a section of samples of length $N$ where $N = \frac{4 F_s}{f_0}$ and compute the $N$-point FFT.  Now you can find the energy of the second harmonic by squaring the magnitude of the $4k+1 = 9$th FFT bin. To account for errors in estimating $f0$ you can sum the squared magnitudes of bins 7 through 11.

Below is some example code that might give an idea of what the algorithm in question does:

%% Generate an example signal
Fs = 5e3;                     % Sampling rate
f0 = 151;                     % Frequency of interest
t = (0:Fs-1)/Fs;              % Time samples
x = cos(2*pi*f0*t);           % Pure sinusoid @ f0
noise = 0.1*randn(size(x));   % Gaussian noise
y = (x>0)-(x<0) + noise;      % Introduce harmonics and noise

%% Compute energy in the harmonics of f0
N = round(4*Fs/f0);           % Determine record length
Y = fft(y(1:N))/N;            % Compute N-FFT (normalized)
Y2 = Y.*conj(Y);              % Magnitude squared
K = floor(N/2/4);             % Number of harmonic "bins"
for k = 1:K
    idx =  4*k+1+(-2:2);      % Five indices centered at 4k+1
    Pyk(k) = sum(Y2(idx));    % Sum of energies around 4k+1
end

%% Plot the results
[Pyy,Fyy] = pwelch(y,[],[],[],Fs);
plot(Fyy,10*log10(Pyy));
hold on;
F = linspace(0,Fs-Fs/N,N);
Fk = F(1+4*(1:K));
plot(Fk,10*log10(Pyk),'ro');
title('Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/rad/sample)')
data2 = sprintf('Energy in Harmonics of %d Hz',f0);
legend('Power Spectral Density Estimate', data2);
hold off;

As far as finding a "signature" in the data, you can see that the algorithm could be very useful for that.  However, it is very important to note that the algorithm presumes some fundamental frequency of the signal in question.  If the signal in question has a fixed period, but the shape of the signal over that period changes with time, then this algorithm would be useful.  However, if the frequency varies with time, then this algorithm is not a good fit for your application.

If you are simply looking to estimate the spectrum of a signal, then maybe this page will be of help: FFT for Spectral Analysis (mathworks.com)


Intuitive explanation of tracking with Kalman filters I would much appreciate an intuitive explanation for (visual) tracking with Kalman filters.
what I know:

Prediction step:

Dynamic system state $\mathbf x_t$: target location at time $t$
Measurement $\mathbf z_t$: the image at time index $t$ (??)

Based on images/measurements $1\rightarrow(t-1)$ I want to predict state $\mathbf x_t$? (using dynamic equation) Is that correct ?

How can I interpret the correction step into those terms (image, target location) ?

First you have to assume a motion model.  Let's say you wish to track a ball flying through the air.  The ball has a downward acceleration due to gravity of 9.8m/s^2.  So in this case the constant acceleration motion model is appropriate.

Under this model, your state is position, velocity, and acceleration.  Given the previous state you can easily predict the next state.

You also have a notion of detection. You have a video of the ball moving, and you somehow have to detect the ball in each video frame (e. g. by using background subtraction).

Your detections are noisy.  Also, the motion of the ball does not exactly fit the constant acceleration model because of air resistance, wind, cosmic rays, etc.  The Kalman filter needs two matrices describing this: one is the covariance of the measurement noise (your detection inaccuracy), and one for the covariance of the process noise (how the motion of the ball deviates from the model you have specified).

If you are tracking a single object, then the Kalman filter lets you smooth out some of the noise, and also predict where the object is when a detection is missing (e. g. if the object is occluded). Here is an example of tracking a single object with a Kalman filter using the Computer Vision System Toolbox for MATLAB.

If you are tracking multiple objects, then the Kalman filter predictions let you decide which detection goes with which object.  A good way to do that is to use the log likelihood of the detection given the error covariance of the prediction. Here is an example of tracking multiple objects with a Kalman filter.


Remove noise introduced while converting a Delta-Sigma Modulated (DSM) signal into a pulse-code modulated (PCM) signal? I'm writing a tool that converts a dsdiff audio file (dsd64) into a flac audio file (DSM -> PCM). The project is both a programming exercise and a signal processing exercise; I'm learning a lot from it.
Currently, I'm doing the following to compute a sample in the PCM stream

Select a window of the DSM sequence that corresponds to the size of the window of a low-pass (FIR) filter.
Replace all 0s in the DSM window with -1s
Convolve (multiply+sum) the elements in the DSM window with the coefficients of the (FIR) filter and rescale the result to fir into the chosen bitrate (24-bit).
Move the window by downsampling factor many elements (16 in my case) and apply repeat this process.

I tested the conversion using a few (arbitrarily chosen) dsdiff files (songs by the band deep purple).
The conversion reproduces the signal (good) but introduces noise (bad) with the following profile (measured during the beginning of the track where there is a bit of silence in the original):

I've been trying to figure out how this is being introduced, and what I need to do to remove it. From my reading of random sources online I understand that this is introduced by integrating the DSM signal, which amplifies the noise floor and maps it to a spectrum similar to the one I see in my plot.
Did I identify the problem correctly? What can I do to verify this and, ultimately, fix it?
Any help is highly appreciated.
Edit:
Added a figure showing the filter's frequency response

Added the Periodogram of the first second of the DSM stream

The spectrum appears to be the noise-shaping associated with Delta-Sigma Modulation. Specifically we see that the noise increase follows a 20 dB/decade slope, indicative of a first-order delta-sigma modulator.
The high end of the noise rolls off consistent with the high frequency roll-off of the filter applied.
Below are some charts I have explaining the first-order Delta-Sigma modulator and how such noise shaping occurs, here specifically with a Delta-Sigma Digital to Analog Converter (but it could easily be a 1 bit digital output just the same).

And equivalently the Delta Sigma Analog to Digital converter diagram is shown below, which will be used to develop a simple loop model for providing an understanding of how noise shaping occurs.

Below is a qualitative loop model, providing an intuitive understanding of noise shaping if one has an appreciation of filtering in simple control loop models. Specifically it is showing an 1 bit A/D Converter as a gain with quantization noise added (this is a highly non-linear process, therefore the model is qualitative versus a quantitative prediction).

The path through the loop from the quantization noise input ($N_q$) to the output follows a high pass function as shown below. Quantization noise would typically be white noise (uniform over all frequency) for a traditional data converter, while with the implementations using Delta Sigma modulation we get this resulting noise shaping resulting in significantly higher precision at a given sampling rate after the high frequency noise is filtered out.

Typically higher order implementations are used, such as the 2nd order Delta Sigma shown below (and higher orders such as 5th order are commonly implemented). This provides much steeper noise shaping (much less noise in the band of interest after filtering) and also serves to significantly reduce pattern noise which would manifest as spurs in the spectrum.


Finding the output of frequency response Suppose the input to an FIR system is

$x[n]=e^{j(2\pi n/7 -\pi/2)}$

Define a new signal $y[n] = x[n]-x[n-1]$ . The signal $y[n]$ can be expressed in the form

$y[n] = Ae^{j( 2 \pi f n + \phi)}$

Find $A$, $\phi$, and $f$.

I really would just like confirmation that I am correct.
I calculated $A = 0.8678$

$\phi = -1.7952$

$f = 0.2244$

Unfortunately, apart from the amplitude $A$, I can't confirm that you're correct. The first thing that you can see without doing any calculations is that $f$ must be the same for input and output, because the system is linear and time-invariant, so the frequency of the complex exponential remains unchanged. Consequently, $f=1/7$. The output signal can be written as

$$\begin{align}y[n]&=e^{j(2\pi n/7-\pi/2)}\cdot\left(1-e^{-j2\pi/7}\right)\\&=e^{j(2\pi n/7-\pi/2)}\cdot e^{-j\pi /7}\left(e^{j\pi /7}-e^{-j\pi /7}\right)\\&=e^{j(2\pi n/7-\pi/2)}\cdot e^{-j\pi /7}\cdot 2j\sin(\pi/7)\\
&=2\sin(\pi/7)\cdot e^{j(2\pi n/7-\pi/7)}
\end{align}$$

So you have

$$\begin{align}A&=2\sin(\pi/7)=0.8678\\
f&=1/7=0.1429\\
\phi&=-\pi/7=-0.4488
\end{align}$$


How to build a simple model for voice-controlled system? I want to build a speaker independent voice-controlled system. If I ask 10 people to say one unique word(for example"hello"), and record their voice and do LPC(linear predicting coding) upon them, I have 10 arrays of LPC coefficients, right?
So how should I combine these 10 arrays to build one array(my model of "hello word" and use it to compare with my system's new input voices? (I like to use DTW algorithm for this comparation).

Thanks!

DTW does not really work in speaker-independent systems, you'd better use a proper speech recognition toolkit.

For 10 recordings of a single speaker, you compare to all one by one and select the best match.


Strong (0,0) pixel among an otherwise normal cross correlation peak I am working on an imagery registration problem via cross-correlation and am seeing a weird phenomena sometimes. First, I am pre-processing each image by subtracting the mean image brightness from every pixel and then applying a cosine window to deal with the fact that my images are not periodic across the edges.

Then I am doing a simple cross correlation in Fourier space and creating the correlation surface. I work in Python so here is some scratch code showing what I am doing:

import numpy as np
from scipy.signal.windows import cosine as wind

im1 = np.load(file1)
im2 = np.load(file2)

rows, cols = im1.shape
hr = np.reshape(wind(rows, sym = False), (rows, 1))
hc = np.reshape(wind(cols, sym = False), (cols, 1))

h2d = np.dot(hr, hc.T)

im1 = im1 - im1.mean()
im2 = im2 - im2.mean()

im1 = im1 * h2d
im2 = im2 * h2d

fft1 = np.fft.fft2(im1)
fft2 = np.fft.fft2(im2)

prod = fft1 * fft2.conj()
surface = np.abs(np.fft.ifft2(prod))

When the relative shifts between the two images are largish, the (0,0) pixel seems to be drowned out by the overall correlation surface peak structure. However, when the relative shift becomes smaller, the (0,0) pixel starts to dominate and is no longer a smooth part of the overall surface. This causes a simple peak finding algorithm like numpys argmax to latch onto a shift of (0,0) even though a faint correlation peak surface still exists.

I am curious what can cause something like this to happen? In general, what about a set of two shifted images can cause the (0,0) pixel to dominate?

Edit 4: I discovered a scipy function called match_template that incorporates a normalization to the correlation surface proposed by JP Lewis. This method requires that template dimensions be smaller than the reference image. I accomplished this by zero padding the reference image to 2x the original image size so that upon convolution of the template, the entirety of the image dimensions contain valid correlation results. I found that there is basically no difference between my code above and this match_template method and the (0,0) peak is still the strongest component when it shouldn't be in the some cases. I think a band-aide solution is to interpolate the (0,0) pixel so that it becomes part of the overall surface. I'm jotting all this down hoping this helps someone in the future since this question hasn't gained much traction.

Edit 3: After some intensive literature review, it has become clear that this method of cross correlation is at best an approximation of normalized cross correlation. It seems like the denominator in NCC has no easily computable equivalent in Fourier space. Additionally, doing cross-correlation for images of the same size, inevitably results in inaccuracies anywhere besides the (0,0) shift component. This is due to pixels in the two images that no longer overlap once they have been shifted relative to each other. Spatially, this can/has to be accounted for by padding the images so that the correlation is defined everywhere. This padding results in inaccuracies but is necessary to search all possible translations between the two. I have another question related to this point. When doing a cross correlation like in the code above, what "padding" is applied? How are the non-overlapping parts spatially handled in Fourier space? I never define a padding, yet the surface I get is the same size as the input images...so what is going on behind the scenes for every pixel except the (0,0) one?

Edit 2: I accidentally stumbled on the realization that the (0,0) pixel of the correlation surface is always equal to the mean of the fft product. That makes me feel like its sort of deterministic and that the (0,0) pixel in the surface doesn't really encode any information about what the actual shift?

Edit to add an example image, notice the underlying structure isn't centered on that central (0,0) pixel yet it is the largest component of the surface. Also this surface has been shifted with np.fft.fftshift() to make visualization easier.

Figures that the same day I post a bounty I answer my own question.

The answer to this question is dead/speckle like pixels. For a fixed number of dead pixels of a given brightness, the smaller the relative shift between the two images, the fewer dead pixels that are needed to cause the (0,0) pixel to spike. It can take fewer brighter dead pixels to have the same effect as more dimmer ones. Lastly, the threshold between a strong (0,0) pixel and normal one is extremely tight in the number of dead pixels. Literally a change in just 10 dead pixels can massively increase the probability that the (0,0) pixel becomes a problem.

This problem can be solved by a combination of dead pixel finding/interpolation and applying a light Gaussian blur to the images to minimize speckle.


Why we take different mfcc components in different programs? Why we take different MFCC components in different programs? Suppose, in a program we take No. of mfcc components = 14.What does it signify?

The number of coefficients is just one of the many parameters to tweak when designing a classifier. For speech analysis it turned out that usually 8-14 coefficients are sufficient, and increasing the number of coefficients beyond 14 does not improve the performance for most applications. The exact number you choose will depend on your task at hand. Most systems I've worked with used 12 coefficients, and that was for speech recognition. If you're not sure, you can simply run some experiments and change the number of MFCCs in a reasonable range to see the effect it has for your application. But I have no doubt that you'll end up in the range that I've mentioned above.