A Working Signal Frequency Example? I've been attempting to work on signal transformation and filtering for the past few months. It's a bit difficult, but I've been making progress.

I was wondering if anyone had any simple examples to show loading in a wave file, editing the volume and frequencies of that file, and re-outputting it to a file? I've already gotten up to volume, but I just can't edit the frequencies.

I thought that by using FFT, I would gain direct access to the frequencies of the sound file. Was that incorrect? I recall looking at frequency bins and such, but I don't see how that ties into wave files.

I know it might sound really complex, and I don't want to ask for a lot,  but could someone post or link to an example? It would really help me out.

I'm using Python, with Windows 64-bit. I'm already using Scipy and NumPy for FFT and wave file input / output.

I want to basically be able to low-pass or high-pass filter a sound by code, as well as be able to read the frequency data in each sample (i.e. see how high- or low-pitched the sample is).

As you have observed, instantaneous signal amplitude is an easy to understand, one-dimensional quantity, whose averaged absolute value is loosely associated with perceived loudness, and which can be easily altered by multiplying your signal by a constant.

Frequency is a very different affair. If we lived in a world where the audio signals were made of sustained pure tones, it would make sense to ask "what is the frequency of this sound?". But things are more complicated:

Stationary speech or music signals can be described as a sum of pure tones. Thus, it makes more sense to ask "what is the signal amplitude at this particular frequency?". This is why we look at spectra - graphs showing signal amplitude as a function of frequency. Note that frequency is not the plotted  quantity, it is the $x$ axis!
Speech or music signals are not stationary: some components of the signal decay over time or are modulated, some appears... Thus, it makes even more sense to ask "what is the signal amplitude at this particular frequency and particular time frame?". This is why we commonly see for audio analysis a representation called "spectrogram" - showing the measured signal energy as a function of time and frequency (again, note that frequency is the axis, not the measured quantity).
In one of your comments you say "as well as be able to read the frequency data in each sample". Frequency is a rate of change, and to observe a change, we need a collection of samples - it makes no sense at all to do any frequency analysis on a single sample, and if you use Fourier analysis (which is just one of the ways of estimating a sum of sine waves from a signal), you'll need more and more samples to get a more and more precise resolution on the frequency axis. But remember that sounds change over time, and that performing a Fourier Transform discards timing information (the magnitude spectrum of two simultaneous 100 Hz and 300 Hz tones is the same as the spectrum of a 100 Hz tone followed in time by a 300 Hz tone). Thus, if you care about the temporal position of events, you will need to perform your analysis on slices of the signal, computing a FFT on overlapping blocks of the input signal - this is called the short-term Fourier transform (STFT). You'll have to deal with resolution problems here - the smaller your windows, the more accuracy you have on the time axis, the less accuracy you have on the frequency axis. The bigger the analysis windows, the less accuracy you have on the time axis, the more accuracy you have on the frequency axis. Again, and again, note that frequency is not what you estimate - your STFT is just answering the question "What is the signal energy in this particular frequency range and this particular time frame?".

I hope these will make you realize that there's no such thing as "the frequency" of an audio signal. Frequency is an axis, "where we look at things", not a measured quantity, "what we measure".

Note that there are some one-dimensional quantities homogeneous to frequencies (and thus expressed in Hz), that can be used to characterize stationary (stable) sounds. For generic audio or music, these quantities will vary over time:

Fundamental frequency, is the inverse of the period of the sound. It relates to the perceived musical pitch (note) of the sound. This quantity is undefined for aperiodic sounds. Pitch-shifting algorithms will alter this quantity, to make a recording sound as if it was played/sung higher or lower on the musical scale.
Roll-off frequency, is the frequency under which a given fraction (say 95%) of the energy of a sound is contained. This quantity can be affected by low-pass / high-pass filtering, to make a recording sound "brighter" or "duller". It is one of the many dimensions of the perceived timbre of sound.

Maybe some of these are the quantities you really care about?

One final word about the FFT - it is the hammer in a signal processing version of "if all you have is a hammer, everything looks like a nail". It is an extremely useful tool, but it is not the best fit for many audio signal processing/analysis problems, and more often than not it is used only as an intermediary processing or computational step. In particular, for some of the questions you might be interested in:

Estimating the fundamental frequency might an involve a FFT, but the spectrum is not the right thing to look at, because of the "missing f0" phenomenon (a signal might be perceived as being a 263Hz C3 note, without having energy at this frequency). Don't make the mistake of trying to estimate the musical pitch of a signal by looking at the peak in its spectrum.
Computing a FFT, messing with the numbers, and performing the inverse FFT is a poor way of filtering a signal. This should be done with digital filters, available in scipy.signal.lfilter. You'll find formulae to compute IIR filter coefficients for a variety of musically useful second-order filters here.


Calculate a DFT of a simple finite signal {+1, −1, +1, . . .}) Hey i was wondering whether approach presented in this video (6:10) is a good approach to solve this task. Cause it seems kinda tricky and i was wondering whether it can be done faster.

Computing the DFT sum

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N},\qquad k=0,1,\ldots,N-1\tag{1}$$

directly can be tedious. A computer would generally use a Fast Fourier Transform (FFT) algorithm.

For the special class of exponential signals

$$x[n]=a^n,\qquad n=0,1,\ldots,N-1\tag{2}$$

$-$ and your example is a special case of Eq. $(2)$ $-$ you can find a closed-form expression for $X[k]$ by using the geometric sum formula. (Just be careful with those values of $k$ for which the denominator of that formula would become zero.)


Where did the k of $a_k$ disappear from Fourier Reverse Transform if $\omega=\omega_0k$? Where did the $k$ of $a_k$ disappear from Fourier Reverse Transform if $\omega=\omega_0k$?

We turn $\omega0$ to be $d\omega$, but $\omega=\omega_0k$, so shouldn’t there be a $k$ in the reverse transform?

Attached reference below:

It is because $\omega_o$ approaches 0 such that $k\omega_o$ becomes a continuous function.

So the discrete function $X(jk\omega_o)$ becomes the continuous function $X(j\omega)$.

This is intuitive if you consider the Fourier Series expansion in that each component of the spectrum is spaced in frequency by $1/T$ Hz, (or $\omega_o = 2\pi/T$ radians/sec). Each of these components are given as $k\omega_o$ with $k$ being the set of all integers. So as $T$ gets larger, $1/T$ gets smaller, and hence $\omega_o$ gets smaller and as $T$ approaches infinity, $k \omega_o$ approaches a continuous function.

The concept behind the Fourier Series Expansion is that any continuous time function defined over time duration $T$ can be decomposed into a sum of sinusoids that each have a fundamental period of $T$ and all higher harmonics such that there are an integer number of cycles in $T$, the the waveform can repeat for all time and start and end the same way at all boundaries $t= NT$ for all integers $N$. For this reason no other frequency components can exist and the non-zero frequency spectrum is therefore discrete.

The Fourier Series Expansion given above in terms of sines and cosines is much more compact using exponential frequency notation (and ultimately much more intuitive once you understand that $e^{j\Omega t}$ is a phasor rotating counter-clockwise on the complex plane at rate $\Omega$ radians per second) :

$$x(t) = \sum_{n=-\infty}^{\infty}c_n e^{jn\Omega_o t}$$
Where

$e^{jn\Omega_o t} = \cos(n\Omega_o t) + j\sin(n\Omega_o t)$

$c_n = a_n - jb_n$

And $\Omega_o$ was used in place of $\omega_o$ just to match the plots I already had (where I was distinguishing radian frequency between analog and digital domains).

This becomes the Fourier Transform when we all $T$ to extend to infinity: the spacing given by $1/T$ between each of the components goes to zero, and the frequency spectrum therefore becomes continuous.

An easy way to observe this effect is by simply zero-padding a time domain waveform prior to taking the FFT: as you add more zeros, meaning $T$ gets larger, more and more samples will be interpolated in frequency as the spectrum approaches a continuous function.


Relationship between wavelet shape and wavelet filters MATLAB has a library of wavelet functions, showing their "continuous forms" as well as the the decomposition and reconstruction filters.

In decimated wavelet transform the filter size remains the same and the data points are downsampled by $2^j$ at every level $j$. A question that I cannot find an answer is how are the wavelet filters and the actual continuous forms are related? I am chemist, exploring wavelets for some applications.

For example, MATLAB has db9. Its continuous estimation, wavelet function, psi or phi do not look similar to the decomposition filter points. One text says that some wavelets start out as filter and later their continuous forms are estimated, but why in the case of db9 and many other, the filter points do not match the shape of the wavelet function?

Thanks.

I will start with an averaging filter analogy. The simplest (non trivial) and cheapest way to average a discrete signal $x[n]$ is to  compute the mean of two adjacent samples, and remplace the two samples by this average.  This amounts to compute the discrete system:

$$\frac{x[2n]+x[2n+1]}{2}\,.$$

This can be recast as the convolution $$y[n] = x[n]\ast \frac{ [1\; 1]}{2}$$

You might remember that $\frac{ [1\; 1]}{2}\ast \frac{ [1\; 1]}{2} =  \frac{ [1\;2\; 1]}{4}$. And that  $\frac{ [1\;2\; 1]}{4}\ast \frac{ [1\; 1]}{2} =  \frac{ [1\;3\;3\; 1]}{8}$

This should remind you of the Pascal or binomial triangle, that evolve to a Gaussian shape:

Here, starting with a crude, and a simple finite-impulse response two-point averaging filters,  via iterations, you can end up with an IIR, not recursive, infinite-support Gaussian filter. This example is more a metaphor of how discrete filters generate a continuous function.

To be  able to add downsampling within the filters, without losing information, contraints the data smoothing picture. In the discrete context, this requires to have a downsampled smoothed version, and a difference  version that records when the data is not smooth enough to look like  its smooth version. This is the root to most practical discrete wavelet transforms. The Multi Resolution Analysis relates an (yet unknown)  scaling function $\phi$ at two-scales, via a linear combination, with weights corresponding to the filter $H$ coefficients. This is explained for instance in Signal Processing: Fourier and Wavelet Representations, 6.1. Introduction, Scaling Function and its Properties, page 191. Finally, the spectrum of the scaling function is given by:
$$\Phi^{(\infty)}(\omega)= \prod_{k=1}^\infty \frac {1} {\sqrt 2} H\left( \frac {\omega} {2^k}\right) \Phi^{(\infty)}(0)$$
where iterated half-cut-off frequencies are apparent. The take-away message is: the spectrum of the scaling function is related to the that of the low-pass filter. And in the same manner, the spectrum of the wavelet function is related to the that of the high-pass filter. Some more pointers are given in:

Why is a wavelet transform implemented as a filter bank?
Discrete Wavelet Transform (DWT) Filter Bank

The limit is that, if you start from continuous wavelets, very few of them can be faithfully implemented as discrete schemes as filters, and those filters, when they exist, don't look like wavelets. Because they don't have too. Their resemblance with wavelets only make sense when they are repeated over and over at different scales, across many dyadic sub-sampling.

