Cutoff frequency of a first order recursive filter I've implemented a very simple first order recursive low pass filter in c using the algorithm:

double a, b, prevOutput;

double lowPassFilter(double input)
{
    double output = a * input + b * prevOutput;
    prevOutput = output;
    return output;
}

where:

a = 1.0 - b

According to The Scientist and Engineer's Guide to
Digital Signal Processing, the cutoff frequency of the filter is defined by the relation:

  $$b = e^{-2πf_c}$$

$f_c$ is the -3db cutoff frequency as a fraction of the sampling frequency. Using $f_c = 0.25$, I calculated that $b = 0.20788$ (approximately) and plugged it into my filter. I wanted to verify my filter, so I passed an impulse function through it and recorded the output:

double result[128];
for (int i = 0; i < 128; ++i)
{
    if (0 == i) result[i] = lowPassFilter(1.0);
    else result[i] = lowPassFilter(0.0);
}

I then performed FFT on the result to find the frequency response:

complex_t fftResult* = fft(result);
for (int i = 0; i < 65; ++i)
{
    double amplitude = fftResult[i].real * fftResult[i].real;
    amplitude += fftResult[i].img * fftResult[i].img;
    amplitude = sqrt(magnitude);
    printf("f = %.3f: %f\n", i / 128.0, amplitude );
}

The result I saw was somewhat expected but also not exactly what I expected. In the DC bin ($f = 0$) I saw a amplitude very close to one and I saw the amplitude attenuate as it approached 0.5. However, I also expected that at $f = 0.25$ the amplitude should have been around $0.70795$ (i.e. -3db). but instead the amplitude was around $0.775$:

 ...
f = 0.242: 0.782501
f = 0.250: 0.774971
f = 0.258: 0.767380
 ...

Aside from the inefficient code, is there something wrong with the method I am using to obtain the frequency response or am I interpreting the cut-off frequency relation incorrectly or something else?

That formula for the cut-off frequency is a very inaccurate approximation. In this answer I derived the exact relation between the coefficient of a first order recursive averaging filter and its 3-dB cut-off frequency. Note that in the quoted answer I used the constant $\alpha=1-b$. From formula $(3)$ in that answer we get for the coefficient $b$

$$b = 2-\cos(\omega_c)-\sqrt{(2-\cos(\omega_c))^2-1)}\tag{1}$$

With $\omega_c=2\pi f_c$ we get for $f_c=0.25$ a value $b=0.26795$.

The frequency response of the filter is given by

$$H(e^{j\omega})=\frac{1-b}{1-be^{-j\omega}}\tag{2}$$

The figure below shows the magnitude response $|H(e^{j\omega})|$, and it shows that the desired cut-off frequency is achieved with the chosen value of $b$:

Note that the 'EDIT' part of the answer quoted above also explains where the inaccurate formula for the cut-off frequency comes from. The following figure shows the resulting cut-off frequency of the filter as a function of the desired cut-off frequency when the exact formula $(1)$ is used and when the approximation is used:

Clearly, the approximation is quite bad and it only kind of works for very small cut-off frequencies. I see no reason why it should be used at all.


Intuition behind the Gaussian Filter in Image Processing I've been using a Gaussian Filter and know how to use it but I don't understand the physical significance behind using the filter. We use it but why? What is the significance? I'm just starting out in image processing, so please no high level responses. Thanks!

As far as your comment suggests, you are using a LP Gaussian Filter.
LP = Low pass.
This means it passes the low frequency component of the image, and prevents the high frequency component.

In image processing, you can imagine the high frequency component as the places in the image where you have a sudden jump in pixel values. Like moving from a black background and suddenly to white one.

These components are associated with edge detection normally, but at the same type they are associated with noise. It's because noise has a high frequency.
How come?
There is for example a type of noise known as "Salt and Pepper" where there appears spikes in the image. You see the it as a weird dot in the image.
These normally have way different value than the surrounding.

With this sort of a filter (Low Pass Gaussian Filter) you can suppress the noise, but you will lose some sharpness in the image due to killing some parts of the edges.
Therefore; this is called smoothing. It's because you tend to make the image smooth in a way that pixel values are not changing rapidly but rather in an elegant way.

You can imagine this as averaging. It's because you are substituting the values in a certain region of the image by it's average. That will help in smoothing the image.


Is there any alternative basis for a fourier-like transform? Fourier transform of a continuous signal is just the projection of the signal on the sinusoidal family for imaginary part and the same family with phase offseted by a quarter of period for the real part.
Beside the magical Euler's formula $e^{ix} = \cos x + i\sin x$ that ease the analysis of the transform and the fact that sinusoidal functions happen to be the solution of a prefect harmonic oscillator, I always wondered where the Fourier supremacy came from. Sinusoidal signals aren't so common in everyday's life. Even with human speech, vocal tract function isn't sinusoisdal.
For example, one transform basis that would seams interesting to me would that with a saw tooth function as signals with fast attack and slow release aren't uncommon...
I guess that the transform need to be inversible and that might not be the case for every basis, but aren't there any alternative basis that provide a bijective forward and inverse transform ? (And if so, what would be the conditions ?)

I will elaborate a little. I'm looking for a basis family of periodicals functions named $B(t)$ for example with B of unitary periods such as : $$ S(f) = T\{s(t)\}(f) = \int s(t) B\left( t * f\right) dt + i \int s(t) B\left((t - \frac{1}{4})*f\right)dt$$ with T inversible in some way...

The most important reason for the "Fourier supremacy", as you called it, has been pointed out by Jazzmaniac in his comment. The complex exponential $e^{st}$ with $s=\sigma+j\omega$ is an eigenfunction of the convolution operator:

$$e^{st}*f(t)=\int_{-\infty}^{\infty}e^{s(t-\tau)}f(\tau)d\tau=
e^{st}\int_{-\infty}^{\infty}e^{-s\tau}f(\tau)d\tau=e^{st}F(s)$$

where the eigenvalue $F(s)$, if it exists, is the Laplace transform of $f(t)$, and for $\sigma=0$, i.e. for $s=j\omega$, $F(j\omega)$ is the Fourier transform of $f(t)$ (again, if it exists). Even though there are several important differences between the Fourier and the Laplace transform, they are deeply related, and for answering your question their differences are not important. By analogy, the same is true for the discrete-time Fourier transform and for the $\mathcal{Z}$-transform in the discrete domain. So when in the remainder of my answer I mention the Fourier transform, you might as well replace it by Laplace transform or $\mathcal{Z}$-transform, and everything will remain valid.

So why should we be concerned with eigenfunctions of the convolution operator? Because convolution describes the input-output relation of linear time-invariant (LTI) systems, and these systems are often very good approximations of practically important systems such as frequency selective filters, integrators, differentiators, etc.

Apart from being an invaluable tool for the analysis of LTI systems, the Fourier transform also has an important physical interpretation as the spectrum, i.e. the frequency content, of a signal. Furthermore, the FFT is a highly efficient algorithm for the computation of its discrete version for finite length signals.

So why don't we just say we're doing fine with the Fourier transform and its siblings? Because there are applications where we are looking for other properties of a transform. One example would be time localization in signal analysis, since the basis functions of the Fourier transform have no time localization at all because of their infinite extent. In such cases other transforms have become widely used, such as the short-time Fourier transform and the wavelet transform. For discrete-time signals, filter banks provide very general signal expansions. They can be designed to match the desired properties of the signal expansion.

Since you mentioned the saw-tooth function as a possible basis function I would like to point out to you the Slant transform (mentioned in this answer), which is used in image coding. Furthermore, there is the related Slantlet transform, which is a special type of discrete wavelet transform with piecewise linear basis functions. You can see the basis function in Fig. 5 of this paper.


Analogue filter analysis (band-pass) The transfer-function:
  $$G(s) = \frac{\beta s}{s^2 + \beta s + \omega_0^2}$$

  is to be used in an application that requires the magnitude of the frequency response to be of the band-pass form.

  Show that this property is present by considering the cases when $\omega = 0, \omega = 1$ and show also that $\lvert G(j\omega_0)\rvert = 1$.
  Determine also the equation that defines the frequencies where the filter magnitude is $−3\textrm{ dB}$.

I use $G(s)=G(j\omega)|_{s=j\omega}$, then
$$
\lvert G(j\omega_0)\rvert= \frac{j\beta \omega_0}{-\omega_0^2 + j\beta\omega_0 + \omega_0^2}=1
$$
when $\omega=0$, the numerator equal to $0$ so $\omega=0$ is be filtered out.

My questions are:

Could I use the $G(s)=\lvert G(j\omega)\rvert$ directly or $G(s)=G(j\omega)=\lvert G(j\omega)\rvert e^{-jk\omega}$?
I have no idea about the case when $\omega=+\infty$?
The last question, I use $-20\log\lvert G(j\omega)\rvert=-3\textrm{ dB}$, so we can solve the $\omega$ from this equation, is that correct?

Even though this may not be homework, it is a basic textbook problem, and for this reason I just give you a few hints, which should nevertheless help you solve the problem if you put in a bit of effort.