And most discrete wavelets (like Daubechies') don't have an analytic formula. But if you look at the bottom line, from afar, the filters and the wavelets look similar in terms of envelope (red line below) and zero-crossings. So the global shape is not so different, as shown below with appropriate scaling of the wavelets and  filters:

This similarity seems common for smooth-enough wavelets. Honestly, I have not investigated this into more details, because the actual ressemblance  of filters and wavelets is  not so important, their properties are way more of interest.


Converting Audacity Filter Curve EQ into transfer function and applying it to a signal via python First of I am very new to Signal Processing and to python in general. I am trying to write a script where I would feed a voice recording into it, internally apply an eq and have the modified signal returned. Researching the topic I found some ways to process a signal via transfer functions in python, however I only found functions like a low-pass filter etc. which would not fit my needs as I have a specific eq manipulation in mind. I constructed the required filter curve eq in audacity:
If my limited understanding of signal processing is correct I would need to convert the continous curve into a discrete transfer function and then apply this transfer function to the input signal. However I have no clue how I would go about converting the curve into a Transfer-function. Furthermore, as it is the first time I am working with python I have very limited knowledge on how to implement this process in a python script. I would greatly appreciate any help you could provide me, thanks!

You can use the manage/export functionality to convert your Audacity filter preset into a text file. The file contains frequency/value pairs plus some side info. It looks like Audacity builds an 8192 FFT grid and uses a B-spline to interpolate the frequency/value pairs over the entire grid. I assume that just do an inverse FFT plus some circular time shifting to build a linear phase impulse response and than do something like fftfilt() to apply the filter.
You can just read the text file into python, parse it, and then follow the same steps.


Spectral window when FFT is not a power of 2 I am trying to better understand zero padding and the FFT. To do so, I am attempting to plot the spectrum of the Hann window function, i.e.

$$ w(n) = \frac{1}{2}\left(1-\cos\left(\frac{2\pi (n+1/2)}{N}\right)\right) $$

To accomplish this, I did the following: 1) I created a series of points $n=0,1,...,49$ that I substituted into the hamming window function to get $w(n)$, 2) I zero-padded the data to $M=512$ to get $x(n)$, 3) I used the FFT function in R to compute the Fourier transform of $x(n)$, 4) I calculated the magnitude squared of the Fourier transform obtained in part 3). This gave me the following result:

I then repeated the above steps, but with $M=500$ instead of $M=512$ (i.e. I zero-padded to 500 instead of 512). This gave me the following result:

Why are these two spectra so different? Everything that I have read about the FFT seems to say that I should choose $M$ to be a power of 2 in order to reduce computation time. But I haven't found anything that explains why the spectrum should change as it does.

Any help would be appreciated... thanks!

500 is an exact integer multiple of the period of the cosine function used in your raised cosine window.  Therefore there are zeros in the spectrum at all multiples of the cosine's frequency (due to those basis vectors being exactly orthogonal).  512 is not an integer multiple of that cosine's period, so none of the correlations against the DFT basis vectors of a 512 length FFT will be zero (just small).  Thus the difference between the two spectrums.


What kind of filter is that? Is it IIR? I am trying to answer the following question:

Is the system described by equation:

$$y[n]=0.5y[n-1]+x[n]-0.5x[n-1]$$

an IIR filter? My answer is yes.

Thank you

This is the FIR filter, although it looks like an IIR. If you calculate the coefficients you get finite impulse response:

$h=[1]$

This happens due to zero-pole cancellation:

$Y(z)-0.5Y(z)z^{-1}=X(z)-0.5X(z)z^{-1}$

$H(z)=\dfrac{Y(z)}{X(z)}=\dfrac{1-0.5z^{-1}}{1-0.5z^{-1}}=1 $

Yes, it can be tricky. Seeing $y[n-k]$ coefficients in LCCDE (Linear Constant Coefficients Difference Equation) doesn't necessarily mean it's an IIR filter. It might be just a recursive FIR filter.


Is the FFT magnitudes of multiple audio sources additive? Is the FFT magnitude of two audio signals when played together the sum of their individual FFT magnitudes across frequency bins? Note I'm talking only about the magnitude, ignoring phase.

For instance, if two pure signals of sine wave of the same frequency but 180 degrees off phase are played simultaneously, would the magnitude of the FFT at said frequency be double that of a single signal or 0?

I believe it's the latter but I want to confirm.

Thanks

The FFT algorithm efficiently computes the complex valued DFT frequency vector $X[k]$ for a given time domain sequence $x[n]$.

The magnitude of the FFT output $X[k]$ for a given bin $m$ is the magnitude of the complex number $X[m]$, so like this :

$$X[m] = X_r[m]  + j ~ X_i[m] \implies |X[m]| = \sqrt{ X_r^2[m] + X_i^2[m] }$$

Now assuming two signals $x_1[n]$ and $x_2[n]$ with their corresponding complex valued DFT vectors $X_1[k]$ and $X_2[k]$, we can see from the magnitudes of complex numbers that

$$| X_1[k] + X_2[k] | \neq | X_1[k] | + | X_2[k] | $$

hence the magnitude of FFT/DFT is not additive...


Mathematical justification for zero padding? This question asks what's the point of zero padding. The accepted answer is certainly very insightful, but I don't understand a big chunk of it:

  Zero padding allows one to use a longer FFT, which will produce a
  longer FFT result vector.

  A longer FFT result has more frequency bins that are more closely
  spaced in frequency. But they will be essentially providing the same
  result as a high quality Sinc interpolation of a shorter
  non-zero-padded FFT of the original data.

  This might result in a smoother looking spectrum when plotted without
  further interpolation.

  Although this interpolation won't help with resolving or the
  resolution of and/or between adjacent or nearby frequencies, it might
  make it easier to visually resolve the peak of a single isolated
  frequency that does not have any significant adjacent signals or noise
  in the spectrum. Statistically, the higher density of FFT result bins
  will probably make it more likely that the peak magnitude bin is
  closer to the frequency of a random isolated input frequency sinusoid,
  and without further interpolation (parabolic, et.al.).

What exactly is the meaning of "resolve" and "resolution" here, and how, mathematically, is it apparent that zero padding does not increase resolution.
Dually, how is apparent mathematically that zero padding means interpolation?

Resolution (in the context of spectral analysis) is the ability to distinguish two or more closely spaced sinusoidal components in a spectrum. If you can't resolve them, this means that you only see one maximum instead of two or more. Resolution is determined by the window length, i.e., by the number of time domain samples (given a certain sampling rate). Windowing will smear out or broaden any narrow band components in the spectrum. Remember that windowing, i.e., multiplication in the time domain, corresponds to convolution in the frequency domain, and this convolution with the spectrum of a standard window is a sort of averaging operation in the frequency domain, hence broadening any (frequency domain) impulses.

Resolution can only be improved by adding more data points, i.e., more information, to your analysis window. Zero-padding obviously doesn't add any information. What zero-padding does is increase the sampling density of the discretet-time Fourier transform (DTFT) of the windowed signal. For a length $N$ signal, the DTFT is given by

$$\tilde{X}(\omega)=\sum_{n=0}^{N-1}x[n]e^{-jn\omega}\tag{1}$$

The discrete Fourier transform is

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{2}$$

Comparing $(1)$ and $(2)$, we see that the DFT is just a sampled version of the DTFT:

$$X[k]=\tilde{X}\left(\frac{2\pi k}{N}\right)\tag{3}$$

Zero-padding increases the number of equidistant samples taken from the DTFT $\tilde{X}(\omega)$.

A detailed analysis of the fact that zero-padding corresponds to interpolation in the frequency domain is given in this answer.


envelope rms implementation by Matlab I am trying to implement the moving RMS by Matlab.
x = randn(50, 1);
xRMS = sqrt(movmean(x.^2, 21));
xRMSref = envelope(x, 21, 'rms');
plot(xRMSref,'DisplayName','xRMSref');hold on;plot(xRMS,'DisplayName','xRMS');hold off;
legend()

Why my estimation differs from Matlab's?

That is because envelope centers the signal before computing the standard deviation and afterwards corrects for the mean. Try this instead:
x = randn(50, 1);
xRMS = sqrt(movmean(x.^2, 21));
xRMS_centered = sqrt(movmean((x-mean(x)).^2, 21)) + mean(x);
xRMSref = envelope(x, 21, 'rms');
plot(xRMSref,'x-','DisplayName','xRMSref');hold on;plot(xRMS,'o-','DisplayName','xRMS');plot(xRMS_centered,'s-','DisplayName','xRMS centered');hold off;
legend()

*edit: since you were asking why envelope behaves this way: I think it's because it tries to return an envelope of the signal. It actually comes with two output arguments, a lower and an upper envelope with the idea that the signal should be between the two. Clearly if that's what we want, we need to account for the mean. See example in the following figure:

If you didn't center, you wouldn't capture the local variation around the mean accordingly.


Which (Fourier transform?) processing method to use for time variant audio processing? (newbie) Edit:
I've since been able to understand the how the relvant parts of Fourier transform work and how they relate to my problem / what I want to do.
And have come to the conclusion that Fourier transform is not the best suited to my problem and that I have an alternative way to do what I want that is much simpler to control / behave in the way I want in practice.

Since I'm new here I don't know what's the right thing to do with my question on this site. Just in case anybody was following this / was preparing an answer, I'm leaving this question up for one more day with this edit and after that I'll delete the question.

I've devloped a basic algorithm for processing audio. I have some programming experience but none with Fourier transforms (but will dedicate all of 2020 to learn this full time if need be).
My question relates to how best to program my algorithm in practice.

What I need to do is the following:

Analyze the frequency content at exact locations of a 44.1kHz PCM
audio stream / file.
Detect and modify the amplitude of specific frequencies at specific
time locations (I do not need to modify phase or have knowledge of
phase to do this).
Then turn this modified version back into audio without audible artifacts from processing errors/limitations.

As I said, I currenty do not know much about Fourier transform. So far my understanding is that an FFT (with windowing) will allow me to analyze the frequency spectrum of a block, detect and modify the amplitude of specific frequencies. But not modify the amplitude of specific frequencies at a specific time in that block and not the same frequency at a different time in that same block? For instance if for 44.1kHz sample rate and block size 8192 as far as I understand it normal FFT will not allow me to detect for instance two seperate 10kHz transients in that block as seperate and modify only one of them?

If my above understanding is correct. Can anybody suggest proper (Fourier transform like?) processing method(s) to use / investigate to achieve my goal?

I've throught of a rough way to do what I want with FFT/IFFT by dividing the frequency spectrum (5Hz to Nyquist for instance) into seperate octave bands all processed seperately, for instance first band 11.025kHz to 22.05kHz under 44.1kHz sample rate and block size 4 for instance not sure of the right block size just an estimate (so I get the transients in one block and still sufficient frequency resolution for that frequency band), then 5.5125kHz to 11.025kHz band under 22.05kHz sample rate (downsampled) also with block size 4 as an estimate, etc with a total of 12 bands down to about 5Hz for instance.
But this will be a lot of work and will take me a lot of time to develop and surely someone must have done something similar and smarter a long time ago to do the type of processing I wish to do :) If anybody can point me in the right direction (and perhaps if I got some assumptions wrong correct me) I will be very thankful!

I believe the answer to your question is Short-time Fourier Transform or Short-term Fourier Transform.  There is the wikipedia article.  I tried to show the essential math in this answer.

Both the Phase-Vocoder and Sinusoidal Modeling are done with the STFT of some form or another.  The two methods sorta merge in concept at the STFT level.  Remember that the analysis window for STFT need not be the same as the synthesis window in constructing the processed audio for output.


Channel capacity with additive Gaussian noise with changing variance As we all know the capacity for a AWGN channel with constant noise variance is just given by $$\frac{1}{2}\log(1+SNR)$$

What if the variance of the noise changes with the type (assuming the variance is known for any time index)? Consider the following discrete channel:
$$
y[n] = x[n] + w[n], \quad n=1,\dots,N
$$
where $w[n]$ is the noise during the $n$-th symbol distributed as $\mathcal{N}(0,\sigma^2_n)$, with known $\{\sigma_n^2\}_{n=1}^N$.

What is the capacity defined in such case?
It should be related to some dynamic rate control, right?
Could anyone give me some reference on related topic?