All questions refer to the filter's frequency response, which is simply obtained by evaluating $G(s)$ for $s=j\omega$.
The first question should read "Show that this property ... $\omega=0$, $\omega={\infty}$" (not $\omega=1$).
Your expression for $|G(j\omega)|$ is wrong; you wrote down $G(j\omega_0)$ (without magnitude). Of course, if $G(j\omega)=0$, or $G(j\omega)=1$, then also $|G(j\omega)|=0$ and $|G(j\omega)|=1$, respectively.
In order to see what happens for $\omega\rightarrow\infty$, write down $|G(j\omega)|^2$ (that's maybe easier than the magnitude), and take the limit.
For the 3 dB cut-off frequencies (there must be two positive solutions), use $|G(j\omega)|^2$ derived above and solve $|G(j\omega_c)|^2=\frac12$. Substitute $x=\omega^2$ and solve the resulting quadratic equation for $x$. From that solution, derive the two positive solutions for $\omega_c$.

All of this involves basic math skills, such as computing the (squared) magnitude of a complex number, and solving quadratic equations. If you find that difficult, read up on these basics because you'll need them all of the time when doing signal processing.


What are the differences between estimation and detection optimal filters (matched and Wiener)? What is the differences between the two of them? I know that the matched filter is a detection filter because it detect what is the peak of the signal and return the time delay of reflected and transmitted signal. Whereas the Wiener filter is an estimation filter because it estimates the mean-square-error estimation of the degraded signal. But I'm not too sure on what is the definition of the detection filter itself, and the estimation filter definition. Any help would be appreciated.

If you want to maximize the signal-to-noise ratio using a linear filter in the presence of additive random noise, then it can be proved that matched filter is the optimal filter. It is usually performed for detection of a known waveform.
If you want to minimize the mean square error of estimation of a desired signal corrupted by additive random noise by a linear filter, then it can be proved that wiener filter is the optimal filter. It is usually deployed to estimate a known stationary target signal.

About what detection and estimation is, consider the radar example. Detection is to find out if there is a plane/target. If it turned out there was actually a target and it was detected, then we would like to estimate its speed, altitude, etc.

A more relevant example for detection is digital communication. Assume we use on-off keying, that is, we send a waveform if we want to represent a bit $1$, and don't send anything if the bit was $0$. During the time interval of communication of each bit (timeslot), a liner filter which is matched to the waveform of bit $1$ can be used to detect if there is a bit $1$. If no waveform was detected, then it is assumed that $0$ is sent.

An example for wiener filtering is signal compression. Assume you want to compress a stationary process $x(t)$. One approach is to estimate it using a wiener filter whose input is just white noise. Given the target signal $x(t)$, the wiener filter can form the spectrum of the white noise such that its output has the minimum mean square error compared to $x(t)$. Now the coefficients of the wiener filter (plus a source of white noise) would be enough to reproduce $x(t)$.


Average bit error and symbol energy How do you draw decision boundaries given this type of 16-qam constellation diagram, I know that there should be no vacant decision boundary but I can't seem to get it.Also, assuming that the correct decision boundary is drawn, how do you derive its total probability of error, if we let the SNR ratios of all symbols are equal?

The optimal decision regions are the Voronoi Regions. I dont know, if this is what you are looking after.

import numpy as np

points = np.array([(1,1), (1,-1), (-1,1), (-1,-1), (3,3), (3,0), (3,-3), (0,-3), (-3,-3), (-3,0), (-3,3), (0,3), (5,0), (0,5), (-5,0), (0,-5)])
from scipy.spatial import Voronoi, voronoi_plot_2d
vor = Voronoi(points)
voronoi_plot_2d(vor)
plt.axis('equal'); plt.xlim((-7,7)); plt.ylim((-7,7));

I actually, consider analytically calculating the error probability of this constellation a very hard task, since the regions are not very regular.


Integrator and the multiplier in the Optimal Binary Detector The optimal binary detector for asymmetric, equiprobable PAM pulse normally have structure as shown below,

But I cannot understand what is the multiplier and the integrator using for.

Why the detector need to care the signal difference and multiply with the received signal?
Why can't I directly put the received signal to the threshold comparsion?

There are several things missing/extra in your diagram.

What you are using is
rectangular PAM pulses of duration $T$ to send data across the channel, and so
you really don't need the multiplier. It is necessary only if $s_1(t)$ and $s_2(t)$
are different from rectangular pulses (though they are still of duration $T$,
and in that case, the input $s_1(t)-s_2(t)$ to the multiplier does not tell
the whole story.
The integrator that you are using is really what is called an
integrate-and-dump correlator. The input is
integrated over intervals of time $[(n-1)T, nT)$. Just before the end of
the integration period, at time $nT^{--}$, the integrator output is sampled
and the sample value is what is used in the threshold device. Following
the sampling at $nT^{--}$, the integrator is dumped (meaning that
the output is reset to $0$) at time $nT^{-}$, and the integrator is
restarted at $t = nT$ for the next data bit from $nT$ to $(n+1)T$.
Note that there is what computer engineers call a critical race involved:
the circuit designer has to make absolutely sure that the sample signal
arrives at the integrator before the dump signal. Getting this wrong
can be a fireable offense.
With the integrate-and-dump correlator working as described above,
the input to the multiplier during $[(n-1)T, nT)$ needs to be
$s_1(t-(n-1)T)-s_2(t-(n-1)T)$, that is, the signal $s_1(t)-s_2(t)$
(which is nonzero only for $t \in [0,T)$) delayed by $(n-1)T$ seconds.
What the integrate-and-dump correlator is doing is computing the
output of the matched filter for signals $s_1(t)$ and $s_2(t)$
in additive white Gaussian noise at times $nT$.  Note that the
correlator output is not equal to the matched filter output
at times other than $nT$.
The reason you cannot simply sample the received signal
and make a decision from the sample is that the SNR is abysmal.
Technically, AWGN has infinite variance
(cf. this question),
but even otherwise this is a big problem.  The matched
filter maximizes the SNR at the sampling instants,
and the correlator is just calculating what the matched
filter will give you at the sampling instants.

For details of all this, see, for example, pp. 85-93
of these ancient lecture notes of mine.


Shift invariant in wavelet I always hear that wavelet transform is not shift invariant, and that there are other types of wavelet, like stationary wavelet and double density dual-tree wavelet transform, that are shift invariant.

Can anyone explain to me, what is the meaning of "shift invariant" .

Let's say you have a signal which is all zeros except for a spike at one point where x(8)=1 (total N=32, for example).  If you perform the DWT on this signal and then calculate the total energy (by taking the square root of the sum of the squares of all the results), you will get a value - call it "E1".

Now, let's take another signal which is still all zeros except at one point, n=12, where x(12)=1 - JUST LIKE THE FIRST SIGNAL.  It is the same signal except it is shifted by 4 points.  If you take the DWT of this signal and calculate its total energy (again by taking the square root of the sum of the squares of the results), you get another value - call it "E2".

You might expect that the total energy of these two would be the same - but they're not!  E1 will NOT equal E2, so therefore the DWT is NOT "shift-invariant" because the energy "varies" whenever you "shift" the incoming signal - even though it's basically the same signal.

So transforms like the dual-tree complex discrete wavelet transform have the property that they are "shift-invariant", which means that you can shift the signal around before calculating the transform and the resulting energy will still be the same.  This is important for all sorts of reasons, but if you think about it, shift-invariance is a really good thing to have because it means that it doesn't matter WHERE in the signal you start calculating - it will still have the same energy coming out.

This page has a pretty good description and demonstrations non-shift-invariant transforms in 1D, 2D, and 3D.

Here's a really simplistic example that (I hope) illustrates the main idea of shift-invariance:

Suppose you are trying to find the wavelet transform of a piece of music you have, but the WAV file you used has 5 seconds of silence before the music starts.  You get the wavelet transform coefficients and graph them nicely on a scalogram (a time/scale graph - which you can think of as a time/frequency graph).  You notice that it indicates 5 seconds of silence (0 values) before the musci starts.

Now you edit the WAV file to cut out the 5 seconds of silence down to 1 second of silence (you are shifting the sample in time by doing this - hence "shift" in shift-invariance).  You take the discrete wavelet transform again and plot it in a scalogram and you expect to see the same thing except without the silence - you expect it to just be shifted over to the left, but...

It won't be the same!  In fact, if you take the total energy of the resultant signal, it won't even be the same, because the DWT by itself is NOT shift-invariant.

But if you use the UDWT (also called "stationary" DWT), then the two images will be the same except the second one will not have all that silence in the beginning.  You could take a screenshot of the first one and crop it in photoshop to get the same image as in the second one.

The dual-tree complex DWT (DTCDWT) is near-shift-invariant, meaning it's very close to being shift-invariant.

The term literally means "does not vary when shifted".

Hope that helps!


Finding Reference Audio Signal in Test Audio Signal and Cropping Accordingly As seen in the diagram, below I have a reference audio and a test audio. I want to find at what part of the test clip the reference audio be heard. Once, that is found I want to crop the test file from the point where the match begins and the point where the match ends. How can I do this matching and cropping of audio signals using Python? Some code will be of help.

The diagram below(from wiki) is for illustration purposes only

My audio clips:
Reference audio:  africaChirp
Test audio:       Africa

I want to find the time in seconds at which africaChirp appears in Africa. africaChirp was cropped from the source file, Africa, found here.

I have also included my code to show my progress:

import numpy as np
#import librosa
import logging
from scipy.io import wavfile

rate2,test = wavfile.read('africa.wav')

#indices
testStart  = 6000;
testEnd    = 200000;

refStart   = testStart + 11000;
refEnd     = refStart + 2000 - 1;

testSig=test[testStart:testEnd]
refSig=test[refStart:refEnd]

refStartRel=refStart-testStart+1
refEndRel=refEnd-testStart+1

no_samples_test=testSig.shape[0]
no_samples_ref=refSig.shape[0]

xcorr=np.zeros((no_samples_test-no_samples_ref+1,1))
xcorr_norr=np.zeros((no_samples_test-no_samples_ref+1,1))

refSig_norm= np.linalg.norm(refSig)
#print(xcorr.shape[1])
for i in range(1,xcorr.shape[0]):
    testSig_samples=testSig[i:(i+no_samples_ref)]
    xcorr[i]=np.sum(np.multiply(testSig_samples,refSig))
#    norrFilt=(filt-np.mean(filt))/(np.std(filt))
    linalgnorm=np.linalg.norm(testSig_samples)
    xcorr_norr[i]=xcorr[i]/np.multiply(refSig_norm,linalgnorm)

xcorr_max_id=np.max(np.abs(xcorr_norr))

print(xcorr_max_id)

If you have a reference signal you want to find in a different signal then your model matches almost perfectly (Up to the environment the signal to be found is in) to Matched Filter.

So basically you need to do cross correlation between the Test Signal and the Reference Signal.
Find the point of maximum correlation and create a cropping zone around it according to the length of the reference signal.

Update

I downloaded the files and wrote a MATLAB Code to do the task of finding a reference signal within a signal.
I generated an equivalent test case by cropping 14,000 samples from the Song File - Toto - Africa as the Test Signal and from them cropped 2,000 samples as the Reference Signal.

Matched Filter in the case above must be tweaked to normalize different volume levels in the signal.
Hence the correct way to so is Matched Filter which is normalized by the Norm of the sections being cross correlated (This is actually correlation in Statistics).

The the result it as following:

As can be seen, the algorithm detects the exact starting index of the reference signal.
Pay attention that the classic cross correlation yields the wrong answer as it is being tricked by higher volume levels of the test signal.
When we take care of different volume level by the normalization then the maximum correlation happens exactly where it should.

The MATLAB code is available at my Signal Processing StackExchange Question 50003 - GitHub Repository.


How to Determine if FIR has a Linear_Phase Response w/o Matlab Given the FIR transfer function:

h(z) = $ .36 + .384z^{-1}  + .1608z^{-2} +.9712z^{-3}  + .352z^{-4} + .18z^{-5} - .2z^{-6} $

How do you determine if this transfer function has a linear - phase response.

Real-valued and causal FIR filters have a (generalized) linear phase response if and only if their coefficients satisfy either

$$h[n]=h[N-1-n]\tag{1}$$

or

$$h[n]=-h[N-1-n]\tag{2}$$

where $N$ is the filter length (number of taps).

There are four types of linear phase FIR filters, depending on the type of symmetry (even as in $(1)$ or odd as in $(2)$), and the filter length $N$ (even or odd). They have different properties, which is explained in this more detailed answer.


matching the mouth shape of an audio file(lip sync) I am working on something like cartoon, but the character's mouth shape is automatically detected and matched by the word they are saying.

my guess is to use speech recognize algorithm to analyze the wav file, and get the content(word) at the exactly moment, for the exactly word, prepare corresponding mouth shape picture.

but it look like a dumb way... is there any convenient way to work this out?

I saw many games support this, how did they achieve that?

The simplest solution I can think of is to use an envelope follower, the higher the amplitude, the more open the mouth is. This is pretty simple but will one give tell you how 'open' the mouth should be, not the shape. This should already be quite effective for a cartoon. To make it more realistic you should probably try to detect vowels (Oh and Ah for instance) by measuring the formants. http://en.wikipedia.org/wiki/Formant gives a list of average frequency's for some vowels. Of course these are only average, so you need to experiment.

I don't think you need to recognize exactly which words are spoken, this will be much more complex, and there are quite a lot of words in each language, and quite a lot of languages... :-)

If you add consonant detection as well, you probably have all the information you need...


Music note of different instruments What do I need to know to difference a note for different musical instruments? for example, if they play the same note, (same frequency), how can I difference them if they have the same frequency?

Any (melody) musical instrument note has a fundamental frequency (think "a sine wave") which our ears tend to interpret as a "pitch", and a number of harmonics (extra sine waves) which are integer multiples of that fundamental frequency. These extra sine waves with double, triple, ... times the fundamental frequency add up to something that our ears tend to interpret as a timbre. The reason for why you find this special relation between fundamental frequency and harmonics is a consequence of the physics of how real instruments make sound.

The difference between an "A4" note on a piano and an "A4" note on a violin is in which harmonics are present in the sound, and how strongly they are present. One way to examine this is by comparing their respective spectrograms. Several audio tools allow you to draw a spectrogram for a given sound (e.g. sonic visualizer or audacity).

For percussive instruments there's no easy relation between fundamental frequency and harmonics and the spectrum typically looks much noisier (many more high frequency components). Also the spectrum at the start of a melody instrument note (the so-called "attack") will typically look noisier.


Removing low frequencies from a signal I am trying to remove low frequencies from a signal and intuitively I chose the high-pass filter, more specifically - a Butterworth filter, Order 4 (because I am not sure how to choose properly the order and 4 seemed as a good choice) and cutoff frequency of 50 Hz. The problem is, that the filter removes the low frequencies, however, the peaks that were caused from those low frequencies are still to be seen.

I have read also about adaptive filtering, but according to the literature, I need a reference signal, which would be used as a desired signal. I tried implementing an adaptive RLS/LMS filter and as a reference signal, I processed my data with a low-pass filter - Butterworth, 4th Order, 50 Hz for the cutoff frequency. This method also did not get me far.

I have provided a a copy of my data, which is sampled at 1500 Hz.

You want to remove the heart beat signal and keep the "noise". We can solve this problem by using a denoising algorithm, and subtracting the denoised signal from the original signal.

Setting frequency cutoffs for a frequency domain filter can get tricky and turn into a game of whack-a-mole because there's "high frequency" components in the heartbeat blip (due to the sudden rise and fall) and also in the wiggly stuff between a heart beat's T wave and the next heart beat's P wave.

Loosely speaking, the requirements in this de-noising problem are as follows:

Remove the little wiggles
Maintain larger jumps that appear in a heart beat PQRST waveform

This sounds like a great place to apply $\ell_1$ denoising or total variation denoising. The idea is to approximate the given signal $y$ with a signal $x$ such that the derivative of $x$ is "sparse" i.e. it doesn't change too frequently, but when it changes, the change is large. The denoised estimate $\hat x$ is given by:
$$
\hat x = \min_x ||y-x||_2^2 +\lambda\sum|x_i-x_{i-1}|
$$

I used proxTV toolbox in Python to solve this optimization problem.

import scipy.io as sio
import numpy as np
import matplotlib.pyplot as plt
import prox_tv as ptv

mat_struct = sio.loadmat('Signal1.mat')
noisy_signal = mat_struct['x'].T[0]

filtered_signal = ptv.tv1_1d(noisy_signal, 50)

time_vec = np.linspace(0, len(noisy_signal)/1500., len(noisy_signal))

plt.close('all')

fig, ax = plt.subplots(3,1,sharex=True)

ax[0].plot(time_vec,noisy_signal)
ax[0].set_title('noisy signal')

ax[1].plot(time_vec,filtered_signal)
ax[1].set_title('filtered signal')

ax[2].plot(time_vec,noisy_signal - filtered_signal)
ax[2].set_title('noise')
ax[2].set_xlabel('time (s)')

plt.tight_layout()
plt.show(block=False)

And here's the resulting plot:

Of course, there's a different kind of whack-a-mole you'll have to play with this technique: the $\lambda$ parameter which I set to $50$ in the code.


How to understand the discrete and continuous convolution in matlab? For example,I want to calculate the convolution of two continuous signal in matlab as the image shows.
Indeed, I chose to discrete the signal to compute in matlab.

Obiviously, the 1st and 3rd are right, I think the difference of discrete and
continous convolution might be the reason, but why?

As Hilmar said you need to be careful of the difference between circular convolution and linear convolution. Linear convolution in the time domain translates to multiplying discrete time fourier transforms (DTFT) or continuous time fourier transforms (CTFT) in frequency and I'm sure you've seen this done.

But as you know, the DTFT and CTFT are infinite sums/integrals so cannot actually be computed by a computer. The discrete fourier transform (DFT) can be computed by MATLAB though and it uses an efficient algorithm called the fast fourier transform (FFT) which runs on the order of NlogN time, where N is the length of the signal. You can look up that multiplying two DFTs is equivalent to performing circular convolution which just "wraps" the signal around (read Oppenheim's DSP book or wikipedia).

Check this out which shows the difference between performing linear convolution, how to use the fft() function in MATLAB to get the result, and how circular convolution is different.

clear,close all;

a = [zeros(40,1); ones(20,1); zeros(41,1)];
b = [zeros(50,1); ones(30,1); zeros(21,1)]/2;
c = conv(a,b);
nfft = length(a)+length(b)-1;
A = fft(a,nfft);
B = fft(b,nfft);
C = A.*B;
d = ifft(C);
e = ifft(fft(a).*fft(b));

figure
subplot(6,1,1)
plot(a)
title('a[n]')
subplot(6,1,2)
plot(b)
title('b[n]')
subplot(6,1,3)
plot(c)
title('Convolution in time domain')
subplot(6,1,4)
plot(d)
title('Convolution by multiplying in freq. domain')
subplot(6,1,5)
plot(e)
title('Circular convolution by multiplying in freq. domain')


Convolution v.s. Cross Correlation in Digital Image Processing I understand that from mathematical point of view, only difference between Convolution and Cross Correlation is that Convolution is commutative, while Cross Correlation is not.

In many articles Cross Correlation is explained as measuring similarity between two signals, and Convolution is explained as calculating relationship between two signals.

My question is this, what is the point of commutative property? How does it make Convolution different from Cross Correlation in practice? Or is it just there for "better" mathematics and easier proofs?

Apart from the formulae, let us go back to their actual meaning, and how they are derived. Talking about convolution: this operation is inherent to Linear-Time-Invariant (LTI) systems. In other words: if

you want to analyse a system that is linear, and time-invariant,
or you want to apply a processing or a filter that does not vary in time or space

then the convolution is the ONLY suitable operation.

While cross-correlation is a measure of similarity between two series, computed as a function of the displacement of one relative to the other. Many choices are plausible, the most common is "bilinear" and compatible with energy preservation.


What is the transfer function $H_{xy}$ of this flow graph? I am trying to get transfer function $H_{xy}$ from this graph. I managed to get the other one but I am struggling with this one. Any help or guidance?

Thank you

If you mark with $v[n]$ the signal between summing nodes for input $x[n]$ and output $y[n]$, then:

$v[n] = x[n] -8v[n-2] -4v[n-3]$

Re-arranging sides:

$v[n]+8v[n-2]+4v[n-3]=x[n] $

Using the Z-transform:

$V(z)[1+8z^{-2}+4z^{-3}]=X(z)$

Now let's calculate the transfer function between $v[n]$ and $x[n]$:

  $H_{x\to v}(z)=\dfrac{V(z)}{X(z)}=\dfrac{1}{1+8z^{-2}+4z^{-3}}$

Now moving to the relation between output $y[n]$ and $v[n]$:

$y[n]=v[n]+2v[n-1]$

Again calculating the Z-transform:

$Y(z)=V(z)[1+1+2z^{-1}]$

Now converting the transfer function $H_{vy}(z)$:

  $H_{v\to y}(z)=1+2z^{-1}$

Final step is to merge these two together. We know that in time domain it is a convolution: $h_{x\to y}[n]=h_{x\to v}[n]\star h_{v\to y}[n]$, thus in Z-transform domain it will be the multiplication, yielding the final result:

  $H_{x\to y}(z)=H_{x\to v}(z)H_{v\to y}(z)=\dfrac{1+2z^{-1}}{1+8z^{-2}+4z^{-3}}$


Why does Matlab spectrogram look stripy when windowed input is zero-padded? I am working on a little project, and stumpled upon something I cannot understand yet. I am learning cognitive science, and don't have DSP knowledge background, only learning from what I find on the internet. I have got an 46380x1 double array, which contains my signal. I would like to draw a nice spectrogram out of it like this:
%%signal = %%too long to put here, 46380x1 double array
winsize = ceil(length(signal)/10); %% =4639
win = rectwin(winsize);
shift = ceil(length(signal)/150); %% =310
nfft = winsize;
srate = 78;

spectrogram(signal, win, winsize-shift, nfft, srate, 'yaxis');
ylim([0.015 2]);

It works fine, and shows this:

But if I would like to zero pad the windowed signal (Matlab docs say, increasing nfft will do that), the output gets these stripes in it, and no matter what factor I scale "winsize" by, It will always show some kind of distortion in the spectrogram. For example I set nfft=round(1.5*winsize) before making the spectrogram, and the its output changed to this:

Why is this happening? I know I cannot increase the amount of information in the signal by zero padding, but it might be more readable to the human eye if not so pixelated, that's why I am trying to apply zero padding here. Can you please help, what I should try to achieve this without getting those lines?
Anyway, there is a huge magnitude output component of the spectrogram, throughout the whole time axis, at the very close-to-zero frequency, if I view the trajectory starting from 0hz. view(-135,65); ylim([0 2]); Why is that?

Thanks in advance!

Since I don't have your original data, I'm going to give the answer from my intuition. Try Hann window function and 50% overlap and see what happens.
I have tried a human voice signal as input and whether zero-padding is applied doesn't change the results.
% signal = ...
winsize = 512;
win = hann(winsize, 'periodic'); % use hann window instead of rectangular window
shift = winsize/2; % 50% overlap
nfft = winsize;
% nfft=round(1.5*winsize);
srate = 48e3;
figure;
spectrogram(signal, win, winsize-shift, nfft, srate, 'yaxis');

Here are the results:

Forgive me for the Chinese xlabel and ylabel, they are respectively "Time (s)" and "Frequency (kHz)", and the colorbar label stands for "Power/Frequency (dB/Hz)". I don't see any visible difference between the two figures.
In addition, you can use a smaller FFT size. Generally, using such a long FFT size means that you need an extremely high frequency resolution.
$$
\Delta f = \frac{f_s}{\mathrm{NFFT}}
$$
where $\Delta f$ is the frequency resolution, which is 0.0168 Hz when $f_s$ = 78 Hz and NFFT = 4639. While a super high frequency resolution means a super low time resulution with the fact that $\Delta T = 1/ \Delta f$. In your case you only have 10 time sample results.

For the second question, 0 Hz means the direct current (DC) component. Maybe you should add some information about what kind of signal are you analysing.


Frequency range for set of all discrete-time sinusoids Consider discrete time sinusoids of the form $$x[n] = \cos(\omega n) \ ,$$ where $n$ is an integer. What frequency range of $\omega$ would constitute all the possible sinusoids? I'm thinking that $0\leq\omega<\pi$ would suffice, but my teacher tells me it's $-\pi\leq\omega<\pi$.

assuming $\omega$ is real, if the form of $x[n]$ is how you've defined it above, your thinking is correct.  doesn't matter if it's $+\omega$ or $-\omega$, the sinusoid behaves exactly the same.  if instead, the sinusoid had form

$$ x[n] = e^{j \omega n}  \ ,$$

then your teacher would be correct.


Yet another fft question - a vocoder Lets assume a pipeline X:

            in_buffer1 * window -> fft -> car_to_pol -> pol_to_car -> ifft -> out_buffer1
          /                                                                               \
in_signal                                                                                out_buffer1 + out_buffer2  -> out_signal
          \                                                                               /
            in_buffer2 * window -> fft -> car_to_pol -> pol_to_car -> ifft -> out_buffer2

where in_buffer1 and in_buffer2 overlap in time by exactly 50% and the window is i.e. a hamming window with an offset matching the buffer.

the out_signal should be just a reconstructed in_signal with a delay depending of the buffers lengths with no frequency leaks.

However lets then consider a following pipeline Y for a single buffer with two sifnal sources and one result

in_signalA(n) ->  in_bufferA(n) * window -> fftA(n) -> car_to_polA(n)  ->     imgA(n)     ->   pol_to_car -> ifft -> out_buffer(n)
                                                                       \                   /
                                                                        realB(n) * realA(n)
                                                                       /
in_signalB(n) ->  in_bufferB(n) * window -> fftB(n) -> car_to_polB(n)

Then apply the same overlapping oerlapping as per the pipeline X where the n is the number of the pipeline, i.e. 0 - no offset, 1- 50% offset,
This should produce a 'vocoder' effect.

However I do get cracks on the output signals which suggest there is a frequency leak I'm not compensating for. Any idea what I'm missing here?

The solution turned to be the right number of overlapping input buffers.
To get the source into fft and then back we need to apply windowing on the sources to reduce frequency leaks. Reconstructing such signal without any modification works just fine. However as @applesoup poined out certain operations in the frequency domain introduce their own discontinuities at the edge of output buffers so the output buffers must be windowed as well before combining into the output signal which in turn introduces inconsistency to the amplitude of the signal as in effect we need to apply same window twice. The solution to this problem is to double the number of buffers and change the overlap to 25% (of looking from the other side, to 75%) so each sample is reconstructed from 4 windowed output buffers which compensates the double windowing of each individual buffer.


Computing response of LTI system in terms of its step response I have solved a problem, Kindly help me in determing if I solved it correctly.
I will post the question and my own working below.

Problem

My working

HINTS:

$\cos(2)-\cos(2)=0$
What remains is $y(t)=\cos(t-1)-\cos(t)$. Can this really be the response of a causal system to an input signal that starts at $t=0$? Maybe you forgot some important detail in your answer.
The correct answer must consist of $3$ different expressions for the time intervals $t<2$, $2<t<3$, and $t>3$. Your answer corresponds to one of those three expressions, so you need to add the other two, and the corresponding time intervals.


EEG Power at specific time points I'm in the process of building features for an EEG machine learning problem and have been doing some reading about what features may be useful. Turns out, that alpha band (8-12Hz) power can make a good distinction between my categories. The nuance is that for 1 category, alpha power rises and then falls over the time series, but for the other category alpha power begins low and then rises

I have a butterworth filtered signal which gives me a fairly clean alpha wave, but is it possible to determine power of this band at every time point over the series?

If I can get away with doing the calculations on just the filtered signal then that'd be ideal, without having to resort to running an FFT (the eventual goal is to do this online and my concern is speed)

Could it simply be the squared amplitude at every time point? Or do I need to run a specific FFT/time-frequency type analysis for this?

I've been trying to find some Python code for this but haven't had much luck because I mostly am not sure the best way forward

The Goertzel algorithm allows you to compute individual terms of the Discrete Fourier Transform, and is more efficient than the FFT. However, if you wish to later compute the spectrum of other bands besides alpha, Goertzel will turn out to be more computationally expensive than the FFT.

Unless you are seriously limited by computational power, in principle, you can compute the Fast Fourier Transform in quasi-real-time, using a sliding window (also known as the Short Time Fourier Transform, or STFT). The Discrete Wavelet Transform is also sometimes used for real-time EEG power estimation. To find the power in the alpha band, you can look at the corresponding frequency bins of the STFT for either the original or the filtered signal.

The output of the STFT will be delayed by the length of your window (as you need all the samples inside the window to arrive before computing the power in that time bin). However, with a window length of 512 samples under 8kHz sampling rate, for example, this amounts to 0.064 seconds delay, which is fairly close to 0.0 seconds!

The STFT (Short-Time Fourier Transform) is effectively a computation of the FFT one 'chunk' or 'time-frame' of the signal at a time. Often these chunks are overlapped by some amount, such as 50%, meaning that if the window length is 512 samples, then halfway through the first frame (at the 256th sample), you already begin to process samples that contribute to the second time-frame.

The purpose of using overlapping frames is to avoid artefacts caused by taking square slices of the signal. (Computing the power of one sample at a time, like you suggest, amounts to convolving your signal with a square window, which is equivalent to convolving the frequency representation of your signal with a sinc function. This causes spectral leakage, or 'ringing' artfeacts. Thus, it is common to use a sliding window with smooth edges, such as a Hamming window, to reduce these artefacts.)

Note that the choice of window length will affect your spectro-temporal resolution. A longer window will have better frequency resolution, but worse temporal resolution, and vice-versa. This phenomenon is known as the Gabor limit, or Heisenberg-Gabor limit.


Audio Signal Noise Filter Problem I'm currently working with audio signals and have a problem:

C = A*B + N, where

C = recorded signal from microphone consisting of:

A = known music file data played on speakers next to microphone

B = some convolution on the recorded A-sound due to speaker->mic roundtrip
(I mean the recorded signal won't be 100% the same as the audio data from the file before it is played to the speaker and recorded by mic. (Is this an impulse response?))

N = some additional noise sounds recorded by microphone

My goal: an approximate estimation as to whether there is a signal N and how loud it is.

I don't have a need for accurate data!

Additional info:
I'm working with Apple's vDSP API. I have cross correlated the signals A and C, so I have the time window in which the signals overlap.

In the overlapping window, I have both signals in the time and frequency domain.

Currently I'm helpless if, for example, a Wiener filter is the right approach and if I'm capable to apply one with my known parameters(Is a known noise required? or the impulse response of the environment?).  I tried to apply a Wiener deconvolution by dividing C/A in the frequency domain with no success.

Once more: I don't need accurate data, just a rough guess how much N is there in the signal C. Actually a SNR like measure would be sufficient.

Wiener filtering is one approach. It might even be the best approach. A Wiener filter is designed to minimize the noise (in the least squares sense) and invert the effect of the impulse response, given a known signal and a signal that is known to be tainted with noise and an impulse response. Once you have a Wiener filter, you can then compare the amplitude of the filtered signal with the amplitude of the direct signal to estimate the noise (I think).

Now I've never implemented a Wiener filter myself, so I will defer to a textbook with chapter dedicated to the subject: Advanced Digital Signal Processing (Electrical and Computer Engineering) by Glenn Zelniker and Fred J. Taylor. It's out of print, but I'm sure you can pick up a used copy for cheap somewhere. I suggest it because it's more mathematical than engineering in it's approach and that might appeal to you. Many other textbooks have info on the subject.

If a more ad-hoc approach suits you, here is a suggestion. I've never tried it, but it might work:

Create a filterbank, F, (a properly windowed FFT will work) to analyze both the known music file A and the input C. You will also need some measurement function, M. I would suggest $M(x) = | x |$, rather than $M(x) = x^2$, but both will work, and you may also want to do some smoothing over time. By comparing the results of $M(F(C))$ to $(M(F(A))$ when the noise is low, you should be able to determine the effect of B. Of course, you may not know that the noise is low, in which case, you'll have to get a bit more clever, creating some estimate based on some statistics. The goal, however, is to get to a point where you can predict the output of $(M(F(C))$ from $(M(F(A)))$, to within some bounds. The lower the bounds the better you will be at detecting the noise.

To estimate the noise, you will compare how far above those bounds the signal has actually gone. This will work well for band-limited noise, but for quieter, broadband noise, it may not work well. In that case, you may want to create an aggregate statistic based not just on how far above the bounds a signal has gone, but how many bounds have been crossed, or simply apply your measurement function directly to your unfiltered signal: $M(C)$ vs $(M(A))$.


If a square wave is a sum of odd harmonic impulses, why is it continuous in the frequency domain? A square wave is a sum of sinusoids so surely it should be represented as individual discrete impulses in the frequency domain, where all other frequencies are 0. Why instead are those intermediate frequency components not 0, and form a continuous sinc shape? I know that the continuous fourier transform of $A\text{rect}(\frac{t}{\tau}) = A\tau \text{sinc}(f\pi\tau)$, but why is this the case when you think about it in the first perspective?

As @Hilmar mentioned I think you get confused between Square wave and  Rectangular function.
In Wikipedia about Square Wave :

A square wave is a non-sinusoidal periodic waveform in which the
amplitude alternates at a steady frequency between fixed minimum and
maximum values, with the same duration at minimum and maximum.

Which its Fourier Transform is only at harmonic frequencies and its value is equal to Fourier series coefficient.
but Rectangular Function:

The rectangular functionis defined as:

which look like this:

and its Fourier Transform is a sinc function that looks like this:


What does the negative sign mean in the image frequency domain? When I transform an image from the spatial domain to the frequency domain, some frequencies are negative. What is the meaning of the negative sign?

For example:

Each DCT output bin is a (weighted) correlation against a cosine function of a certain frequency.  A negative value would represent a negative correlation, e.g. something in the image data is of the opposite phase to that cosine (e.g. maybe dark when the cosine is 1, and light when the cosine is -1, instead of vice-versa for a positive correlation).


How to obtain round peaks instead of triangular peaks when smoothing an input? I have a signal (in blue) that I would like to smooth out (in red, moving average).

I would like to get red peaks to look more round than triangular and at the same time their width would become closer to the original ones.

Do you know / can you suggest a filter that would yield such result ?

Background:

The input is the zero-crossing rate from a song, for building a colored waveform as I've tried in this question : Coloring a waveform with spectral centroid or by other means

Here's the result of the new approach:

(unprocessed)

(ideal result, cheated somehow here using a post, expensive gaussian blur)

Update

Here's the result, using @Laurent Duval answer:

Also, I still need to try all of your suggestions again as my input/initial take was buggy, surprisingly now the ZCR yields better output than using FFTW out of the box (pics on the right):

Update 2

Simple moving average (green: 1-pass, gold: 2-passes)

Use two or three passes of moving average (can be pipelined).


Confused about applying Routh Hurwitz to $s^2 +s + k$ I have a closed loop transfer function
$$G_{\rm cl}(s) = \frac{k}{s^{2} + s + k}$$
I am trying to find the critical gain of the system. From using the Routh Hurwitz criterion I get a $k = 0$.

Second order Routh Hurwitz:

$$a_0s^2+a_1s+a_2 = 0$$

$$\begin{cases}
s^2 & a_0 & a_2\\
s^1 & a_1 & 0\\
s^0 & a_2\\
\end{cases}$$

In your example

$$\begin{cases}
s^2 & 1 & k\\
s^1 & 1 & 0\\
s^0 & k\\
\end{cases}$$

Hence, $k$ should be positive.


Streaming windowed sinc interpolation/resampling: trying to understand a Rust implementation I'm working on a fork of the Rust dasp library, which is intended to be a DSP toolkit that abstracts over samples/frames/signals, and contains a number of functions and tools that work with these. One such bit of functionality is a signal interpolation feature, this supposedly allows for on-the-fly lazy resampling of a signal using sinc resampling. The Rust code that does the heavy sinc work is found here: https://github.com/RustAudio/dasp/blob/master/dasp_interpolate/src/sinc/mod.rs
However, I'm a bit confused. Looking at previous asked questions on this site and other resources on sinc interpolation and resampling, I'm seeing some yellow flags:

I can't seem to find any other "streaming" implementations of sinc resampling: they usually require a fixed-size buffer of samples as input, produce a fixed-size buffer of interpolated samples as output, and don't yield interpolated values one at a time.
Other sinc implementations are significantly more complex in terms of their required parameters. For example, the r8brain resampler requires a cutoff frequency and an attenuation level in the constructor (link), both of which are nowhere to be seen in dasp's resampler/interpolator.
I'm not seeing how the resampled result is lowpass filtered.

I'm admittedly new to the world of DSP, so I'm not sure if this code is correct, and if not, I'd like to understand why and how to fix it as part of my fork.

This looks partially similar to the BASIC code in my Q&D resampling cookbook here:  http://www.nicholson.com/rhn/dsp.html#3  which computes both the Sinc and the window (von Hann) on the fly (per sample).
But it appears that the .rs code is not determining the lowpass cutoff frequency (phi upper and lower) correctly.  It needs to know whether the resampling is up or down in frequency, so that the Sinc computation can use something below half the lesser of the two frequencies.
Computation of one new resampled sample at a time works just fine, as long as the rate change is fixed (or increases slightly), and the initial and ending transients (related to the Sinc window width (or depth?)) are ignored (or accepted as transients).


awgn function in MATLAB I'm using 'awgn(input_signal,snr_db,'measured')' function in MATLAB to add white noise to the input signal.
Here is the simplified script.
Does zero padding the input signal affect the behavior of awgn function?
clc;clear all;close all;

data = randi([0 3], 1, 2^12); % Generate data
qam_mod = qammod(data, 4);  % Modulate data
qam_mod_u = upsample(qam_mod,64); % Upsample

ch_out_interplt = awgn(qam_mod_u, 5, 'measured');
ch_out_no_interplt = awgn(qam_mod, 5, 'measured');

ch_out_downsample = downsample(ch_out_interplt,64); % Downsample

scatterplot(ch_out_downsample); title('with up, down sampling, 5 dB SNR');
scatterplot(ch_out_no_interplt); title('no up, down sampling, 5 dB SNR');

The scatterplots are different with and without zero padding.

Does that mean awgn function behavior depends on the zero padding (length of the signal)?
How does the 'awgn(input_signal,snr_db,'measured')' function calculates the input signal power?
How to make the power calculation independent of the length of the input signal?
Is there any better way to simulate the awgn noise in MATLAB?

Does that mean awgn function behavior depends on the zero padding (length of the signal)?

Yes, it does depend on the length of the signal.

How does the 'awgn(input_signal,snr_db,'measured')' function calculates the input signal power?

For an input signal of length $N$, $x[n]$, the awgn function calculates the power as $\frac{1}{N} \sum_{n=1}^N |x[n]|^2$.

How to make the power calculation independent of the length of the input signal?

The above definition defines $N$ as the length of the signal, including the zero valued samples. If you don't want to count those towards your signal length, then you needs to have an effective signal length $N_0$ that only counts non-zero samples.
sigPower = mean(abs(qam_mod_u(qam_mod_u ~= 0)).^2);
noisePower = sigPower/10^(SNRdB/10);
noise = sqrt(noisePower/2)*(randn(size(qam_mod_u)) + 1j*randn(size(qam_mod_u)));
ch_out_interplt = qam_mod_u + noise;

Is there any better way to simulate the awgn noise in MATLAB?

The awgn function is fine as long as you understand what it is doing. For situations when the signal is turning off and on, it might be easier to do the noise on your own for exactly the reason of your question. You can also make it easier by putting this all in your own function too.


The FFT of a sine wave with varying frequency and constant amplitude does not show a plateau, why? I want to use python to see what the audio response of a speaker is. This would be done by comparing the output of the speaker measured with a microphone to the given input. That means that for both of these signals that will be in the time domain I need to perform a FFT. Eventually I want to play a target curve close to a Harman curve. I tried to generate the signal to put it in a wav file and when checking if the signal does what I ask with the FFT I got some weird behaviour. So then I tried to perform an FFT on a sine wave that has a constant amplitude and goes from 21-19700 Hz in 5 seconds with a sample rate of 44100 Hz. I would like to use the welch method where I took the window to be the full length of the sample as in this case the signal I generated does not have any noise. This got me the following result:

I used the following python code
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch

samplerate = 44100
duration = 5
time = np.linspace(0,duration,samplerate*duration, endpoint= False)
frequency = np.linspace(21,19700,len(time),endpoint= False)
Const_wave = np.sin(np.pi*2*frequency*time)
f, check = welch(Const_wave,samplerate,nperseg=len(Wave))
plt.plot(f,check)
plt.xscale('log')
plt.xlabel('Frequency [Hz]')
plt.yscale('log')
plt.ylabel('Power')
plt.grid()
plt.show()

Afterwards, I tried the scipy.fft module which yielded a different but also not expected result:

for this I used the following code:
f = rfftfreq(samplerate*duration,1/samplerate)
check = rfft(Const_wave)
plt.plot(f,check)
plt.xscale('log')
plt.xlabel('Frequency [Hz]')
plt.yscale('log')
plt.ylabel('Power')
plt.grid()
plt.show()

I would expect to see a flat plateau of equal power on all frequencies, but that is not happening. I also tried to do the same with a linear sweep wav file generated by REW and got similar results, so is my expectation of a flat plateau incorrect or is there something else at play here?

Your Python code has two weak spots, rather typos, but these are not the source of your failure to "show a plateau" in PSD graph. First we review these "typos".
The line
f, check = welch(Const_wave,samplerate,nperseg=len(Wave))

must read
f, check = welch(Const_wave,samplerate,nperseg=len(Const_wave))

or else, the variable Const_wave must be Wave. Without correcting this "typo", which might appear when you copied the code to the post, your code does not compile.
Then,
time = np.linspace(0,duration,samplerate*duration, endpoint= False)
frequency = np.linspace(21,19700,len(time),endpoint= False)
Const_wave = np.sin(np.pi*2*frequency*time)

results in a frequency sweep up to 19700*2, not to 19700, but maybe you did it intentionally. If not, notice that the instantaneous frequency is the derivative of the phase. To reach a frequency of $f_{max}$ in the frequency sweep, the formulas are
$$
phase(t) = f_{min}·t + (f_{max} - f_{min})/duration·t^2/2 \\
sweep\_wave(t) = \sin(2π·phase(t)) \\
instantaneous\_freq(t) = f_{min} + ((f_{max} - f_{min})/duration)·t
$$
and you should write the code
from scipy import signal
...
fmin = 21
fmax = 19700
time = np.linspace(0,duration,samplerate*duration, endpoint= False)
Const_wave = np.sin(np.pi*2*(fmin*time + (fmax-fmin)/duration*time*time/2.0))

Back to your question why "a sine wave with varying frequency and constant amplitude does not show a plateau": the answer can be found in your choice of the Welch method's parameters. You apply the Welch method to the signal in the entire range processed with a Hann window (see the scipy.signal.welch reference), and your Hann window "leaks" into the PSD. You have two options to reach the desired "plateau".
1st: Specify the named parameter window with the array-like value, like this
f, check = signal.welch(Const_wave,samplerate,window=np.ones(len(Const_wave)))

2nd: Retaining the default (Hann) window (omitting the named parameter window), decrease the value of the named parameter nperseg to a lower value commensurate with a maximum wave period (minimum frequency) in the sweep, like this
f, check = signal.welch(Const_wave,samplerate,nperseg=samplerate/frequency[0])

you can also try to adjust the noverlap parameter value (defaulted to half the nperseg value).

In both cases, the graphs would not yet become ideal "plateaus", and to understand discrepancies--quite substantial in the second case--you should delve deeper into the Welch's method working in particular and into the study of signal processing algorithms in general.
Any way, you have to, since the audio engineering requires a quite solid background in signal processing. In the future, your involvement with more advanced AE projects may (and often will) require a subscription to the Journal of the Audio Engineering Society.
In the meantime, let me recommend you references related to the issues you raised in your question (including Harman target curves):
Transfer-Function Measurement with Sweeps by Swen Müller and Paulo Massarani
Simultaneous measurement of impulse response and distortion with a swept-sine technique by Angelo Farina
Github's project AutoEQ


Is camera magnification a constant? Magnification $M$ is defined as $\frac{-d_i}{d_o}$, where $d_i$ and $d_o$ denote the image distance and object distance, respectively. For a given camera, is the value $M$ typically the same for all values of $d_o$? That is, is it true that $\frac{d_{i_1}}{d_{o_1}} = \frac{d_{i_2}}{d_{o_2}}$ for all $d_{o_1}, d_{o_2}$ where $d_{i_1}$ and $d_{i_2}$ are the corresponding image distances?

Yes for a given camera with a fixed focus lens, sensor size and sensor to lens distance the optical magnification;
$$ M = |{ \frac{ h_{obj} }{h_{img} } }|$$
which is object vertical length divided by its image length on the sensor (a number less than unity, hence a reduction,  unless a macro image is taken with  a macro lens wihich has a true magnification) will be a constant for most of the proper object distance range.


DFT - Removing window effect in spectral domain with convolution I was thinking about the DFT windowing subject and a thought came to my mind. A DFT will yield the spectrum of a signal convoluted with spectrum of the window used, therefore having a main lobes and side lobes.

I figured it would be possible to remove the window effect on the spectrum of the signal by convoluting again both the signal and the window  spectrum magnitude, and it did indeed works as you can see on the following image.

Left is the original spectrum generated with a hanning window. Right is the spectrum convoluted by the DFT of a hanning window. Top is the Spectrum itself, bottom is MATLAB findpeaks result.

I never read anything regarding this technique, but I am pretty sure I haven't invented anything there. So I am wondering if there is a benefit of doing this processing on the spectrum or if there is a downside to it that I am not seeing.

From what I see, this could help peak detection as we can see on the previous image. Also, it looks like the spectrum is a little bit distorted as we can see on the 2 following images. :

Where the blue graph is the spectrum and the red graph the post-convoluted spectrum.

Any thought about this?
Is there a problem that could arise from this post-FFT convolution?
Any paper that treats the subject?

EDIT

You can find a script here that will generate the following graph:

There is, in fact, a downside to your suggestion: the signals you showed are all very clearly divided into their frequency components, but in general, real-life signals tend to be more noisy.

Depending on the application, you'd want as much leakage dampening (higher main lobes/smaller side lobes of a signal frequency in the transformed spectrum) or, as another example, the narrowest possible main lobe.

In your plots, it is visible that smoothing the magnitude spectra with a window does the opposite of these: the main lobes get wider and smaller, while leakage products from transforming a finite time signal gain in power. Applied to noisy signals, this would yield a significant disadvantage.

Your suggestion, however, is still pretty useful for peak identification!


Help with a product of a sum of infinite series i would like some help with this:

given:
$$
h[n]=2^{-n} u[n]
$$
calculate the series:
$$
 \sum_{r=- \infty }^ \infty  h[n+4r]
$$

let me start from the general form
$$ \sum_{ k = M}^{ N } r^{k} = \frac{r^M - r^{N+1}}{ 1 - r}$$

then if you multiply $r^k$ by a unit step $u[n]$ the above series become

$$ \sum_{ k = M}^{ N } r^{k} \cdot u[k] = \sum_{ k = 0}^{ N } r^{k} $$

As you can notice the lower index become zero, because the unit step is equal to one only if $k\geq0$, otherwise is equal to zero.

Therefore in this specific case we have $u[n+4r]$ that is equal to one if the index $n+4r$ is greater or equal to zero, then:
$$ r \geq - \frac{n}{4}$$

At this point I think that you have all the necessary information to resolve the series


Non Uniformity in MFCC Training Data I'm currently working on Person Identification, and using MFCC as the feature to my classifier. I have a database full of tagged(with labels) voices, and I use a CSV file to access the respective voice files. I generate the coefficients by dividing each voice file into equal frames of 1024 bytes & then, I find the coefficients of each frame. My problem is that all the training data isn't of the same size. Thus, the total number of coefficients vary from file to file. This dramatically affects the performance of my system. Is there any way around it? I've tried treating the whole file as a single frame, but the efficiency hasn't been very good. Any help would be appreciated.

Definitely don't use the whole file as a single frame.
As Alexander says, you want to use 25 millisecond frames with a 10 second overlap/hop (the vocal tract is stationary over this duration, so this is a common choice in speech recognition).
Each frame of audio makes a certain number of MFCCs (say 20). This number is constant for each frame of audio. But a file of different duration will, of course, have a different number of frames (and thus MFCCs). The size of your features (20) should always be uniform, but the number of observations (how many frames) you get from each file usually will not be, and that is okay.
Once you have extracted these acoustic-features, you need a machine learning algorithm that can "learn" what MFCC features are characteristic to a particular speaker. This is a very broad topic, but I can point to you this project: https://github.com/ppwwyyxx/speaker-recognition where you can see some of the machine learning algorithms used.

Good luck!


Spectrogram, time vs frequency localization, and length of signal So I'm having some trouble conceptually and am hoping somebody can clear something up for me. This is long and I apologize for that but I want to provide as much information as possible so that the source of my confusion might be exposed.

I have a signal with a sampling rate of $74767 { samples \over sec}$. I want to take the spectrogram in matlab of this signal but I'm really only interested in the content in the first $1\over2$ second. Now, if I throw out everything after 0.5 seconds and attempt to compute the spectrogram using the time-frequency tool box I've posted the link to below. (it's required I use this one) I run out of memory. To solve this I decided to keep only every $X$ samples so that I'm working with less data. This will effectively, I assume, lower my sampling frequency by a factor of $X$. That is to say then that I first throw away everything after 0.5 seconds, and then keep only ever $X^{th}$ sample. My total signal length is then equal to: $$N = {74767 \times 0.5 \over X}$$

My new sampling frequency should therefore be twice this value as that would be, hypothetically, the number of samples I would have collected in one second.

$$ f_s = 2 \times ({74767 \times 0.5 \over X}) $$

I'm having two issues however with this.

The first is that I'm slightly confused about the frequency resolution and number of bins for my spectrogram.

I have to pass to the spectrogram calculation the number of frequency bins I want to calculate. I always pass it the length of the signal, which as discussed above is $N = {74767 \times 0.5 \over X}$. Now, I know that to plot this with the proper scale I need to know what my frequency resolution is per frequency bin. This is going to always be 2 because the frequency resolution is: $$ freq_{res} = {f_s \over N} = {1 \over T}  $$

As I've shown above, I reasoned this out, but it seems that if I increased $N$ and thereby sampled at a higher rate, then I would have more samples per second in a half-second time frame and my frequency resolution should be higher! I mean after all, if I'm looking from 0 to 0.5 seconds, and didn't have to divide by $X$ because of the memory limit, then why wouldn't I have a higher frequency resolution? I would have more samples!

The second issue I'm having trouble with is understanding the resolution of the spectrogram with respect to the number of frequency bins given by $N$. Clearly I can choose $X$ to vary $N$. As a test I included $N=3000$ and $N=8000$ below. As $N$ gets larger I start to lose precision in the frequency domain. This makes sense as I know the uncertainty principle tells us that there is a balance between precision in the frequency vs time domain. However, how do I know what a good value of $N$ is? My highest frequency content I really care about is the 25th harmonic ($60 Hz \times 25 = 1500 Hz$). Does that mean that I should perhaps choose a value of $X$ that makes my sampling frequency twice the Nyquist (in my case that would be $3000 Hz$). Doing that I get:

$$ 3000 Hz = 2 \times ({74767 \times 0.5 \over X}) $$

$$ \therefore X = {{74767 \times 0.5} \over 1500 } \approx 25$$

However, when I plot that my localization in the frequency domain is great, but my localization in the time domain is pretty terrible, I've included a plain-text link screenshot of it below as well.

What then, is a good metric to decide on my localization in the time domain? And how do I reconcile the seeming contradiction above? Thank you for any help and/or advice.

N = 3000

N = 8000

N = 1495. X = 25

Plain-text Link for image found here: i.stack.imgur.com/1HLTL.png

(I don't have enough Karma to post more than 2 links).

TFT Toolbox: tftb.nongnu.org/

The frequency resolution stays the same because increasing the sampling rate increases the maximum frequency that can be resolved. For example, doubling the sampling rate would double the maximum frequency.
Adjust the size of the spectrogram window, and set the number of FFT points based on that. The command is

[S,F,T] = spectrogram(X,WINDOW,NOVERLAP,NFFT,Fs)

X is the signal
WINDOW is the size of the window, try 1000 for starters
NOVERLAP is the amount of overlap as the window is slid along. signal length / (WINDOW-NOVERLAP) is approximately how many time points you will have.
NFFT is the number of FFT points. The maximum I would use is half the WINDOW value, say 500.
Fs is the sampling frequency. Input this and F will come out as the frequency points and T as the time points with no need for your own calculations.

You should not be running out of memory - I routinely calculate spectrograms on data with millions of points, you just have to tweak the input parameters. For example, putting WINDOW = 1000 and NOVERLAP = 999 will give you one time point for every sample and the output spectrogram will be, in your case, 74767 wide! Try the above settings and see how it goes.

If you make WINDOW large your frequency resolution will increase but you time resolution will decrease.


simulating multipath via tap delays in FPGA I have a Xilinx FPGA running at 200MHz with an ADC and DAC on board.  I am taking a 2462MHz sine wave (pretending it is 802.11b chan 11), mixing it down to 40MHz putting it through my FPGA, mixing it back up and looking at it on a Spectrum Analyzer.

I have 4 taps (255 clock delays deep) to play with in my FPGA and am currently only using 2 of them.  I worked out the phase delay on each tap index by saying that the phase for a particular index (i) is phase=i*(2*pi)(2.462e9/200e6).  So index 1 has a phase of 0.13(2*pi) and index 250 has a phase of pi.

Armed with that, I put my signal through with no delay on the first tap, and then delayed the signal on my second tap by 250 clocks (so it should be 180 degrees out of phase).  Yet the signal doesn't disappear on the specAn.  I haven't found a tap that does make the signal cancel out, so that was where I was curious what I was doing wrong.

I was hoping there was something simple I am missing, as opposed to finding out it can't be done without implementing some sort of fractional delay filters (or something to that effect).

Perfect destructive interference is possible in theory, but very difficult to do in practice without feedback to dynamically adjust.  The reason for this is that there are three things you have to get just right to eliminate the tone: frequency, time (or equivalently, phase), and amplitude.  If you get any of those wrong there will still be a tone (or rather, two tones) left over.

Frequency  No matter what you do at the baseband level (i.e. in the FPGA), the frequency will not match if your upmixing does not exactly match your downmixing.

Time/Phase  This is a lot more complicated than you realize.  First, your phase-per-tap calculation requires that the sinusoid be exactly 2.462GHz and that your FPGA frequency be exactly 200MHz.  Neither one is.  You can never get an exact frequency, there is always some error and jitter.  If you take the derivative of the phase-per-tap with respect to the sinusoid frequency (2.462GHz), you will find that the phase per tap is very sensitive to sinusoid frequency errors.

You have also not taken into consideration the time delays through your ADC, DAC, and connectors, passive analog components, etc.  They all matter too.  Long story short, I would not try to figure out the correct delay, I would empirically see how the jammer's phase relates to the original sinusoid and adjust from there.  Yes, you will almost certainly have to use fractional delays.

Amplitude  Even if you get the frequency and phase right you will only be changing the tone's magnitude unless you match its amplitude.

I hope that see that you have taken on a very difficult task.  Since you are only trying to test simulated delays, I would find an easier way to do that.  I would look at using chirp signals, and then vary the delays to move the responding chirp in time.


What is the principle behind my voice echoing perpetually when 2 separate devices are used for playback and recording on video calls? I had to use 2 separate devices for a video call, desktop for video/audio output and smartphone for recording, since the desktop microphone was not working properly. Unfortunately for me, when I started speaking, my first few words got echoed perpetually, and my following words were ignored by the microphone. What could be the reason behind this phenomenon?
More importantly,
Why does this phenomenon not occur in video calls if the same device is used for recording and playback? How is this phenomenon avoided?

On each individual device, the speaker output can get subtracted from the microphone before it gets sent to other locations. This prevents others from hearing themselves through your microphone. When using two devices in within audible range of each other, the devices cannot subtract the speaker audio from the microphone audio because the information path does not exist. The delay time is the determined by the network delays between the two devices and if the device with the speaker audio is very loud, you can get a gain of more than 1 to the device with the microphone. If the device with the microphone is using automatic gain control and you are not the loudest thing in the room, then you will eventually only hear the loudest thing in the room. In your case, the echo must have been very loud.


Image processing and the Fourier Transform I am trying to understand how and why the fourier transform is used in image processing / computer vision. Below is what I have gathered so far. Would my understanding of it be correct? If not, could somebody explain it to me in simple, plain english? Or, does anybody have anything to add to it? Last but not least, could somebody explain the "discrete fourier transform"?

The fourier transform decomposes an image into its sine and cosine components. Put simply, sine and cosine are waves starting at a minimum and maximum respectively. In the real world, we can't tell whether a wave that we observe started at a maximum or minimum point, and therefore we can't really distinguish between the two. Therefore, sine and cosine are simply referred to as sinusoids.

When applying the FT to an image, we transform it from its spatial domain into a "frequency domain", which in essence is the image represented in terms of its variation in colour and brightness over time (well, not time, but space. That is, over a number of pixels).

EDIT: Why would I use the Fourier Transform? And what are its benefits over other methods? For example, one application in literature is in shape recognition or noise elimination. In basic terms, how could one go about shape recognition using the FT?

At a conceptual level, the Fourier Transform tells you what is happening in the image in terms the frequencies of those sinusoids. For example, if you have a picture of a plain wall, the values of the pixels change very little as you go from left to right or from top to bottom.  In the frequency domain that means that your image contains low frequencies, but no high frequencies.

On the other hand, if you have a picture of a picket fence, then the values of the pixels change all the time as you go from left to right.  So in Fourier domain you have high frequencies in the X direction, but not in the Y direction.

Finally, if you have a picture of a checkerboard, then the pixel values change a lot in both directions.  Thus the Fourier transform of the image will have high frequencies in both X and Y.

Because the Fourier transform tells you what is happening in your image, it is often convenient to describe image processing operations in terms of what they do to the frequencies contained in the image.  For example, eliminating high frequencies blurs the image.  Eliminating low frequencies gives you edges.  And enhancing high frequencies while keeping the low frequencies sharpens the image.

FFT is used extensively in image processing and computer vision. For example, convolution, a fundamental image processing operation, can be done much faster by using the FFT.  The Wiener filter, used for image deblurring, is defined in therms of the Fourier transform.  But more importantly, even when the Fourier transform is not used directly, it provides a very useful framework for reasoning about the image processing operations.

Steve Eddins, one of the authors of "Digital Image Processing with MATLAB", has a whole series of blog posts on the Fourier transform and how it is used in image processing.


Plotting Real Amplitude of Sine Wave in FFT Consider a sine wave having an amplitude of 9 and a frequency of 0.3 Hz. I wish to perform an fft and plot the frequency on the x-axis and the real amplitude (i.e. 9) on the y-axis. I am using the following MATLAB code:

Fs = 10;
data = xlsread('Signal.xlsx');
Y(1)=0;
Y = fft(data);
N = length(Y);
freq = 0:Fs/N:Fs/2-Fs/N;
freq = freq';
amplitude = abs(Y(1:floor(N/2)))/floor(N/2);
plot(freq,amplitude);
xlabel('Frequency [Hz]');
ylabel('Amplitude');

When running the code, an amplitude close to 9 is obtained at a frequency of 0.3 Hz. In reality, I have a sinusoidal signal of a rotating motor shaft. I am measuring the angle using an encoder. The data can be found at:

http://www.dropbox.com/s/pe3g44b48ugauqp/Signal.xlsx?dl=0&m=

I know that the amplitude of the sine wave should be about 0.6 m but when conducting an fft, it is indicating that the amplitude is around 0.4 m. I am sampling at 10Hz and have about 185 000 data points. When trying to reduce the data, the amplitude of the sine wave is approximating the desired one i.e. 0.6 m. As far as I know, the more data points you have, the better will the resolution be when conducting an fft. So why am I getting the opposite? Can someone please help me? Thanks.

Judging from your post, you have over 5 hours of data, your are using to calculate the spectrum. This is a lot of data points to process with FFT. Let's focus on your code and modify it a little bit to generate the $0.3\; \mathtt{Hz}$ sinusoid with an amplitude of $9$:

Fs = 10;
t = 0:0.1:1000;
A = 9;
f0 = 0.3;
data = A*sin(2*pi*f0*t);
N = length(data);

Y(1)=0;
Y = fft(data);

freq = 0:Fs/N:Fs/2-Fs/N;
freq = freq';
amplitude = abs(Y(1:floor(N/2)))/floor(N/2);

subplot(211)
plot(t, data)
xlabel('Time [s]')
ylabel('Amplitude')
title('Time domain signal')
grid on

subplot(212)
plot(freq, amplitude);
xlabel('Frequency [Hz]');
ylabel('Amplitude');
title('Amplitude spectrum')
grid on

This gives us what you expect to obtain.

Now let's assume that your $0.3\; \mathtt{Hz}$ signal is not always present (sometimes there is only noise present). I suspect this is how your data looks like, especially over so many hours. I will approximate this situation by adding the random noise after half of a sinusoid duration: data(5000:10000)=randn(1,5001);
 (consider it as someone turned off the machine halfway through the experiment). Obviously you can make it even more intelligent by having signal for 1 second, followed by 1 second of noise, and so on.

How it will affect your amplitude spectrum? Obviously you won't get $9$ as an amplitude because there is only half of the energy present before at $0.3\; \mathtt{Hz}$. Thus you will observe something like:

If you want to detect a true amplitude, then I suggest you to use the Time-Frequency analysis tools. Easiest one is the Short-Time Fourier Transform. You will be able to observe the amplitude of your component at a given time, and extract it later on.

After analysing your signal indeed you can expect the peak with an amplitude of $0.62$ at $0.2\; \mathtt{Hz}$ (not as you claiming $0.3\; \mathtt{Hz}$):

You problem relates to the leakage as you get lot's of frequency bins and likelihood of your frequency falling exactly between get's lower. Indeed, to many is also not a good idea. Below you can observe two spectra for two synthesized sinusoids with same amount of samples as your dataset and amplitude $0.62$. Frequencies are $f=0.2\; \mathtt{Hz}$, and $f=0.2003\; \mathtt{Hz}$. You can observe that for a perfect match, amplitude is estimated correctly, although for $0.2003\; \mathtt{Hz}$ you get $0.41$.

One workaround to deal with that problem in your case, would be to apply windowing to your signal. Although it will decrease amount of leakage, but still it won't remove it absolutely. Below you can observe example for the Hamming window.

In my opinion better solution would be to use less samples. For example if you expect your frequency to be around $0.2\; \mathtt{Hz}$, then you can choose $100$ samples (I know - very crude), which will give you resolution of  $0.1\; \mathtt{Hz}$ and component will fall exactly inside. Here is the example:

One last remark, regarding your code. I do not understand what your declaration Y(1)=0 is suppose to do? If it's initialization, then it is unnecessary. If you are trying to remove the DC component, then it's not a good place nor way.

In case of using other window than rectangular one, i.e. Hamming you must normalize it properly. In case of DFT you must divide your output samples by sum of all window samples. In case of 'no windowing' which is in fact rectangular window you get sum of all samples being equal to N. Putting it into code:

win = hamming(N); % rectwin(N), hann(N), ...
data = data .* win;
Y = fft(data)/sum(win);

For more theory please refer to other posts on DSP SE or take a look into this great paper.


How to Calculate the Local Gradient of an Image in MATLAB I know that the gradient function is supported by MATLAB, but the local gradient is not, so how do you calculate it?

  Local gradient is another texture descriptor, which is defined as:

  $$I_g=\sum_{z=-r}^{r} | \nabla I(x+\tau) |^2 $$

Ok,
First of all, pay attention it is a calculation per pixel using a sata from a blog.
Basically summing the Gradient norm over a block / windows.

To calculate what you submitted above do the following:

1.Define the Finite Differences filters to approximate the gradient. Moreover build the "Win Summation" kernel which sums all elements within a window.

mLocalSumFilter = ones(filterWinSize, filterWinSize);
vGradXFilter = [-0.5, 0, 0.5];
vGradYFilter = [-0.5; 0; 0.5];

2.Calculate the image gradients:

mInputImageGradX = imfilter(mInputImage, vGradXFilter, 'replicate', 'same', 'corr');
mInputImageGradY = imfilter(mInputImage, vGradYFilter, 'replicate', 'same', 'corr');

3.Calculate the norm of the gradient vector. Afterwards sum the norm of the gradient within the window.

mInputImageGradNorm = sqrt((mInputImageGradX .^ 2) + (mInputImageGradY .^ 2));
mInputImageGradNorm = imfilter(mInputImageGradNorm, mLocalSumFilter, 0, 'same', 'corr') ;

This is it.
Enjoy.


Binarization, and then thinning/skeletonization I have pictures of networks, like the one below, and my goal is to obtain the network's skeleton from processing those images.

My approach lies in two steps, first I convert grayscale image to binary image using local thresholding or Otsu method, and then a medianfilter (python function medfilt). Result is shown below.

Then, I use a thinning algorithm to extract the network's skeleton.
Here is the result of my implementation of the Zhang-Suen thinning algorithm.

Considering the quality of the first image, I'm pretty sure it's possible to do a lot better than that. I thus have two questions :

1) Taking the last image, how would you remove all the small lines perpendendicular to the ridges, and the small gaps, that are artefacts ?

2) What algorithm would you actually advise me to use for the above mentionned steps (conversion to binary image and thinning) ?

I'm mainly working in python.

Your implementation of the ZS algorithm is wrong. The correct output should be:

Here you can see the output superimposed on the original binary image:

Removing the small branches would require a little more effort and is not exactly off the shelf, but basically you could tag every segment between junction points and end points with its length and prune the ones which are too short. Then you can also take into account the "main" direction and not prune in that direction.


Custom ringtone sounds good on computer, can't be heard on phone I did a ringtone using Audacity and then used this website to convert to mp3. I sent this mp3 file to a friend who said he couldn't hear it well. I have an iPhone so I converted it to m4r. Both the mp3 and m4r files play nicely in my computer. However, when I imported it to my iPhone SE, only certain parts could be heard. Here it is if you want to try it

https://soundcloud.com/jpproducoes/test-ringtone/s-TmLZ7

Open it on your computer, then open the same link on your phone (you can use soundcloud, the result will be the same). On the phone it will sound much more distorted, and some parts won't be audible. I believe this might be due to sampling rate problems, but I am not sure. Any ideas?

Plug your phone into some large speakers.  If you can hear the tone much better, then the problem might be the small size of the speakers on the iPhone, which roll off low frequencies.

To prevent this, make sure any low pitched notes in your ringtone are made with timbres rich in strong higher overtones.

Added: If you create a ringtone from a stereo file, it's probably best if both channels contain exactly the same data.


Best reference for object size estimation in an image I'm working on a project where I need to measure the size of a 2D object on a flat surface from an image. I have total freedom to create a reference with any shape and size.

The camera I am using is the raspberry pi cam.

My question is simple: what is the best reference to use?

EDIT:

Ended up using this thanks to Peter K..

The usual way to do this is to form a calibration plate.  One way to form a calibration plate is to place known-shaped "fiducials" in a known-sized array on the surface where you want to do the measurement.

Then, knowing the measurements of the calibration plate, you can count pixels and do a map from pixels in each direction (possibly in each area of the image) to real-world measurements.

You also need to know that your camera will always be in a fixed location relative to the measurement surface, and that the lighting is good enough for you to pick up the external shape of the object.

An example of one possible calibration plate is below.


How to choose order and cut-off frequency for low-pass Butterworth filter? When using IMU data (accelerations and angular velocities), many papers use a 4th order low-pass butterworth filter with a cutoff frequency of 5 Hz. However, none of them really justify why those filters have been chosen.

Can anyone please explain how the order and cutoff frequency are chosen for a filter of this nature? My data has a sampling frequency of 125 Hz

an $n$th-order Butterworth filter in the $s$-domain has a magnitude function that is two straight lines connected with a soft corner at $\omega_0$ and -3 dB.  the straight line at the left is a flat 0 dB but the line on the right (assuming you're looking at log frequency and dB) has a slope of $-20n$ dB/decade drop.  That is the same as $-6n$ dB/octave drop, if you like your log frequency in units of octaves.  So you have to figure out if the drop is steep enough too whack the frequency components you want whacked.

$$ \Big| H_n(j \omega) \Big|^2 = \frac{1}{1 + \left(\frac{\omega}{\omega_0}\right)^{2n}} $$

$$\begin{align}
 20 \log_{10} \big| H_n(j \omega)\big| &=-10\log_{10}\left(1 + \left(\tfrac{\omega}{\omega_0}\right)^{2n}\right) \\
\\
&\approx  \begin{cases} 0 \qquad & \omega \ll \omega_0 \\ -10\log_{10}\left( \left(\tfrac{\omega}{\omega_0}\right)^{2n}\right) & \omega_0 \ll \omega \end{cases} \\
\\
&=  \begin{cases} 0 \qquad & \omega \ll \omega_0 \\ 20 n \log_{10}( \omega_0) - 20 n \log_{10}( \omega) & \omega_0 \ll \omega \end{cases} \\
\\
\end{align}$$


Finite impulse response FIR filters If H(Z) is linear phase FIR filter, then what can 1/H(Z) represent?
  Can it be causal and stable?
  Can it be stable if it is not required to be causal?

I think it represents non linear phase FIR
 filter since linear phase FIR filter have  poles at Z = 0 , but 1/H(z)
 has poles different from zero and that what would introduce non linear
 phase.
Concerning the stability, It is causal and it is stable if all
 poles lie inside the unit circle, it can be stable also and anti
 causal if the unit circle lies is inside the region of convergence .

Is it correct?

In order for H(z) to be a linear phase filter, it must zeros both on the inside of the unit circle and at the complementary locations (1/z) which are outside the unit circle. Therefore a linear phase circle has no stable causal inverse (since this would necessitate poles outside of the unit circle.)

Note that a linear phase filter can be decomposed into the cascade of a minimum phase filter (all zeros inside the unit circle) and a maximum phase filter with the same magnitude response (all zeros outside the unit circle). For this reason, only causal minimum phase FIR filters have a stable causal inverse.

Further interesting detail, the reverse filter (all the coefficients in reverse order) of a minimum phase filter IS the maximum phase filter. This is all now intuitive: the coefficients for a minimum phase filter will be dominant toward the start of the filter (resulting in the minimum delay as the signal will emerge from the filter sooner). Reversing this results in the same magnitude response but the filter will have the maximum delay as the signal will emerge from the filter later. When you cascade two filters, you convolve their coefficients, so we can also see how the cascade of a linear phase filter with a maximum phase filter results in symmetric coefficients. And it is well understood that any FIR filter with symmetric coefficients is a linear phase filter!

Example [2 1 1] is a minimum phase filter, and [1 1 2] is the reverse maximum phase filter. The cascade of the two results in the convolution of [2 1 1] and [1 1 2] which is [2 3 6 3 2] which is a linear phase filter.

If the filter has real coefficients, the all zeros must also have a complex conjugate. (But whether the filter has real or complex coefficients does not change the answer--- both real and complex linear phase filters will not have a stable causal inverse, nor will they have a stable anti-causal inverse).

Regarding causality and the ROC, for causal FIR filter the ROC is $|z|>0$ (all the poles are at $z = 0$), and for anti-causal filters it is $|z| < \infty$ (and all the poles are at $z = \infty$). In either case the requirement for zeros to be at reciprocal locations holds in order for the result to have linear phase.


Mechanics of a matrix Interleaver I am having trouble understanding exactly how the matrix interleaver. I have read the following page from MathWorks. In it, it gives the following example where "123456" is interleaved as "142536." Basically, it split every pair. Now, send the first packet as "123456" and 2nd packet as "142536" will allow one to correct any single deletions in each packet with the lower bound being 6 undecodable combinations(i.e. the same character is deleted in both packets).

Now, this is fine for a single deletion, but what if there are more than one deletion in each packet, will the matrix interleaver take this into account and generate a different pattern?
Is there a benefit to having one pattern over another? It seems interleaving is just increasing the entropy.

So, consider the following 2 patterns: 142536 and 531642. I would argue that the 2nd pattern has more entropy as every 3 bits have no adjacent characters(e.g. 531, 316, 164, 642).

We have different interleaving techniques, and matrix interleaving is one of them. But at the end all of them do one thing: interleaving is a technique to protect against burst errors (no matter how we do it).

To make it more clear, you should consider the reason a packet cannot be decoded (and is failed at the receiver). Each packet usually contains a number of codewords (for example of length $N$). These codewords are actually generated by encoding messages of length $K$ by an $[N,K]$ forward error correction (FEC) code. A given FEC  code has a limited correction capability (usually denoted by $t$) meaning that it cannot correct more than $t$ errors in each codeword.

When a packet experiences fading, a large segment of packet might get corrupted. Eventually, there might be some codewords untouched by the fading while some others are fully influenced (which is called burst error). In such case, the codewords that experience fading might have more than $t$ errors. Since the packet is only accepted if all codewords have less than $t$ errors, this results in a packet failure.

The objective by performing itrerleaving before FEC encoding and after FEC decoding is to distribute the errors evenly among different codewords. Hence, it becomes more likely for a packet to get accepted even if it experiences a "deep" fade. So Interleaving by itself is not used for error correction. It cannot change anything regarding the entropy either (since it only applies a permutation to the same random variables). It only helps to get a better efficiency from the FEC code when it faces burst errors.


confusion understanding FFT butterfly diagram?where are 12 , two input butterflies? https://www.alwayslearn.com/DFT%20and%20FFT%20Tutorial/DFTandFFT_FFT_Butterfly_8_Input.html

I was reading above link and then  i came across a confusing sentence, i highlighted it in attached snapshot. It says there are 12 , 2 input butterflies ,while i am only able to see only 4, 2 input butterflies at left most of snap in stage1  but unable to see any 2 input butterfly in stage 2 or stage 3

Please see below diagram to see the Stage 2 and Stage 3 butterflies.


Generating a timeseries with an arbitrary power spectrum I need to generate random numbers with a given PSD. To do so, I found this recipe. In case it's behind a paywall, here's now it works:

Given your power spectrum $P_k$ with $0\leq k<N$ points generated using a timeseries with timestep $h$, you define

$$A_k=\sqrt{\frac{P_k}{2Nh}}\exp{(ia_k)}$$

where $a_k$ are uniformly randomly chosen in $[0,2\pi)$. The desired timeseries is then the Fourier transform of $A_k$. The idea here is that since the PSD throws away the phase information, we simply have to generate random phases for the transform to get a possible realization of the timeseries that generated our PSD.

However, to ensure that all of the resulting time points are real, the author gives the condition that $A_{N-k}=A^*_k$, or equivalently $a_{N-k}=-a_k$.

When I implement this, however, I don't get a real signal, so I seem to be doing something wrong. I see one problem right away: for $k=0$, the FFT assumes a periodic signal, so the only way to satisfy the condition is if $a_0=0$. Then you can generate the rest of the phases as follows (in python):

def gen_phase(N):
    phase = np.zeros(N)
    randoms = np.random.random(size=N/2)*2*np.pi
    phase[1:N/2+1] = randoms
    phase[N/2+1:] = np.flipud(-randoms)
    return phase

The other problem is that it is not possible to satisfy the given condition is $N$ is even, as far as I can tell. This is a pain practically, because the fft is very slow for odd numbers (particularly ones of the form $2^n+1$...)

Can someone suggest what I might be doing wrong? The signal I get looks more or less right ion terms of the amplitude, but it has roughly equivalent real and imaginary parts, so something is wrong. The rest of the code is too simple to be the culprit, I think:

def main():
    psd = read_psd()
    N = len(psd)
    phase = gen_phase(N)
    h = 0.24e-6
    A = np.sqrt(psd/(2*N*h)) * np.exp(1j*phase)
    timeseries = np.fft.fft(A)
    pl.plot(np.real(timeseries))
    pl.show()

if __name__=='__main__':
    main()

EDIT: I misinterepreted. Because $A_k$ depends on $P_k$, I don't have the freedom to choose $A_{N-k}=A^*_k$ in the first place, because I don't have control over the amplitude. Any ideas?

You are reversing your phases okay.  You should be able to use an even value of N if you set the phase for N/2 to zero.

The problem you have is that your PSD spans the entire N.  Unless $ P_k = P_{N-k} $ you are not generating a DFT with complex conjugate symmetry which is what you need to produce a real valued signal.

Also, you should be using the "ifft" instead of the "fft" to do an inverse FFT.

Hope this helps,

Ced

Edit:  I didn't see your edit before I posted.  It seems you figured this out.

======================

Followup:

For even N:

def gen_phase(N):
    phase = np.zeros(N)
    randoms = np.random.random(size=N/2-1)*2*np.pi
    phase[1:N/2] = randoms
    phase[N/2+1:] = np.flipud(-randoms)
    return phase