I think you are looking for the information rate,

$$ I(x_{1},x_{2},\cdots,x_{n}; y_{1}, y_{2}, \cdots, y_{n}) = \frac{1}{n} \sum_{i=1}^{n} I( x_{i};y_{i} ) $$

when input achieves the capacity of the channel $ x \sim N(0,P) $. Then,

$$ C = \frac{1}{2n} \sum_{i=1}^{n} \log \left( 1 + \text{SNR}_{i} \right) $$

I think every time we use the channel we have a new random variable $ y $ since we have a memoryless channel and $ w $ changes. Then, we have a sequence of related random variables $ \left\lbrace x_{i} \right\rbrace $ and $ \left\lbrace y_{i} \right\rbrace $. With the information rate we compute the mutual information between this two sequences.

A good reference about this Elements of Information Theory.

I hope this could help you.


Transients in marginally stable system I know that if a system is asymptotically stable then after infinite time its output  contains only characteristic of input signal and  it is also known as steady state and there are no transient component left in output, and if system is asymptotically unstable then whatever will be input is its output will become infinite after infinite time so we can say that it doesn't reach at steady state and it's transient components become infinitely large,
but what about transient and steady state in marginally stable systems?
Because in marginally stable systems even after infinite time it contains both inputs characteristics and systems characteristics so what exactly transient and steady state components will be?
I don't know whether I understand transient and steady state conditions properly
(whatever I explained above)  or it is just poorly defined in books?

I think your understanding is correct. The steady state response of an LTI system is the part of the response that is caused by a steady excitation at the input, like a sinusoidal signal, or any other periodic signal. The transient response is caused by changes at the input, like switching on a signal, or changing the parameters of a periodic input signal, such as amplitude, phase or frequency.
As an example of a marginally stable system, take an ideal integrator with impulse response $h(t)=u(t)$, where $u(t)$ is the unit step function. You can show that the response to a switched sinusoid at the input
$$x(t)=\sin(\omega_0t)u(t)\tag{1}$$
is
$$y(t)=\frac{1}{\omega_0}u(t)-\frac{1}{\omega_0}\cos(\omega_0t)u(t)\tag{2}$$
The first term on the right-hand side of $(2)$ is a non-decaying contribution caused by switching on the sinusoid at the input at $t=0$, and the second part is the steady state response to the sinusoidal input. The non-sinusoidal part of the output signal is what would be a transient if the system were a stable system, in which case it would decay. Since for this example it is non-decaying due to marginal stability, the term transient would of course be a misnomer.


How to detect and remove ringback-tone and ivrs voice etc from the beginning of an audio call recording In an audio recording (say a telephone conv. b/w two people), how would I programatically detect and remove the dial-tone at the beginning of a call using python.
Ex :
As you can see the first 15 seconds or so is just a dial tone like tring-tring-tring-tring.

Are there any audio analysis libraries in python that could  help me achieve this?

If this is not the right forum, kindly point me to the right place.

You can use librosa and scikit-learn to create a machine learning classifier. It would work roughly like this:

Training

Get training signals of (A) just phone ringing, and (B) no phone ringing, e.g. ordinary conversation.
Segment the training signals with a frame size of ~50-500 milliseconds.
Extract features from each frame, e.g. MFCCs.
Train a scikit-learn classifier, e.g.

classifier.fit(X, y)

where X is a ndarray of feature vectors, and y are the target labels, e.g. "ring" (1) and "no ring" (0).

Prediction

classifier.predict(X)

where X is an ndarray of feature vectors extracted in the same way from a test signal.

The latest frame which returns a positive "ring" label is where to truncate the signal.


What are the differences between these two Gabor filter functions I need to enhance the visibility of veins on dorsal hand vein images in my project. I use two different even-symmetric Gabor filters bank improve vein visibility.

First bank consists of these gabor functions:
$$G^\mathit{e}_\mathit{mk}(x,y)=\dfrac{\gamma}{2\pi\sigma^2}\exp\Bigg\{-\frac{1}{2}\left(\dfrac{x_\mathit{\theta}+\gamma^2y_\mathit{\theta}^2}{\sigma^2}\right)\Bigg\}\times \left(\cos(2\pi f_\mathit{0}x_\mathit{\theta})-\exp(-\dfrac{\upsilon^2}{2})\right)$$

Second bank consists of these ones:

$$G^\mathit{e}_\mathit{mk}(x,y)=\exp\Bigg\{-\frac{1}{2}\left(\dfrac{x_\mathit{\theta}+\gamma^2y_\mathit{\theta}^2}{\sigma^2}\right)\Bigg\}\times \cos(2\pi f_\mathit{0}x_\mathit{\theta})$$

where $m$ is the scale index, $k$ is the orientation index, $f_\theta$ is the filter center frequency, $\sigma$ is the standard deviation (often called scale), $\gamma$ is the aspect ratio of the elliptical Gaussian envelope, $\upsilon$ is the factor determining DC response, $x_\theta=(x\cos\theta+y\sin\theta)$ and $y_\theta=(-x\sin\theta+y\cos\theta)$ are rotated versions of the $x$ and $y$ coordinates.

I have coded these filters in MATLAB, I do not have any problem on coding. But I cannot understand the underlying difference between these two gabor functions.

Depending on the location of the peak and the scale of the two axes of the Gaussian envelope the filter may have a large DC response. A popular approach to get a zero DC response is to subtract the output of a low-pass Gaussian filter, which is what the first of these two does. In the case of images, if the DC response is not removed, the filter will respond to the absolute intensity of the image.

This tutorial gives a bit more detail.


Parks-McClellan - compiling C++ code from Iowa Hills page I am trying to compile C++ code for FIR filter design using Parks-McClellan algorithm from http://www.iowahills.com/Example%20Code/NewParksMcClellan.txt . However, I cannot find the following header files: "FIRFilterCode.h" and "Main.h" referenced in "NewParksMcClellan.cpp". Does anyone have experience compiling this code? Thank you.

You need to back track on that site.
There is a complete set of source code and header files available that you will need.

This is all about the FIR-Filters.

It links to this page with examples.
You need the "Code Kit Download."
This will get you the FIRFilterCode.h.
Main.h will be missing because that belongs to the GUI  - which they don't seem to provide code for.
However, it should be fairly simple to fix.
Remove the include Main.h, then remove all references to TopForm.
You'll probably want to replace them with some kind of logging since all they do is provide feedback to the user.
The involved lines are 273+274 and 278+279.
Notes to that effect are included in the comments in NewParksMcClellan.cpp


integral with vectoric varying coeficient i have a basic function exp(-x.^2).*log(x).^2 which i integrate in a certain interval.
i want to multiply my basic function with a vectoric coefficient called coef_vec that varies with the interval.
so if the integral is at 5 my basic function whould be multiplied with coef_vec(5).
i know that its som how turning the integral into a loop.
Is it possible in matlab?
Thanks.

coef_vec=linspace(1,10,100)
fun = @(x)coef_vec*exp(-x.^2).*log(x).^2;
q = integral(fun,1,10);

If you have a closed-form expression for the weighting function, you can use it  directly. Here, the weight is linear, so you can replace coef_vec by xand be ok.

fun = @(x)x.*exp(-x.^2).*log(x).^2;
q = integral(fun,1,10);

If not, an option is to define coef_vec on a larger interval, interpolate it with a suitable closed-form interpolating function, and use this under the integral.

Or, simply, evaluate the integral with a function like trapz. Here is a comparison, with close results

coef_vec=linspace(1,10,100);
fun0 = @(x)exp(-x.^2).*log(x).^2;
fun1 = @(x)x.*exp(-x.^2).*log(x).^2;
q1 = integral(fun1,1,10);

q2 = trapz(coef_vec,coef_vec.*fun0(coef_vec));
disp([q1 q2])

>>    0.0245    0.0245


Two type of calculating R in steepest-descent modeling algorithm I have wrote this algorithm for steepest descent:

function [e2,w_opt]=SteepestDescent_Alg2(x,d,itr,w)
N=length(w);
sample=length(x);
x=[zeros(N-1,1); x];%padding

x1=zeros(sample,N);
for i=1:N
    x1(:,i)=x(N-(i-1):end-(i-1));
end

R=x1'*x1/sample;
P=x1'*d/sample;
lambda=eig(R);
mu=1/(min(lambda)+max(lambda));

w_opt=zeros(size(w));
dd=(d.'*d)/(sample);
e2=zeros(1,itr);
for i=1:itr
   grad=2*R*w_opt-2*P;
   e2(i)=dd+w_opt'*R*w_opt-2*P'*w_opt;
   w_opt=w_opt-mu*grad;

end
wc=R\P
end

These lines forming R:

x1=zeros(sample,N);
for i=1:N
    x1(:,i)=x(N-(i-1):end-(i-1));
end
R=x1'*x1/sample;

But if we use other way to calculating R:

function R=corlnm2(h,N)
    %This function obtains the correlation matrix of the output
    %of an FIR filter which is derived by a white noise, and its
    %coefficients are given in vector "h".
    %"N" specifies the dimensions of "R".
    %
    %   R=corlnm2(h,N)
    %j=sqrt(-1);
    a=size(h);
    if a(1)>a(2)
       h=h';
    end
    m=length(h);

    for k=1:N
        if k<=m
            R(k)=h(k:m)*h(1:m-k+1)';
        else
            R(k)=0;
        end
    end
    R=toeplitz(R);
end

The code have written for this problem, also assume there is coloring filter with tap numerator coefficient of h:

Both two way of estimating R, getting value, really close together, but the second way gives weird learning curve.
Also the input noise set to 0.1 variance, thus minimum kesi(Mean Square Error) must be 0.1, calculating R in first manner give correct learning cure(plot of MSE vs itteration) but the W's are extremely close to Wo
even with more than 0.001 precision, why? And I think both way must give same answer why but this isn't ...?

So let's have a look at your calculations:

x1(:,i)=x(N-(i-1):end-(i-1));

Let's say N = 10 and sample = 100 so that end = 110.

For i=1 that means we sum from 10 to 110.

For i=2 that means we sum from 9 to 109.

etc.

Comparing that with:

R(k)=h(k:m)*h(1:m-k+1)';

where m is the same as sample:

For k=1 that means we sum from 1 to 100.

For k=2 that means we sum from 2 to 99.

I think all that's happening is your first calculation is assuming x is zero-padded and the second one is assuming that it's not.

That there are differences between the two is hardly surprising, given that they are summing different (possibly random) values.


Laplace Transform of Cosine, Poles and Mapping to Frequency Domain I am trying to understand the connection between Laplace transform ($s$-plane), and frequency domain calculation.

Let's take the Fourier transform of $\cos(\omega_0t)$, which equals to $\pi[\delta(\omega - \omega_0) + \delta(\omega + \omega_0)]$. So clearly the frequency domain has only two non-zero values at two particular frequencies, and others are zero. Fine!

Now lets do the Laplace transform of the same function: $\cos(\omega_0t)$, which gives me $F(s) = \frac{s}{s^2 + \omega_0^2}=\frac{s}{(s+j\omega_0)(s-j\omega_0)}$.

So there are two poles on the $y$-axis ($j\omega$ axis) as in above plot; $\beta$ is $\omega_0$ here. Now if I move upwards from the origin, I am staying on the frequency axis, and for every point except the pole locations, I am getting a non-zero value [from $F(s)]$, whereas only the pole locations should give me non-zero values (as can be seen from the fourier transform), and zero for all others points.

But this is not happening here. Where exactly I am going wrong here.

Your guidance will be greatly appreciated.

You're comparing the transforms of two different functions. You consider the Fourier transform of the function $x_1(t)=\cos(\omega_0 t)$, but you took the Laplace transform of the function $x_2(t)=\cos(\omega_0t)u(t)$, where $u(t)$ is the unit step function:

$$X_1(j\omega)=\int_{-\infty}^{\infty}x_1(t)e^{-j\omega t}dt\\
X_2(s)=\int_{0}^{\infty}x_2(t)e^{-st}dt=\int_{0}^{\infty}x_1(t)e^{-st}dt$$

Note the difference in the lower integration limits.

The (bilateral) Laplace transform of the function $x_1(t)=\cos(\omega_0t)$, $-\infty<t<\infty$, does not exist, whereas the Fourier transform of the function $x_2(t)=\cos(\omega_0t)u(t)$ does exist:

$$\mathcal{F}\{\cos(\omega_0t)u(t)\}=\frac{\pi}{2}[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)]+\frac{j\omega}{\omega_0^2-\omega^2}\tag{1}$$

However, it cannot be obtained by simply replacing $s$ by $j\omega$ in the expression for $X_2(s)$ because $X_2(s)$ has poles on the imaginary axis, and replacing $s$ by $j\omega$ only gives the correct expression for the Fourier transform if the imaginary axis is inside the region of convergence.

Note that in general for causal functions there are three cases concerning the relationship between the Laplace transform and the Fourier transform:

If the region of convergence (ROC) contains the $j\omega$-axis (i.e., all poles are in the left half-plane), the Fourier transform is simply obtained by the substitution $s=j\omega$.
If there are poles on the $j\omega$-axis, but no poles in the right half-plane, then the Fourier transform contains Dirac delta impulses (as in the example above) plus a term that can be obtained from the Laplace transform by setting $s=j\omega$. Note that in $(1)$ the right-most term is simply $X_2(j\omega)$.
If there are poles in the right half-plane, the Fourier transform of the corresponding causal signal does not exist.


Question about ASK, correlation, and frequency domain So I'm currently studying ASK (amplitude shift keying) and how it is modulated and demodulated. The modulation take out data signal, $M(t)$, and multiplies it with our carrier frequency, lets say $\cos(2\pi f t)$, and that's our modulated signal.

For the demodulation, we our taking our received signal and multiplying it with our carrier frequency, which in the frequency domain shifts our received signal to the left and right by our carrier frequency. So we end up with out signal both at the origin and at 2 times our carrier frequency. After that we can do a low pass filter to recover out message signal $M(t)$.

So my question is, what happened to our correlator? I thought in digital communications you take your received signal and pass it through a match filter, which can also be implemented as a correlator. In discrete time this is a multiplication with a signal and then a summation, right? What it seems like we're doing is just multiplying just a single point from the received signal with a single point of our carrier frequency. I've simulated this in MATLAB, and I can see what happens to the signal, I'm just confused why this simple multiplication be being called a correlator.

That being said, can this multiplication and low pass filter be done in the frequency domain by any chance? I know a convolution (or correlation) in the time domain is just a multiplication in the frequency domain, but i'm not sure how to do what we're doing for the ASK demodulation in the frequency domain.

Thanks for your time and input!

You are mixing different notions.

What happens in basic digital transmission is :

Digital signal generation (0's and 1's)
Bit to symbol association (digital modulation), in our example ASK.
Shape filtering, which allows the signal's spectrum to be contained within a standardized shape, since we can't just throw a digital signal with very sharp transitions in the air, in would pollute the electromagnetic spectrum.
Finally, we mix the signal with a carrier frequency in order to shifts it's spectrum around the appropriate transmitting frequency

What happens at the receptor is the same operations in anti-chronological order :

Analog demodulation (the signal is mixed with a carrier of the same frequency)
Matched filtering, which is the same operation as in the emitter, but with a specific filter, which is ideally the complex conjugate and temporally reverted of the filter used in emission ($g^*(-t)$). Why is it this one and not another ? See Nyquist ISI criterion.
Symbol estimation
Digital demodulation of estimed symbol

About the filtering operation specifically : time domain and frequency domain operations are equivalent. If you do a time domain filtering, or any mathematical operation (under trivial notations) :
$$
y(t) = g(t)*x(t)
$$
It is equivalent to frequency domain multiplication :

$$
Y(f) = G(f)X(f)
$$

Analog amplitude modulation of a signal $r(t)$ is, like you say, the operation $r(t)\cos(2\pi f_0t)$ (multiplication) which in frequency domain is equivalent to $R(f)(1/2)*(\delta(f-f_0)+\delta(f+f_0))$ (convolution).

Demodulation takes $r(t)\cos(2\pi f_0t)$ and multiply it again :
$$
s(t) = r(t)\cos(2\pi f_0t)\cos(2\pi f_0t)
$$
Which gives after simplifications :
$$
s(t) = r(t)\frac{1}{2}\left(1+\cos(4\pi f_0t)\right)
$$
That's when your LPF comes into action : it filters the unwanted $2f_0$ component and keeps only $r(t)$ which is the demodulated signal, and to do this mathematically, we write $s(t)*lpf(t)$, $lpf(t)$ being the LPF's impulse response.

In short : analog demodulation is either a multiplication in time domain ($r(t)\cos(2\pi f_0t)\cos(2\pi f_0t)$) followed by a convolution in time domain ($s(t)*lpf(t)$)

OR

It is a convolution in frequency domain ($R(f)*M(f)*M(f)$) (with $M(f)=(1/2)*(\delta(f-f_0)+\delta(f+f_0))$ followed by a multiplication in frequency domain ($S(f)LPF(f)$)

BUT both methods outputs are the same, physically speaking. We work in frequency or time domain according to which domain is best suited to give easy calculations, as you may have noticed.


Frequency-domain deconvolution: "Direct" filtering vs "Wiener" filtering Can someone help with clarifying the difference between two approaches
to frequency domain "deconvolution:
For the frequency domain problem:
We want to find a filter $F(\omega)$ which will convert some observed signal $B(\omega)$
into a desired signal $D(\omega)$, which is known
i.e. $D(\omega) = B(\omega)F(\omega)$

The intuitive / direct approach is to just calculate $F(\omega)=D(\omega)/B(\omega)$

Application of the least-squares criterion leads to $F(\omega)=R_{DB}(\omega)/R_{BB}(\omega)$ i.e. $D(\omega)B^*(\omega)/B(\omega)B^*(\omega)$

Q1:
If the direct division can be done successfully (i.e. no zeros in the denominator),
the two approaches seem to always give identical answers in practice.
Is that because the direct division is actually handled in code (e.g. Python) using the conjugates.
Q2:
Does the distinction between the two approaches only arise when the division
requires the addition of some stabilizing noise in the denominator?

The direct approach is noise sensitive and the second one (which is also known as $H_1$ estimator) is somehow noise resistant.
$H_1$ estimator assumes that there is no noise at the input signal, $B(\omega)$ in this case, and the noise exists only at output $D(\omega)$ and is uncorrelated to input $B(\omega)$. Therefore, the observed output signal is
$$
D'(\omega) = D(\omega)+N(\omega)
$$
where $N(\omega)$ is the uncorrelated noise. Under this assumption, the $H_1$ estimator gives a better result
$$
F'(\omega) = \frac{D'(\omega)B^*(\omega)}{B(\omega)B^*(\omega)} = \frac{D(\omega)B^*(\omega)+N(\omega)B^*(\omega)}{B(\omega)B^*(\omega)} = F(\omega) + \frac{N(\omega)B^*(\omega)}{B(\omega)B^*(\omega)}
$$
Since the noise is uncorrelated to the input signal, the second term can be eliminated by time averaging.

Here is some code and you can try some other SNR values. In practice for real data, you cannot obtain the clean signal without noise. All signal you measure is kind of noisy.
clear

% Input signal
x = mls;

% Comb filter (unknown system)
b = [1, zeros(1, 50), 1];
y = filter(b, 1, x);
fvtool(b, 1)

% Add white noise to the output
SNR = 10; % try some lower value and see what's happening
yn = awgn(y, SNR);

% H1 estimator
pxy = cpsd(x, yn);
pxx = cpsd(x, x);
H1 = pxy./pxx;
figure; subplot(211); plot(mag2db(abs(H1)))

% direct division
X = fft(x);
Y = fft(yn);
X = X(1:ceil(length(X)/2));
Y = Y(1:ceil(length(Y)/2));
H0 = Y./X;
subplot(212); plot(mag2db(abs(H0)))


How does window length influence cut-off frequency of filter and group delay? I'm attempting to interpret the following function. To my understanding, this aims to implement an FIR low-pass filter with a Nuttall window. Then, it filters the signal by a simple application of the convolution equivalence theorem.
static void GetFilteredSignal(int half_average_length, int fft_size,
    const fft_complex *y_spectrum, int y_length, double *filtered_signal) {
  double *low_pass_filter = new double[fft_size];
  // Nuttall window is used as a low-pass filter.
  // Cutoff frequency depends on the window length.
  NuttallWindow(half_average_length * 4, low_pass_filter);
  for (int i = half_average_length * 4; i < fft_size; ++i)
    low_pass_filter[i] = 0.0;

  fft_complex *low_pass_filter_spectrum = new fft_complex[fft_size];
  fft_plan forwardFFT = fft_plan_dft_r2c_1d(fft_size, low_pass_filter,
      low_pass_filter_spectrum, FFT_ESTIMATE);
  fft_execute(forwardFFT);

  // Convolution
  double tmp = y_spectrum[0][0] * low_pass_filter_spectrum[0][0] -
    y_spectrum[0][1] * low_pass_filter_spectrum[0][1];
  low_pass_filter_spectrum[0][1] =
    y_spectrum[0][0] * low_pass_filter_spectrum[0][1] +
    y_spectrum[0][1] * low_pass_filter_spectrum[0][0];
  low_pass_filter_spectrum[0][0] = tmp;
  for (int i = 1; i <= fft_size / 2; ++i) {
    tmp = y_spectrum[i][0] * low_pass_filter_spectrum[i][0] -
      y_spectrum[i][1] * low_pass_filter_spectrum[i][1];
    low_pass_filter_spectrum[i][1] =
      y_spectrum[i][0] * low_pass_filter_spectrum[i][1] +
      y_spectrum[i][1] * low_pass_filter_spectrum[i][0];
    low_pass_filter_spectrum[i][0] = tmp;
    low_pass_filter_spectrum[fft_size - i - 1][0] =
      low_pass_filter_spectrum[i][0];
    low_pass_filter_spectrum[fft_size - i - 1][1] =
      low_pass_filter_spectrum[i][1];
  }

  fft_plan inverseFFT = fft_plan_dft_c2r_1d(fft_size,
    low_pass_filter_spectrum, filtered_signal, FFT_ESTIMATE);
  fft_execute(inverseFFT);

  // Compensation of the delay.
  int index_bias = half_average_length * 2;
  for (int i = 0; i < y_length; ++i)
    filtered_signal[i] = filtered_signal[i + index_bias];

  fft_destroy_plan(inverseFFT);
  fft_destroy_plan(forwardFFT);
  delete[] low_pass_filter_spectrum;
  delete[] low_pass_filter;
}

There are two major points that I don't understand in this algorithm:

Half average length. What is it? It is used to calculate the filter order and the comment says it influences the cutoff frequency. I realise that the window length limits resolution when applying to a signal. My confusion is that all DSP guides I look up usually assume some filter order N as a design parameter, rather than selecting it based on a desired cut-off frequency. Is there a relationship between N and then cutoff frequency then in this case/generally?
Compensation of the delay. I'm not familiar with the details of filter design, but after reading up, I found that it is essential to shift a time-domain signal after FIR filtering. Filtering usually has frequency-dependent group delay, but FIR filters have the advantage of a constant group delay being $ \frac{N-1}{2} $, where N is the number of taps or order of the filter. Assuming $ \text{half_average_length} \cdot 4 = N $, it follows that $ \text{index_bias} = (\text{half_average_length} \cdot 4 - 1) /2 $ . I realise this wouldn't make sense as it's not a whole number, but then there is something wrong with my interpretation of $ \frac{N-1}{2} $. What is wrong with my interpretation of this formula?

Half average length. What is it?

Primarily something the author made up. They design a lowpass more or less as a weighted moving average filter. They simply use the window itself as the impulse response of the filter. That's NOT the same as designing an ideal low pass filter at the cutoff frequency you want and then windowing it down to the filter order. This seems to work reasonably well, you get frequency response of a nuttall window (which isn't bad) and the cut off frequency is roughly the sample rate divided by the window length.
Ir your sample rate is 44.1kHz and your window length is, say, N=32, the cut off will be around 1320 Hz or thereabouts.

Compensation of the delay. .... What is wrong with my interpretation of this formula?

Your interpretation is correct. If the window length is even, you can't correct for the group delay by simply shifting samples. The code only "approximately" compensates for group delay and is off by half a sample.
There are a few more things wrong about the code:

If someone calls it with half_average > fft_size/4 , bad things will happen
Multiplication in the frequency domain is circular convolution (not linear convolution) in the time domain. The input must be properly zero-padded before FFT to avoid time domain aliasing.


How is the ends of a bandpass filtered signal reliable after applying the cosine taper? A taper function decaying to zero at ends such as a cosine window is routinely applied to the signal before bandpass filtering to reduce end effects. Obviously the window alters the waveform.

So, as for the filtered waveform, how is it reliable near the beginning and end?

how does the the decay ratio of the window (e.g. 5%, 10% ...) influence? It's selected arbitrarily or by some criteria?

In the filtering functions of matlab and scipy, some extrapolation is applied to the ends of the signal to match their slope. In such situations, is the windowing still necessary?

A taper function is applied for spectral estimation if the artifacts from a non-tapered rectangular window provide less useful results that the artifacts due to the tapered window.

Not tapering a finite length sub-sample vector is equivalent to rectangular windowing, which convolves Sinc shaped artifacts onto the spectrum.


How to convert the mean and variance of a processed received signal into a SNR or BER? So I'm trying to eventually get some sort of signal quality metric, and to do that, I'm trying to take the mean and variance of my bit determination signal and use that to come up with a Signal to Noise Ratio, but I think I'm missing something. What I'd like to do is take my calculated SNR value, determine the Probability of error, $P_e$,  from this, and then transmit a bunch of known data, and at the receiver see how many bits actually are received in error to see if it matches with my calculated $P_e$.

What's going on is when we get our signal, we multiply it with the previous bit to get a correlation, and from that we sum it with the correlation from the previous bit to get our detection signal. This detection signal is what I'm taking the mean and variance of. If the mid-bit point of this detection signal is above a threshold (0 for now) then we say the bit is a 1, if it's below the threshold, we say it's a 0. I'm taking the mean and variance of all the 1's and the mean and variance of all the 0's, and am trying to come up with an average SNR type quantity. My thought was I could just say $\mu = Mean$, $\sigma ^2 - variance$ so say my $SNR = \frac{\mu}{\sigma^2} $ since the mean of this decision signal should be directly related to my signal power, and the variance of these points (which should relate to the spread of my decision points for 1's and 0's) will be directly related to my noise.

When transmitting Binary symbols, I know that if I increase my data rate, then the power per bit will decrease, and I see this when calculating $\mu$ and $\sigma^2$. For instance the $\mu$ is cut in half if I double the bit rate, but $\sigma^2$ stays the same, which I'm interpreting as the power/bit is being halved, but the noise is staying the same so the mean should be half and the variance of the noise should be roughly constant.

Our received signal is received in Volts, and then put through and ADC so I realize there is probably something to take into consideration there... I'm just not quite sure what yet. Am I thinking about all of this correctly?

I also feel that the question needs a little bit more clarity but with what is stated so far, I think I can have a go at a response and edit it later if more information is added to the question.

  What I'd like to do is...

I do not intend to 'put you on the spot' with my question here: Have you checked how other people have tackled this? It seems that what you are trying to do is get an estimate of received signal quality (?). One of the simplest ways to achieve this is Received Signal Strength Indication (RSSI).

The estimation of both Signal To Noise Ratio (SNR) and Bit Error Ratio (BER) imply the use of a 'reference signal'. We will come back to this.

  What's going on is when we get our signal, we multiply it with the previous bit to get a correlation, and from that we sum it with the correlation from the previous bit to get our detection signal. This detection signal is what I'm taking the mean and variance of. If the mid-bit point of this detection signal is above a threshold (0 for now) then we say the bit is a 1, if it's below the threshold, we say it's a 0. I'm taking the mean and variance of all the 1's and the mean and variance of all the 0's, and am trying to come up with an average SNR type quantity.

This sounds like a differential codec (e.g. this one). Two questions emerge naturally here: The first one is, how do you know that the bit your are comparing against was decoded correctly, in order to re-use it to guide the decoding of your next bit (if I am getting the scheme correctly here) and the second question takes us back to RSSI and is "If you are trying to estimate statistics of the received signal timeseries with the objective of estimating the SNR, then why don't you use the RSSI?".

OK, so, let's see how are we going to do this. What is Signal to Noise Ratio?. It is the ratio between the strength of a "desired" signal versus the strength of an "undesired" signal (obviously, in the same context!). What is our "desired" signal here? It is the output of the modulator which is coupled to our medium somehow (antenna, speaker, light, whatever). What is the "undesired" signal here? It is the level of background "noise", undesired perturbations depending on the medium, which distort reception.

So, a way of guessing how bad our reception is going to be would be to "listen" to the background noise when there is NO reception (for some time) and create the N and then compare this with the strength of the signal S when we know that the transmitter is transmitting. Here, we form a ratio between a reference (our background noise) and a signal (received signal strength). The closer these two numbers are, the harder it is to accurately demodulate and eventually decode the signal.

If you absolutely must derive an SNR and a BER estimate, then you can examine the use of pilot signals or sequence. Communication systems that use a carrier signal, usually have a Phase Locked Loop (PLL) to lock on the carrier and track it. To allow for the PLL to achieve "lock", prior to transmitting the actual "information" signal, you can transmit a "dummy" sequence of some length, for the PLL to lock (e.g. listen to this clip). After the PLL "locks", the local oscillator (of the demodulator) is synchronised with the remote oscillator (of the modulator, at the transmitter side) and we can now start demodulating symbols. BUT we now do not know, in this bitstream, where does the actual information starts and where does it end. The pilot sequence may have been 100 symbols long. Maybe the PLL locked after 4 cycles, in which case, we have to throw away 96 "dummy" symbols, but maybe it locked after 50 symbols, in which case we have to throw away 50 "dummy" symbols. The point is that we don't know this 4 or 50. All we know is a line coming out of the PLL that is raised to tell us "I have now achieved lock".

So what we do, after the pilot tone, is to insert a start sequence. For example 01101001100101100. Once the system detects that, it knows that it can start decoding the actual information.

So, what you can do, is insert a long pilot sequence the structure of which you know (reference signal). For example 10101010101010101010101010101010101010101010101010101[BLOCK_START_SEQUENCE]. Once your PLL achieves lock, you don't throw away the symbols but you store them until you find the start sequence at which point you stop. Now, the pilot sequence could be 100 symbols long, the PLL achieved lock after 20 which means that we have decoded 80 symbols until the block start sequence. You know that these 80 symbols should be 10101010101010101010 (reference signal) and you count how many of them have been decoded in the right way. For example, 80 symbols decoded, 40 of them decoded right, that is a bit error rate of 0.5. For more information, you might want to see this paper.

Presumably, you can then use the BER in some adaptive scheme to select a strong or weak code. In addition, you might also want to look at equalisation and/or adaptive coding

Hope this helps.


Fast Fourier Transform using numpy I am a computer science student and didn't really have signal processing as a subject. Maybe I should be clear on the concepts of sampling rate and frequency of the signal but I am a little confused.
I was trying to code and understand the process of finding MFCC filter banks on audio files. One intermediate step is to find FFT. I computed and plot the fft results of my audio but wasn`t sure that I was doing things correct as I didnt know the component frequencies that may have been present on my audio.

So I decided to form a sample wave and find and plot the fft results of the test signal.(The np.fft.fftfreq functions return the frequencies corresponding to the fft computed by np.fft.fft()).So my questions are

1) I found out value of y for time at a separation of 1ms seconds( 0 to 1, 1000 values). So my sampling rate should be 1000 right?. Then according to the np.fft.fftfreq() function, the second parameter should be 1/sampling rate. Then it gives correct values for 1/1000. However if I give sampling rate as 100(T=1/100) and give the second parameter as 1/100. It gives wrong results. Why??

2) If I give the frequency of the sine wave as 1080 instead of 80 for T=1/1000. Then also it gives same results as for 80. Why is it so?

# EXAMPLE SIGNAL TAKEN TO VERIFY ABOVE FFT PROCEDURE IS RIGHT
NUM=1000
T=1/1000
tt=np.linspace(0,NUM*T,num=NUM)
#So I take readings every T seconds = 1ms

print(1/T)
y=np.sin(2*np.pi*1080*tt)#+np.sin(2*np.pi*90*tt)
#plt.plot(tt,y)
NFFT=1024
mag_frames = np.absolute(np.fft.rfft(y, NFFT))
xf = np.linspace(0.0, 1.0/(2.0*T), NUM/2)
freq=np.fft.rfftfreq(NFFT,1/1000)
plt.plot(freq,mag_frames)

3) What should be value for the second parameter when I perform fft on a audio data imported from a .wav file by librosa library with sampling rate = rate

I am sorry if the problem is very simple or basic. I am really a beginner in this area, and could not get much help from results in Google search

Answers to Your Questions

You are right that your time instants are spaced at an interval of $T$. Hence, the sampling rate is $1/T$. ($T$ is defined in line 3 of your code.) You are also right that the second parameter of np.fft.fftfreq() is the sampling interval $T$. I am, however, not sure why you get wrong results when setting $T=100\,\text{Hz}$.
This phenomenon is called aliasing. Nyquist's sampling theorem states that frequencies below $f_\text{s}/2=\frac{1}{2T}$ can be captured correctly by the sampling process. Frequencies above half of the sampling frequency are mirrored back to frequencies below half of the frequency. Therefore, a frequency of $1080\text{Hz}$ is too high to be captured correctly.
The choice for the second parameter of the FFT, the FFT length, depends on the specific application. You have to experiment what leads to the best results. Still, your choice of NFFT = 1024 is typically a good start.


Implementing a Band-Pass Filter from Scratch I am trying to implement a band-pass filter from scratch. Here is an article which I thought explained the nuts and bolts of how to build one; the author combines a low-pass with a high-pass filter (convolving both) to create a band-pass filter. I applied it to a mixture of sine waves, but the result was not satisfactory.

import numpy as np
from numpy import sin, linspace, pi
from pylab import figure, plot, show, title, xlabel, ylabel, subplot, savefig, xlim, ylim
from scipy import fft, arange

def plotSpectrum(y,Fs):
    n = len(y) # length of the signal
    k = arange(n)
    T = n/Fs
    frq = k/T # two sides frequency range
    frq = frq[range(n/2)] # one side frequency range

    Y = fft(y)/n # fft computing and normalization
    Y = Y[range(n/2)]

    plot(frq,abs(Y),'r') # plotting the spectrum
    xlabel('Freq (Hz)')
    ylabel('|Y(freq)|')

Fs = 150.0;  # sampling rate
Ts = 1.0/Fs; # sampling interval
t = arange(0,30,Ts) # time vector

ff = 0.2;   # frequency of the signal
y = sin(2*pi*ff*t)
ff = 0.8;   # frequency of the signal
y = y + sin(2*pi*ff*t)

figure()
subplot(2,1,1)
plot(t,y)
xlabel('Time')
ylabel('Amplitude')
subplot(2,1,2)
plotSpectrum(y,Fs)
savefig('out1.png')

fL = 0.1
fH = 0.3
N = 33
n = np.arange(N)

hlpf = np.sinc(2 * fH * (n - (N - 1) / 2.))
hlpf *= np.blackman(N)
hlpf = hlpf / np.sum(hlpf)

# Compute a high-pass filter with cutoff frequency fL.
hhpf = np.sinc(2 * fL * (n - (N - 1) / 2.))
hhpf *= np.blackman(N)
hhpf = hhpf / np.sum(hhpf)
hhpf = -hhpf
hhpf[(N - 1) / 2] += 1

hp = np.convolve(hlpf, hhpf)
filtered = np.convolve(y, hp)

t2 = np.linspace(0,30.,len(filtered)) # time vector
figure()
subplot(2,1,1)
plot(t2, filtered)
xlabel('Time')
ylabel('Amplitude')
subplot(2,1,2)
plotSpectrum(filtered,Fs)
savefig('out2.png')

For the original signal:

I get an output like this:

To confirm I could filter this signal, I used a butterworth filter,

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

filtered = butter_bandpass_filter(y, fL, fH, Fs, order=2)

That worked fine (although I had to play around with the order a little bit, higher orders of this filter does not work either). What is wrong with the first filter?

Regards,

A simple way to make a real bandpass filter (assuming you know how to make a low pass filter) is multiply the coefficients of your prototype low pass filter by $cos(2\pi f_c nT)$ where $f_c$ is the desired center frequency of your passband and T is the sampling interval for your coefficients.


Effect of low pass filter on the time domain Consider an ideal (rectangular profile in frequency domain) low-pass filter (LPF) which simply sets to zero all the frequencies beyond a certain cut-off frequency. Let $x$ be a discrete-time signal of duration $N$ in time domain, i.e., $x=[x_1,x_2,\ldots,x_N]$. Apply $x$ as input to the LPF and let $y=[y_1,y_2,\ldots,y_N]$ be the corresponding output signal in time domain. I have checked the spectrum of input and output signals and I can confirm the filter works well and also the input and output signals make sense.

From simulations, I can see that the following condition holds:
$$
\sum_{t=1}^{N} x_t = \sum_{t=1}^{N} y_t.
$$
My question is why this condition should hold. Is it there a general rule? I can think of an example where it holds. If I filter every frequency except the DC component, then the output signal in time is the mean value of the input signal and indeed the condition above is verified since
$$
\frac{\sum_{t=1}^{N}x_t}{N}=\bar{x},
$$
where $\bar{x}$ denotes the mean value of signal $x$. If $y_t=\bar{x}$, $\forall t$, then
$$
\sum_{t=1}^{N}y_t = \sum_{t=1}^{N}\bar{x} = N \bar{x} = \sum_{t=1}^{N}x_t.
$$

For example, I know that by Parseval's Theorem the energy of the output signal will be different from the energy of the input signal (and this occurs) but I am not sure why the condition above should hold.
Thank you very much in advance.

If you sample a pure sinewave for an integer number of periods, the samples will sum to zero (it will have equal "ups" and "downs").  If you zero an FFT bin, that's the same as subtracting a sinewave of the same magnitude and phase as was in that bin, or identically, adding a pure sinewave of the same magnitude but opposite phase.  So zeroing bins is the same as adding sinewaves that are exactly periodic in FFT aperture (as those are the basis vectors of a DFT).  So zeroing any set of bins (other than the DC or 0 Hz bin) adds samples than sum to zero, thus having no effect on the total sum of the time domain samples after the modification.


Bandwidth of a trimmed sinus wave I generate a sinus which duration is 15ms and frequeuncy of 13KHz. Is its bandwidth is (1/15ms) ? Or in other words, BW = 66.6Hz?
I ask that because I have two freqs in my system: 13K & 13.1K.
I know that they are seperated by 100Hz. And as the tone I generate have a lower duration I think it spreads into its friend area.

The spectrum of your windowed sinusoid is just a shifted sinc function (assuming a rectangular window). If $T=0.015s$ is the length of the signal then the magnitude of the spectrum at positive frequencies is given by

$$S(f)=T|\text{sinc}((f-f_0)T)|,\quad f>0$$

with $f_0$ the frequency of the sinusoid, and where I use the definition $\text{sinc}(x)=\sin(\pi x)/(\pi x)$. So with $\Delta f=100Hz$ (the distance to the other sinusoid), you get a relative attenuation of

$$|\text{sinc}(\Delta f\cdot T)|=0.21221\doteq 13.5\,\text{dB}$$

So at the frequency of the neighboring sinusoid, the spectrum of the other (finite length) sinusoid is about $13.5\,\text{dB}$ lower.


Normalize Fourier Transform For Video Game I am in the process of developing a video game where I would like to sync a character's mouth with an audio clip of their speech that is playing.  I found this great website that shows the profiles of different vowel sounds in a frequency spectrum: http://hyperphysics.phy-astr.gsu.edu/hbase/music/vowel.html

After a little (actually, a lot of) effort, I am able to show in real time the spectrum of the audio that is playing.

The next step I am planning to do is find the 4 largest peaks and try to match it to one of the vowel profiles from the website above.  However, I have noticed that the relative power of the frequency domain coefficients is always strong lower than higher.  The picture I showed is in silence. During an SSSSSS sound, the spectrum does show peaks at higher frequencies:

If I ignore the first few points in the silence FT and then relatively graph the rest, the spectrum has the same shape.  If I ignore a few more, again the lower frequencies are much stronger than the higher frequencies.

I assume that this is the result of the math which I can barely get through to code the basic FFT.  Should I be normalizing the result somehow to be able to see peaks in the higher frequencies?  What process should I use for this?

Thanks in advance for any help you can provide!

To add info requested in the comments, by "in real time", I mean that every frame (at about 60 frames per seconds), I get a set of 1024 samples (sampled at 44100 Hz) that is centered at the point of audio that is currently playing).  I then compute the FT of these samples and draw the spectrum.

Not sure if your approach is the right path. Anyhow, hearing aid research has produced many "long term average speech spectrum" (LTASS) curves that you can use to whiten the spectrum by dividing it with the LTASS. Have a look at http://www.dydaktyka.cba.pl/AvSpeech.pdf (Byrne et al. 1994 "An international comparison of long‐term average speech spectra" J. Acoust. Soc. Am. 96, 2108) Just don't divide by the LTASS dB values but linearize them first to $10^{LTASS\_in\_dB/20}$ if you are dealing with magnitude frequency spectrum or to $10^{LTASS\_in\_dB/10}$ if you are dealing with power spectrum.

Here's the worldwide average LTASS from Byrne et al. 1994:

Freq.   Male    Female  Combined
(Hz)    (dB)    (dB)    (dB)
63      38.6    37.0    38.6
80      43.5    36.0    43.5
100     54.4    37.5    54.4
125     57.7    40.1    57.7
160     56.8    53.4    56.8
200     58.2    62.2    60.2
250     59.7    60.9    60.3
315     60.0    58.1    59.0
400     62.4    61.7    62.1
500     62.6    61.7    62.1
630     60.6    60.4    60.5
800     55.7    58.0    56.8
1000    53.1    54.3    53.7
1250    53.7    52.3    53.0
1600    52.3    51.7    52.0
2000    48.7    48.8    48.7
2500    48.9    47.3    48.1
3150    47.0    46.7    46.8
4000    46.0    45.3    45.6
5000    44.4    44.6    44.5
6300    43.3    45.2    44.3
8000    42.4    44.9    43.7
10000   41.9    45.0    43.4
12500   39.8    42.8    41.3
16000   40.4    41.1    40.7


Matching two very similar images I have two images for the same object from the almost same perspective (no more than 1-2 degrees or 1-2 mm is acceptable).
I want to change the perspective of one of them to match exactly the other so I can make subtraction between them successfully to get a complete black result if there is no change of the content of the object.
many thanks

I think all these answers fail to provide a descent methodology.

You should consider "the closeness" of the images as the availability of the initial pose and cast this as a registration problem. The solution can typically be obtained through a least-squares adjustment, i.e., by minimizing the distances of the model points (from first image) to their corresponding image points. In case of a perspective distortion, you could optimize the 9 parameters of the homography matrix. Methods like Lukas Kanade, ESM or hyperplane tracking are some ways to do that. In essence, the procedure is similar to ICP, where an iterative procedure is applied.

Here is one work, extending ESM and comparing a few methods. If your camera is calibrated, you could benefit from some alternative methods.

Note that, to begin with, one could always resort to a simple approach, where you try different combinations of rigid motion $[-2,2]$ degrees for rotation and $[-5,5]$ pixels for translation. Then transform your template with all possible poses (in this case ~$40-50$). Then interpolate the sum of squared differences or normalized cross correlation scores in the pose space (6DOF freedom or 9DOF homography) to get a subpixel accurate refinement. This should give you a rough yet a better alignment. Keep in mind that you could also use this scheme iteratively.

If you have absolutely no clue on how to do any of these, check out this link.


The difference between MSK as CPFSK and OQPSK I am working on MSK simulation. I have read the Chapter 5 of book entitled Digital Modulation Techniques by Fuqin Xiong. It is said that MSK can be viewed as a Sinusoidal Weighted OQPSK or CPFSK. If we write the MSK signal as

$$
s(t)=A\cos\left(2\pi f_c t + d_k \frac{\pi t}{2T} + \phi_k\right), \quad kT\leq t\leq (k+1)T
$$

From the first perspective (MSK as a Sinusoidal Weighted OQPSK), we have
$$
d_k=I_kQ_k , \quad \phi_k=\frac{\pi}{2}(1-I_k)
$$

From the second perspective (MSK as CPFSK), we have
$$
d_k=d_{k-1}\newcommand*\xor{\mathbin{\oplus}}a_k,  \quad \phi_k=\frac{\pi}{2}\left(\sum_{i=0}^{k-1}d_i-kd_k\right)
$$
where $\{a_k\}$ represents the original data stream and $\newcommand*\xor{\mathbin{\oplus}}$ is XOR. I have the following question:

Are the two methods for computing MSK and specifically $\phi_k$ equivalent?
Do they result in the same MSK signal?

In my simulations, the two MSK signals that I get from these two perspectives don't have the same phase.

The answer to your question is yes, these two representations are equivalent. The first important thing to note is that in your formula, $d_k$ must be either $1$ or $-1$, because the frequency difference for MSK must be $\Delta f=1/2T$.

I'll add some information showing the equivalence of the two representations. Expanding your first formula gives

$$\begin{align}s(t)&=A\cos\left(2\pi f_ct+d_k\frac{\pi t}{2T}+\phi_k\right)\\&=A\left[\cos(2\pi f_ct)\cos\left(d_k\frac{\pi t}{2T}+\phi_k\right)-\sin(2\pi f_ct)\sin\left(d_k\frac{\pi t}{2T}+\phi_k\right)\right]\\&=A\left[\cos(2\pi f_ct)\cos\left(d_k\frac{\pi t}{2T}\right)\cos(\phi_k)-\sin(2\pi f_ct)\sin\left(d_k\frac{\pi t}{2T}\right)\cos(\phi_k)\right]\\&=A\left[\cos(2\pi f_ct)\cos\left(\frac{\pi t}{2T}\right)\cos(\phi_k)-\sin(2\pi f_ct)\sin\left(\frac{\pi t}{2T}\right)d_k\cos(\phi_k)\right]\\&=A\left[\cos(2\pi f_ct)\cos\left(\frac{\pi t}{2T}\right)I_k-\sin(2\pi f_ct)\sin\left(\frac{\pi t}{2T}\right)Q_k\right]\tag{1}\end{align}$$

with

$$I_k=\cos(\phi_k)\quad\text{and}\quad Q_k=d_k\cos(\phi_k)\tag{2}$$

where I've used $\cos(a+b)=\cos(a)\cos(b)-\sin(a)\sin(b)$, $\sin(a+b)=\cos(a)\sin(b)+\sin(a)\cos(b)$, the fact that $\sin(\phi_k)=0$ (because $\phi_k\in\{0,\pi\})$, and the fact that $d_k\in\{1,-1\}$.

Equation $(1)$ shows that the MSK signal $s(t)$ can be represented as a QPSK signal with in-phase component $I_k$ and quadrature component $Q_k$ with sinusoidal pulse shaping. In this blog post, it is shown that $I_k$ and $Q_k$ can never change at the same time, and both of them change with a rate $1/2T$. That's why this type of modulation is called offset QPSK (OQPSK). All this is related to the phase recursion shown in Eq. $(3)$ below.

In order for the MSK signal to have a continuous phase, the phase $\phi_k$ must satisfy the following recurrence relation:

$$\phi_k=\phi_{k-1}+(d_{k-1}-d_k)\frac{k\pi}{2}\qquad(\text{mod } 2\pi)\tag{3}$$

This relation is also derived in the above mentioned post (Eq. $(7)$). Note that the formula for $\phi_k$ in your question satisfies $(3)$. However, it's probably easier to directly use the recursion $(3)$. From $(3)$ it follows that

$$\phi_k=\begin{cases}\phi_{k-1}\pm\pi,&k\text{ odd and }d_k\neq d_{k-1}\\\phi_{k-1},&\text{otherwise}\end{cases}\tag{4}$$

From $(4)$ it is clear that $\phi_k$ changes with a rate $1/2T$.

Finally, inverting $(2)$, i.e., expressing $d_k$ and $\phi_k$ in terms of $I_k$ and $Q_k$, is straightforward:

$$I_kQ_k=d_k\cos^2(\phi_k)=d_k\tag{5}$$

because $\cos(\phi_k)$ can only be $1$ or $-1$. Furthermore, $1-\cos(\phi_k)$ equals $0$ for $\phi_k=0$ and $2$ for $\phi_k=\pi$. Consequently,

$$\phi_k=\frac{\pi}{2}(1-I_k)\tag{6}$$

Both Equations $(5)$ and $(6)$ are correctly stated in your question.

In sum, all your formulas appear to be correct. You haven't shown your formulas for $I_k$ and $Q_k$; they should be equivalent to Eq. $(2)$. I have the suspicion that your $d_k$ are either $1$ or $0$. As mentioned above, they should be either $1$ or $-1$, otherwise $s(t)$ is no valid MSK signal. Also note that the mapping between the information bits and $d_k$ is totally irrelevant for all derivations shown above.

EDIT: In analogy with Eq. $(4)$, you can derive recursions for $I_k$ and $Q_k$ from Eqs. $(2)$ and $(3)$:

$$I_k=\begin{cases}-I_{k-1},&k\text{ odd and }d_k\neq d_{k-1}\\I_{k-1},&\text{otherwise}\end{cases}\tag{7}$$

$$Q_k=\begin{cases}-Q_{k-1},&k\text{ even and }d_k\neq d_{k-1}\\Q_{k-1},&\text{otherwise}\end{cases}\tag{8}$$

Choosing an initial phase $\phi_0=0$ results in the initial values $I_0=1$ and $Q_0=d_0$.
From $(7)$ and $(8)$ it is obvious that $I_k$ and $Q_k$ never change at the same time, hence offset QPSK.


Determining Power Spectral Density from Autocorrelation Function and Time-Series Data So I think I've confused myself a bit here. I know that taking the FFT of time-series data (wind velocities, for example) yields a power spectral density. However, I just read that taking the Fourier Transform of an autocorrelation function also yields the power spectral density. So, if one were to take the inverse Fourier Transform of power spectral density, would you get the original time series data (with the mean subtracted out) or would you get the autocorrelation function?

Taking the FFT of a time-series does not give you the power spectral density (PSD). What you can do is take the squared magnitude of the FFT to get an estimate of the PSD. And the latter is equivalent to taking the Fourier transform of the (deterministic) autocorrelation function of the data. See also this answer to a related question.

So the inverse Fourier transform of the PSD is the autocorrelation function.


Unit impulse response $h[n]$ of a discrete time system with multiple characteristic roots 0 I am trying to obtain the unit impulse response of a system in the form:
$$y[n+N]+...+a_{N-1}y[n+1]+a_{N}y[n]=b_0x[n+N]+...+b_{N-1}x[n+1]+b_{N}x[n]$$
$$Q[E]y[n]=P[E]x[n]$$
(where E is the unit advance operator)
Using the following formula:
$$h[n]=\frac{b_N}{a_N}\delta[n]+y_c[n]u[n]$$
Or if $a_N=0$
$$h[n]=A_0\delta[n]+A_1\delta[n-1]+y_c[n]u[n]$$
Where $A_0$ and $A_1$ are determined through iterative calculation of initial values of $h[n]$ and $y_c[n]$ is a combination of the systems characteristic roots that satisfies $Q[E]y_c[n]u[n]=0$
The system in question is:  $y[n]=0.5y[n-1]+x[n]-(0.5)^4x[n-4]$
Upon converting the system to using the operator E notation to obtain the characteristic roots, I obtain 4 roots 3 of which are repeated zeros.
$$(E^4-0.5E^3)y[n]=(E^4-0.5^4)x[n]$$
$$\gamma_{1,2,3} = 0, \gamma_4=0.5$$
Using the relations described above and by simply discarding the repeated zero roots I then obtain the answer:
$$h[n] = (0.5)^nu[n]$$
However, when replacing $x[n]$ and $y[n]$ with $\delta[n]$ and $h[n]$ respectively in the original system difference equation, I obtain:
$$h[n] = (0.5)^n(u[n]-u[n-4])$$
Which is matching with my model answer.
I figured that maybe there is a mistake in my original solution due to the discarding of the repeated zero roots. How should I approach solving systems with zero repeated roots?

If you have a difference equation
$$\sum_{k=0}^Na_ky[n-k]=\sum_{k=0}^Mb_kx[n-k]\tag{1}$$
and if $M>N$, the general expression for the corresponding causal impulse response is
$$h[n]=A_0\delta[n]+A_1\delta[n-1]+\ldots+A_{M-N}\delta[n-(M-N)]+u[n]\sum_{k}c_k\gamma_k^n\tag{2}$$
where $\gamma_k$ are the roots of the characteristic equation, and where I've assumed simple roots. Note that roots $\gamma_k=0$ do not contribute to $(2)$ (even if there are multiple roots equal to zero).
For the given example we have $M-N=3$ and one non-zero root $\gamma=\frac12$. Hence the solution has the form
$$h[n]=A_0\delta[n]+A_1\delta[n-1]+A_2\delta[n-2]+A_3\delta[n-3]+c\left(\frac12\right)^n\tag{3}$$
The constants $A_i$, $i=0,1,\ldots,3$, and $c$ must be determined by iterating through the difference equation with $x[n]=\delta[n]$, $h[n]=y[n]$, and $h[-1]=0$.
If you do that you'll find that $c=0$, which means that the given difference equation describes an FIR system. In terms of transfer functions this means that there is a pole-zero cancellation, leaving only poles at the origin of the complex plane.


Is Image Noise Removal Image Enhancement or Image Restoration Operation? I am working in digital image processing field, recently I am studying about Image Noise, I just want to know that whether this noise removal is part of image enhancement or image restoration.

I have read some papers on it, at some points authors are saying that image enhancement from a given noisy image, and on the other hand in some papers authors are saying that noise removal is image restoration process.

Please explain what is correct?.

It depends on the final goal and the model used.
If we assume we have an image - $ I $ which is a result of a sum on of an image $ W $ with noise $ n $: $I = W + n$ and our goal by reducing the noise is estimating $ W $ than this process is restoration.
Yet I could try achieving a different goal by denoising, such as making it look better (Subjective target) or use prior such as Piecewise Constant image which means the denoising process if more of an enhancement to the image.


Faded line in spectrogram for a chirp I've generated this spectrogram using wavelet convolution.
The chirp is generated like so:
framerate = 10000
DURATION = 10
N = DURATION*framerate
k = np.arange(N)
chirp = sig.chirp(k, 0.01, N, .49)

Each row (linearly representing a frequency) is generated like so:
$$row_f = |IFFT(FFT(chirp) * FFT(wavelet_f))|$$
where $f$ represents the frequency.
When I use matplotlib's built in spectrogram tool, I get this image:

L = 100  # length of one segment
overlap = L/2  # overlap between segments
plt.figure(figsize = (10, 5))
plt.specgram(sounddata, NFFT=L, Fs=2, noverlap=overlap, sides='onesided')

In this plot, the line is seemingly cleaner - not blurring over time. Why is the line blurring? I thought wavelet convolution was supposed to have higher resolution than STFT.

If you think of the spectrograms as filter banks it might be easier to conceptualize. With a STFT, each bin in the bank has the same time/frequency resolution. With the Wavelet transform, the bins at lower frequencies have higher frequency resolution, which means they must also have lower time resolution.  There is an inherent trade off between time and frequency resolution such that one cannot be increased without decreasing the other.  This can be thought of intuitively as an impulse in the time domain has an infinitely wide frequency response, or infinitely high resolution in the time domain and zero resolution in the frequency domain.


Compressive Sensing - Sparse in frequency example I am learning about compressed sensing. I have a question regarding a common MATLAB "sparse in frequency" example that can be find online, for example here and here.

What confuses me in these examples is that they do not take the subset of time samples from the time signal directly, but by inverse transforming the spectrum of the signal. Is this correct? Does this not already involve knowing the frequency spectrum of the signal ? I thought that the goal was to reconstruct the time signal from a subset of time samples, without knowing the exact shape of the spectrum. Or am I missing something?

I think I understand your confusion. In your first example, it would have been more realistic to take

y = x(q)

for example.

A possible explanation for the approach you see here is that authors of such examples would like to ensure they obtain the exact sparse solution they generate the example from to compare against. Therefore they generate the measurements directly from it. It is not "cheating" per se in a mathematical sense, but you can say that it does not reflect the practical reality of taking measurements in the time domain.

What could happen when taking the measurements in the time domain is that if you generate sinusoids with arbitrary frequencies, the signal will appear less sparse in the frequency domain when the frequencies do not match the frequencies of the dictionary. If for example you generate a sinusoid at 1.5 Hz and your DFT has atoms corresponding to the frequencies 1 Hz, 2 Hz, 3 Hz, etc., your signal will not be 1-sparse in this dictionary but will have to use multiple atoms at adjacent frequencies to form the sinusoid at 1.5 Hz. A simple way to avoid this is to generate a signal of the prescribed sparsity directly in the frequency domain and then convert it to the time domain as part of taking the measurements.

Technically, it does not involve already knowing the frequency spectrum of the signal, because X is not known in the reconstruction and the sufficient frequency information in the classic sampling-theoretical sense is lost in the multiplication

y = A*X

because A is a wide matrix.


Calculate autocorrelation of a sinus I am a beginer in signal processing.
I am currently studying a signal $s(t) = \cos(h(t))$, and I have to calculate its autocorellation.
Since $h(t)$ is periodic of period $T_s$, I calculated the value of

$$R_s(\tau) = \frac{1}{T_s} \int_0^{T_s} \cos(h(t)) \cos(h(t-\tau)) dt$$

Honestly, this calculation is not easy (since $h(t)$ is also a cosinus), and in order to simplify my final document, I would like to use complex representation by writing $s(t) = Re\left\{ e^{ih(t)} \right\}$. However, This cannot work for the autocorrelation since $$Re\left\{ e^{ih(t)} e^{-ih(t-\tau)}  \right\} = s(t)s(t-\tau) + Re\left\{ \sin(h(t)) \sin(h(t-\tau))  \right\}$$

So, do you know a way to calculate $R_s$ by calculating $\frac{1}{T_s} \int_0^{T_s} e^{ih(t)} e^{-ih(t-\tau)} dt$ ?

I was thinking about a relation like $R_s(\tau) = \frac{1}{2} Re \left\{\frac{1}{T_s} \int_0^{T_s} e^{ih(t)} e^{-ih(t-\tau)} dt \right\}$, but I cannot prove it :(

It's way easier to use the complex representation of cosine not involving the real part operator. In your case this translates into expressing

$$ \cos(h(t))=\frac{e^{jh(t)} + e^{-jh(t)}}{2} $$

Note that this yields an autocorrealtion of
$$ R_s(\tau)=\int_{0}^{T_s} \frac{e^{jh(t)} + e^{-jh(t)}}{2}\frac{e^{jh(t-\tau)} + e^{-jh(t-\tau)}}{2} dt $$


How to simulate a stereo recording? I have a single channel audio recording. I want to simulate the stereo recording when it is (hypothetically) played on a speaker in front of stereo mics (distance between them - say 10 cm).

Are there any tools which can do this?

Edit- I am looking for something more precise - as in for specific distance between speaker and mics and also for given room dimensions. I tried Audacity 'Reverb tool' which sort of does stereofication but there is no way to specify certain distances (speaker-mic, mic-mic etc.)

Check a music production software (commercial: Logic/ Ableton/ Pro Tools. Open Source: Ardour, Mu Lab) or otherwise audacity. Add a stereo widening tool plugin.
For that check KVR Audio
stereo/ imaging/ widening
stereo imaging more close search


What are the practically relevant differences between various image resampling methods? Mathematica's ImageResize function supports many resampling methods.

Not being familiar with this area, beyond nearest neighbour, bilinear, biquadratic and bicubic (which are obvious from the name), I am lost.

Can you point me to some source that will explain the basic (mathematical) differences between these methods, and in particular point out the practical differences (e.g. by showing sample images where the choice of method really matters and introduces noticeable differences)?

I don't have a signal processing background, so I'd prefer a "gentle" and concise introduction :-)

I'll copy here the list of ImageResize methods for those "lazy" to click the link:

  "Nearest"  nearest neighbor resampling
  "Bilinear" bilinear interpolation
  "Biquadratic"  biquadratic spline interpolation
  "Bicubic"  bicubic spline interpolation
  "Gaussian" Gaussian resampling
  "Lanczos"  Lanczos multivariate interpolation method
  "Cosine"   cosine interpolation
  "Hamming"  raised-cosine Hamming interpolation
  "Hann" raised-cosine Hann interpolation
  "Blackman" three-term generalized raised cosine
  "Bartlett" triangular window interpolation
  "Connes"   squared Welch interpolation
  "Welch"    Welch quadratic interpolation
  "Parzen"   piecewise cubic interpolation
  "Kaiser"   zero-order modified Bessel interpolation

Given an image $I(m,n)$ with $m,n$ integers, the interpolation of that image at any arbitrary point $m',n'$ can be written as

$$\tilde{I}(m',n')=\sum_{m=\left\lfloor m'\right\rfloor-w+1}^{\left\lfloor m'\right\rfloor+w}\ \sum_{n=\left\lfloor n'\right\rfloor-w+1}^{\left\lfloor n'\right\rfloor+w}I(m,n)\ f(m'-m,n'-n)$$

The result $\tilde{I}$ is still only an approximation to the true underlying continuous image $\mathcal{I}(x,y)$ and all that different interpolating functions do is to minimize the approximation error under different constraints and goals.

In signal processing, you'd like the interpolating function $f(m,n)$ to be the ideal low-pass filter. However, its frequency response requires infinite support and is useful only for bandlimited signals. Most images are not bandlimited and in image processing there are other factors to consider (such as how the eye interprets images. What's mathematically optimal might not be visually appealing). The choice of an interpolating function, much like window functions, depends very much on the specific problem at hand. I have not heard of Connes, Welch and Parzen (perhaps they're domain specific), but the others should be the 2-D equivalents of the mathematical functions for a 1-D window given in the Wikipedia link above.

Just as with window functions for temporal signals, it is easy to get a gist of what an image interpolating kernel does by looking at its frequency response. From my answer on window functions:

  The two primary factors that describe a window function are:

  Width of the main lobe (i.e., at what frequency bin is the power half that of the maximum response)
  Attenuation of the side lobes (i.e., how far away down are the side lobes from the mainlobe). This tells you about the spectral
  leakage in the window.

This pretty much holds true for interpolation kernels. The choice is basically a trade-off between frequency filtering (attenuation of sidelobes), spatial localization (width of mainlobe) and reducing other effects such as ringing (Gibbs effect), aliasing, blurring, etc. For example, a kernel with oscillations such as the sinc kernel and the Lanczos4 kernel will introduce "ringing" in the image, whereas a Gaussian resampling will not introduce ringing.

Here's a simplified example in Mathematica that let's you see the effects of different interpolating functions:

true = ExampleData[{"TestImage", "Lena"}];
resampling = {"Nearest", "Bilinear", "Biquadratic", "Bicubic",
   "Gaussian", "Lanczos", "Cosine", "Hamming", "Hann", "Blackman",
   "Bartlett", "Connes", "Welch", "Parzen", "Kaiser"};
small = ImageResize[true, Scaled[1/4]];

Here, true represents the image which I assume to be the discrete equivalent of the "exact" image $\mathcal{I}(x,y)$, and small represents a smaller scale image $I(m,n)$ (we don't know how it was obtained). We'll interpolate $I(m,n)$ by 4x to give $\tilde{I}(m',n')$ which is the same size as the original. Below, I show the results of this interpolation and a comparison with the true image:

You can see for yourself that different interpolating functions have different effects. Nearest and a few others have very coarse features and you can essentially see jagged lines (see full sized image, not the grid display). Bicubic, biquadratic and Parzen overcome this but introduce a lot of blurring. Of all the kernels, Lanczos seems (visually) to be the most appealing and one that does the best job of the lot.

I'll try to expand upon this answer and provide more intuitive examples demonstrating the differences when I have the time. You might want to read this pretty easy and informative article that I found on the web (PDF warning).


Cochleagram vs STFT for CNN-based speech segregation I’m working on a project (python-based) that would use ideal ratio masks (IRMs) as a basis for cleaning noisy speech in various environments. Specifically, this will be accomplished through the use of a convolution neural network, where input features would be either a cochleagram or STFT data matched with its respective IRM label.

(This is the article I’m basing my project on --> http://web.cse.ohio-state.edu/~wang.77/papers/CWYWH.jasa16.pdf)

As I’m (relatively) fairly new to learning about digital signal processing, I was hoping someone would be able to give me an explanation hilighting some of the differences between cochleagrams and STFTs, and what their individual use-cases would be.

Many thanks in advance!

I was hoping someone would be able to give me an explanation
  hilighting some of the differences between cochleagrams and STFTs, and
  what their individual use-cases would be

STFT (short time fourier transform) is a technique that is used to obtain a time-frequency representation of audio signals. Spectrograms and phase vocoders are arguably the most important applications of the STFT.

The cochleagram, (which is really just a variant of the spectrogram), is also used to obtain a time-frequency representation of audio signals but has one major difference compared to STFTs (and spectrograms).

Regular spectrograms have a constant channel (frequency bin) bandwidth, which is given simply by $f_s / N$ (the sampling rate divided by the number of samples ). In contrast to that, cochleagrams' channels have different bandwidths, which are chosen to mimic human hearing  (the cochlea).
There're more than just one (empirical) equation to compute the cochleagram bands and one of them can be found here  ( https://en.wikipedia.org/wiki/Critical_band )

$$ERB(f) = 24.7 * (4.37 f / 1000 + 1) (eq 1)$$

In a cochleagram, the bandwidths grow increasingly larger with increasing frequency. Note also, that it's possible to create a cochleagram by using a normal STFT and then averaging the energy found in several frequency bins so that they match the bandwidth(s) given by eq #1. Alternatively, you can use gammatone filters, which are basically band-pass filters of different bandwidth.


Estimating the Nyquist rate when simulating data I want to find the power spectrum $\tilde{Q}(\omega)$ of a quantity $Q(t)$. The values of $Q$ are computed at discrete time steps so I have a series $Q(\Delta t)$, $Q(2 \Delta t)$, ... of values, but the calculation at each time step is independent rather than (for example) needing to know the value of $Q(\Delta t)$ in order to compute $Q(2 \Delta t)$. I then compute the power spectrum of $\tilde{Q}(\omega)$ by an FFT of the signal.

My main question is what is a good way to choose $\Delta t$ i.e. estimate the Nyquist rate of $Q(t)$ without wasting a lot of computing time?

I'm also wondering, given that the calculations at different times are independent, is there a smarter or more efficient way to find $\tilde{Q}(\omega)$ than generating a time-series and performing an FFT.

You must know the method of computing Q(t), and if Q(t) = Q[f(t)], meaning that Q is a function of another function, like the output of a filter. Then, knowing the autocorrelation or power spectrum of f(t), one can use various methods to find the power spectrum of Q. If the function Q is linear with respect to f(t), then the solutions are simple.

However, it seems like you should clear up your needs. An FFT of the signal is simply its spectrum, not is power spectrum or its power spectral density. They are slightly different.

To answer your last question, it seems like you are formulating this in a probability/statistics sense when you say independent? If they are truly independent (i.e. one value of Q gives absolutely no clue about the next value), then your power spectral density represents white noise and it is flat (a constant). Again, it would be best if you clarified your definition of "independent". For example, if Q is a relatively smooth function of time, then it is not statistically "independent" as you say, and to find the spectrum or power spectrum density, you should analytically look at other functions that Q depends on first.


2D DWT computation order In 2D the discrete wavelet transform (DWT) of an image using lifting based 5/3 filter, if I perform a row-wise operation first then perform column-wise operation then I will get 4 sub-band LL, LH and HL and HH.

During reconstruction, that is IDWT, should I maintain the order in reverse direction that is column first and then row?
Or any order I can maintain?
In JPEG 2000 there is no order information found in the header. Then how I can main the order during reconstruction?

For each single level of a isotropic or square DWT, the row/column order is not important.  Note however that there exists an alternative decomposition, termed anisotropic, rectangular, that does all rows first, then all columns (or the converse) which results in a different subband organization. This is described in the overview of 2D wavelets: A panorama on Multiscale Geometric Representations (arxiv) in Section 3: 2.3.2. Separable orthogonal wavelets.

[EDIT: Added references to related questions]

2D DWT Image Issue
What is the correct order of operations for a 2D Haar wavelet decomposition?
2D DWT computation order