Linear-phase crossover filters Is it possible to design linear-phase filters that sum to a flat frequency response?  If so is it practical to use them in real-time audio processing for as many as 10 bands?
My experience has only been with Linkwitz-Riley IIR filters, but I would like to explore the possibilities of linear phase or minimum phase filters.
From my initial research it looks like the frequency sampling method would result in ripples and wouldn't sum to a flat response (especially across several bands).

Is it possible to design linear-phase filters that sum to a flat frequency response?

Yes, of course.

If so is it practical to use them in real-time audio processing for as many as 10 bands?

That depends on your definition of practical. The problem at audio is that FIR filters at low frequencies get really long. It depends on the lowest frequency, desired steepness and amount cross band rejection. but a few thousand taps is typical. These can be fairly efficiently implemented using an FFT based scheme like Overlap Add, but that introduces a lot of latency: typically it's twice the filter length. That's prohibitive for many real time applications.
If you want to build a standard 10-band audio equalizer: this can be done very efficiently with cascaded biquads where each biquad is a "peakingEQ" filter and simply recalculate the filter coefficients whenever a band changes gain.
Parallel IIR bandpass filters are difficult to phase manage but it's possible. Odd Butterworth and Linkwitz Riley filter sum to a flat magnitude response but NOT a flat phase response, so you need to put compensation allpass filters in the parallel paths


Given a -3dB/ octave filter that makes Pink Noise, how can I make a +3dB/ octave filter that makes Blue/ Azure Noise? I learnt from here that the following filter can turn white noise into pink (-3dB/ octave).
  b0 = 0.99765 * b0 + white * 0.0990460;
  b1 = 0.96300 * b1 + white * 0.2965164;
  b2 = 0.57000 * b2 + white * 1.0526913;
  pink = b0 + b1 + b2 + white * 0.1848;

How would I turn white noise into blue/ azure noise (+3dB/ octave)?  Is it possible to do some simple 'inversion' on this filter or can you please give me a tip as to how I can calculate the filter coefficients.
EDIT
I probably should have mentioned that I need a C++ implementation but Robert's 2nd answer has me close.
Right now I do this in the header:
float *state = nullptr;
The in the implementation file:
state = new float[0.0]; in the constructor then inside the actual loop I take my white noise and do this:
float first = first_order_filter(whiteNoise, 0.99572754, 0.98443604, state);
float second = first_order_filter(first, 0.94790649, 0.83392334, state);
float third = first_order_filter(second, 0.53567505, 0.07568359, state);
out1 = third;

I am expecting pink noise as I have not yet swapped the Poles and Zeros (if I had then I would expect blue/ azure noise) but currently I get what looks and sounds like white noise.

You need to invert the filter, i.e. flip the poles and zeros.
The implementation that you reference is fairly awkward and will require a decent amount of math work to invert: you need to write the Z-transform for each first order section, add all the fractions into a single fraction and calculate the zeros of numerator polynomial.
An easier way would be to use this implementation https://ccrma.stanford.edu/~jos/sasp/Example_Synthesis_1_F_Noise.html and simply swap $A$ and $B$


How to extract the Freezing Index from Accelerometer data? I'm working on a project to detect Freezing of Gait episodes for Parkinson's Disease patients, and according to this paper and other ones, extracting the freezing index which is "The power in the “freeze” band (3-
8Hz) divided by the power in the locomotor band (0.5-3Hz)" will be a good feature.

Here is the explanation: One standard feature which is extracted from the raw signals is the Freezing Index (FI), defined as the ratio between the power contained in the so-called freezing and locomotion frequency bands (3-8 Hz and 0.5-3 Hz respectively). This feature is convenient since it requires only FFT-computation.

But I don't know how to implement it in Python.

I have a dataframe like this:

Then, I did something like this to extract features from the sensors time-series data:

win_size=200
step_size=40
for col in df.columns:
    if col != 'Label':
        df[col + '_avg'] = df[col].rolling(win_size).mean()[step_size - 1::step_size]

Now, I would like to extract the Freezing Index, How can I do this? and I need someone to explain it to me because I don't fully understand it!

I found this article is very useful. There is a function in this article that computes the average power of the signal in a specific frequency band. And since the Freezing Index is the ratio between the power contained in the so-called freezing and locomotion frequency bands (3-8 Hz and 0.5-3 Hz respectively), we can use this function to get the power in each frequency band and divide them.

Here is the function:

def bandpower(data, sf, band, window_sec=None, relative=False):
    """Compute the average power of the signal x in a specific frequency band.

    Parameters
    ----------
    data : 1d-array
        Input signal in the time-domain.
    sf : float
        Sampling frequency of the data.
    band : list
        Lower and upper frequencies of the band of interest.
    window_sec : float
        Length of each window in seconds.
        If None, window_sec = (1 / min(band)) * 2
    relative : boolean
        If True, return the relative power (= divided by the total power of the signal).
        If False (default), return the absolute power.

    Return
    ------
    bp : float
        Absolute or relative band power.
    """
    from scipy.signal import welch
    from scipy.integrate import simps
    band = np.asarray(band)
    low, high = band

    # Define window length
    if window_sec is not None:
        nperseg = window_sec * sf
    else:
        nperseg = (2 / low) * sf

    # Compute the modified periodogram (Welch)
    freqs, psd = welch(data, sf, nperseg=nperseg)

    # Frequency resolution
    freq_res = freqs[1] - freqs[0]

    # Find closest indices of band in frequency vector
    idx_band = np.logical_and(freqs >= low, freqs <= high)

    # Integral approximation of the spectrum using Simpson's rule.
    bp = simps(psd[idx_band], dx=freq_res)

    if relative:
        bp /= simps(psd, dx=freq_res)
    return bp

Then, I created this simple function to return the FI:

def freeze_index(data, sf, band1, band2, win_sec=None, relative=False):
    fi = bandpower(data, sf, band1, win_sec, relative) / bandpower(data, sf, band2, win_sec, relative)
    return fi

And here is how I called it in the rolling window function:

df[col + '_fi'] = df[col].rolling(win_size).apply(freeze_index, args=(fs, [3, 8], [0.5, 3], win_size,))[step_size - 1::step_size]

I hope it's the right solution, and I wish it helps.


Differences between low-pass, band-pass, notch filters I'm studying image processing and I know that there are different filters you can apply. I don't understand the difference between the low-pass/bandpass/notch filter (or high-pass/band-reject/notch). I tried to draw them but I don't know if I did it correctly:

From what I have realized it is that the difference between a low-pass filter and a band-pass is that the first is symmetrical with respect to the y axis while the second is not. Therefore, a low-pass filter passes the lower frequencies and resets those higher than a certain threshold (cutoff frequency). Instead, a band-pass filter passes the frequencies included in a certain range and cancels the others. It's correct?

But I have not really understood the purpose of the notch filters. I realized that there are notch band-pass and band-reject notch. In addition the book I'm studying (Digital Image Processing - Rafael C. Gonzalez) says that a notch filter is composed of the product of two high-pass filters. Just I don't understand how they work, how they are made and what they do.

The book also talks then about optimal notch filters that is used where there are several interference components in the image in the Fourier domain. Again I haven't understood how it works.

Could someone help me? Thank you very much!

So the difference between a notch filter and a band-pass is that I use the notch when I don't know what are the frequencies to pass (or delete) but I get it visually analyzing the image of the Fourier transform, right?

However, this is what there are written in my book:

First,note that your drawn filter responses are one dimensional filters.

Notch filter, selectively suppresses some frequency bands that are not of interest. One of the well known applications of notch filter in image processing is image denoising. As an example, consider the following image, which is corrupted by periodic noise (the tiled parallel bars).

.

If we take Fourier transform of the image, it will be something like this:

Based on our prior observations and knowledge we know this image should not have larges peaks other that central section, so those other peaks are the noise and are responsible for the parallel bars in the image. To remove those, we can apply a notch filter which allows all frequencies to pass, except for frequency components of the periodic noise.

  a notch filter is composed of the product of two high-pass filters

I would say it is product of a low-pass and a high-pass filter.


Given the Graph of a Fourier Series $\sum c_k e^{2\pi ikx}$ Find the Graphs of $\sum c_{3k} e^{2\pi ikx}$ and $\sum (c_k)^2 e^{2\pi ikx}$ Define a 1-periodic function on $\mathbb{R}$ by:

$f(x) :=$
$\left\{\begin{matrix}
1 & if & 0<x<\frac{1}{10}\\
0 & if & \frac{1}{10}<x<1
\end{matrix}\right.$

with Fourier series $f(x) = \sum_{-\infty}^{\infty}c_k e^{2\pi ikx}$.

I'm trying to find the graphs of the following two Fourier series:

$\sum_{-\infty}^{\infty}(c_k)^2 e^{2\pi ikx}$ and $\sum_{-\infty}^{\infty}c_{3k} e^{2\pi ikx}$.

Basically my strategy with the others has been to find a way to get these Fourier series into the form $a\sum_{-\infty}^{\infty}c_k e^{2\pi ik(bx)}$ for $a,b\in\mathbb{C}$.  But I haven't been able to do it with these, hopefully someone can help, thanks.

This doesn't require any complicated computations! Look up some tables of identities about the Fourier Series. Which operation would you have to apply to $f$ so that its Fourier series would become $\sum_{-\infty}^{\infty}(c_k)^2 e^{2\pi ikx}$?

For example, we know that derivating in the time domain is equivalent to multiplying the Fourier series coefficients by $i k$ ; so we know that the graph of the Fourier series $\sum_{-\infty}^{\infty} ik (c_k) e^{2\pi ikx}$ is the graph of the derivative of $f$.


How do I attenuate a tidal signal in depth time series data I work with depth time series data from electronically tagged fish. When the fish spend time on the bottom we get a prominent tidal signal (approx 1.96 cycles per day). This interferes with our analysis of the data, so I have been experimenting with removing it using FFT - i.e. run FFT attenuate the bins around 1.96 cpd then inverse FFT. It nearly works, but the tidal signal varies over time (the time series can be up to a year long, with varying sample rates between 10 s and 2 mins). So in some parts the signal is removed, in others it is exaggerated. I have attempted to divide the time series into blocks and deal with each individually, but it doesn't seem to work. I work in C#. Any help and advice would be much appreciated!

The answer is "No, you cannot attenuate the impact of tide on pressure derived depth data you have already captured without extraneous information".

The reason for this is that the hydrostatic pressure derived depth measurements from the tracked fish ($depth_{fish}(n)$), is only a component of the captured signal.

The illustrative example here is a "freakishly intelligent fish" which manages to stay at one point (with respect to the bottom of the sea) and adjusts its buoyancy with respect to pressure. When the tidal potential piles up the water in a water column above the fish, it ascends. When the water recedes, the fish descends. The pressure derived depth of this fish is constant, when clearly, it has been going up and down.

The hydrostatic pressure signal on the fish body is $P_{recorded}(n) = P_{tidal}(n) \cdot P_{fishDepth}(n)$ where $n$ is some variable that indicates time.

You are FFTing the $P_{recorded}$, aiming at suppressing the $P_{tidal}$, but, multiplication in the time domain is convolution in the frequency domain. Therefore, just by having $P_{recorded}$, you have both $P_{tidal}, P_{fishDepth}$ intertwined. That's not entirely bad news but it certainly is not the best news ever.

You need to somehow fix one to recover the other and to an extent, whether with freakish frightening accuracy as above or simply casual preference, fish do regulate their buoyancy, but perhaps faster than the speed by which the tidal potential moves over the globe. Therefore, your task might become as easy as applying a low pass filter to attenuate the tidal component which is more likely to exist at the low frequencies. The speed at which the fish are moving versus tidal variation is probably your best "weapon" here. For a bit more background information, you might want to see this link which is not exactly what is asked here but is relevant. (Note: It depends very much on the scale of observation too. Fish might have specific life patterns that are modulated by the time of day. So again, their signal might "hide" in the tide because it just so happens to be driven by the celestial objects that also drive the tide.)

Now, the reason I am asking you about your accuracy specifications is because you might be able to construct some kind of model by which you could isolate the movement of the fish from tide and/or depth.

The simplest way to think about this is to say, "My sensor says 23 meters. Where is the fish? It's at XYZ location. When is the fish? It's at some hh:mm:ss time. Right, what is the depth (from the datum) at the where? It's 50 meters. What is the depth at the where and the when? (Due to tide), it's 52 meters.", Ah! So, the tide has added two meters, therefore the fish is now at 21 meters from the bottom.

It might be possible to work out a probabilistic model too. That answers the question of what is the probability of a depth value, given the prior formed by the sum of bathymetry and tidal data and the data obtained by the sensor.

In any case, you either need extraneous information (bathymetry,tide) or to make an assumption (tidal signal much slower, therefore at very low frequencies, apply a low pass, done).

Finally, please note that bathymetry does play a significant role. For a trivial example, imagine a basin that fills up with water (and fish) during the tide and then preserves water (and fish) even after the tide has receded. For those fish, the readings will always be (for instance) between 4-5 meters when around them, the depth could well be zero. So, your $P_{tidal}$ is actually a sum of bathymetry and tide.

I suspect that this is probably "bad news" (?) but I hope it helps anyway.


IIR high pass filter I'm trying implement an IIR filter. I know the IIR filter come from that equation:

and i found that for high pass filter the coefficients are:

where fc = cut freq.

here my matlab code:

fs = 1000;
fc=0.3;
r=exp(-2*pi*fc);
a = [(1+r)/2 -(1+r)/2];
b= [1 -r];
n = 0:1/(fs):20;
n1 = 1024;
x1 = cos(0.3*2*pi.*n);
x2 = cos(0.8*2*pi.*n);
%plot(n,x1);
%hold on
%plot(n,x2);
%hold on
%plot(n,x1+x2);
%aux = fft((x1+x2),n1);
%aux = aux(1:n1/2);
%f = (0:n1/2-1)*(fs/n1);
%length(f)
%legend('x1[n]','x2[n]','x1+x2')
%plot(f,abs(aux));
%hold on

y = filter([0.5759 -0.5759],[1 -0.1518] , x1+x2);
plot(n,y,'r-',n,x2,'b--',n,x1,'g-.');
grid
%axis([0 100 -1.2 4]);
ylabel('Amplitude');
xlabel('Time index n');
legend('y[n]','x2[n]','x1[n]')

the commented parts are from FFT part of program and it's working well but the filter function isn't working. Help me plz and sorry about my english.

You have the following problems:

First your $a$ and $b$ seems to be mis-named; replace b with a and a with b. Where $b$ is the numerator coefficients of the filter and $a$ is the denominator coefficients. Also note that the filter command works like this y = filter(b,a,x)...

Then also your $f_c$ parameter is not (directly) controlling the $-3$ dB cutoff frequency (in Hz) of the high pass filter that you intend to obtain. But it's controlling th eradian per sample cutoff frequency. So you bettwe use teh notation wc for that instaed of fc.

In any case, your input signal frequencies are so low that they are killed by the filter and hence you see no output. Just increase their frequencies to see them at the output. Incidentally you can also decrease $f_c$ well below $0.3$ to get a tighter highpass (a DC-notch) filter..

Based on your comment, since you want a cutoff frequency in radians/sample $\omega_c$  instaed of a analog Hz frequency, then the parameter $f_c$ yields the correct discrete-time cutoff frequency. You can view the frequency response of your filter by the line figure,freqz(b,a) . So the only problem remains is that your signakl frequencies ar too low. replace it like so

x1 = cos(0.1*fs*pi*n);
x2 = cos(0.3*fs*pi*n);

note that in your code the index $n$ has a scale of $T_s$ so the frequencies should be scaled by $f_s$ to get normalized frqeuencies in the discrete-time frequency.


Why does the addition of noise improve coherence between two signals? Taking a simple example of the coherence between two signals, each composed of two sinusoids of frequency 100 and 200 Hz, I expect the coherence to show two peaks in the spectrum located at precisely these two frequencies. The script below can be used to study the effect of adding noise An =0.1 or no noise An = 0. before taking the coherence. This produces in the case of added noise

and in the case of no noise

It is clear that with no noise there are many false peaks in coherence.

Question: Can someone show me mathematically speaking why this is the case?

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

t0, t1 = 0, 2.
fs = 2000
dt = 1./fs
t = np.arange(t0,t1-dt,dt)

#f1, f2 = 50, 50
#x1 = np.sin(2*np.pi*f1*t)
#x2 = np.sin(2*np.pi*f2*t+2*np.pi*0.)

An = 0.1            # Noise amplitude
f1, f2 = 100, 200   # Signal frequencies
p1, p2 = 0.25, 0.5  # Phase shifts (times pi)

x1 = np.cos(2*np.pi*f1*t) + np.sin(2*np.pi*f2*t) + An*np.random.randn(len(t));
x2 = 0.5*np.cos(2*np.pi*f1*t-np.pi*p1) + 0.35*np.sin(2*np.pi*f2*t-np.pi*p2) + An*np.random.randn(len(t));

fcoh, coh = signal.coherence(x1, x2, fs, window='hamming', nperseg=400, noverlap=0)

fig2, ax2 = plt.subplots(1,1)
ax2.plot(fcoh, coh**2,"k-")
ax2.set_ylim([-0.1,1.1])
ax2.set_xlabel(r"$f$ [Hz]")
ax2.set_ylabel(r"Coherence")
ax2.grid("on")

To elaborate on @user28715's answer:
Even with a minuscule amount of noise say An = 0.0000001 you will get a clean graph like your first graph. Python's implementation of the Welch method certainly uses finite-precision floating point numbers and quantization errors in the input and from the calculation steps such as windowing and Fast Fourier Transform (FFT) propagate to the results. From SciPy docs:

Cxy = abs(Pxy)**2/(Pxx*Pyy), where Pxx and Pyy are power spectral
density estimates of X and Y, and Pxy is the cross spectral density
estimate of X and Y.

The normalizing division in the equation means that the values of Pxx and Pyy can be arbitrarily close to 0 while still giving large Cxy.


How is time or sampling rate accounted for in this model? A basic model of coupled strings (eg. piano) is provided here as:

The principle is that it has two identical string simulations each formed by a delay line and LPF. The outputs of these are summed at the "bridge" in the center. The sum is multiplied by $-Hb(z)$ to get the bridge output, which is then subtracted from each of the string loops to simulate damping from energy lost to the bridge.

The bridge equations are:

Where $Rb(s)$ is the impedance of the bridge and $R$ is the impedance of the string.

Assuming I got that right, here's what I don't understand:

If you are taking the summed output of the two string sims and multiplying them by a certain constant $Hb(s)$, then subtracting this from the string loops every sample, how is time or sample rate being accounted for? Ie. If you run it at 96 kHz wouldn't you do the subtraction twice as often as if you run it at 48 kHz doubling the damping? Or otherwise at what interval are you supposed to perform the subtraction?
What does it mean when an expression uses $(s)$ vs. $(z)$? Eg. $Hb(s)$ vs. $Hb(z)$?

Thanks.

Units of s represent the Laplace transform while units of z represent the z transform. It is often much easier to transform a time domain (units of t) signal to Laplace given that it translates integro-differential equations into simple algebra. The z transform does the same thing, but takes advantage of the repetition from sampling to make the transform even easier.

"s" is any complex variable (with real and imaginary values), usually expressed $s=\sigma + j\omega$. When you limit s to only be imaginary ($s=j\omega$), the result is the Fourier Transform, which is a single-valued function versus frequency. Thus for a system, $H(s)$ would be the general transfer function of s and $H(j\omega)$ is the frequency response.

Similary "z" is any complex variable, but due to the mapping from analog to digital in the construction of the z -transform, what is the frequency axis in the s-plane (the $j\omega$ axis) becomes the unit circle in the z plane ($e^{j\omega}$), where $\omega$ here is typically normalized radian frequency (see What is normalized frequency)  with a unique range of $0$ to $2\pi$. Thus a digital system would have a general transfer function given as $H(z)$ while the frequency response would be $H(e^{j\omega}$). Conveniently, the z transform of one sample delay is simply $z^{-1}$, which hints at the great utility of using the z-transform for discrete time systems rather than the Laplace transform (which you could certainly still do instead if you want to take a mathematical beating!).

(Note: The mapping from $j\omega$ to the unit circle applies to many but not all mapping techniques: see z-Transform Methods: Definition vs. Integration Rule)

That said, to model the complete system digitally as a discrete system, you would translate the blocks that are functions of s to be functions of z. It is the process of translating, or mapping, from s to z where the sampling rate comes into play. There are several mapping techniques to choose from and a few popular ones are the matched-z transform, the method of impulse invariance and the bilinear transform. (For more on that see How/why are the $\mathcal Z$-transform and unit delays related? )   A simple example of this is the mapping of integration and differentiation operations from s to z using the matched-z transform (which by the way is identical to the method of impulse invariance for an all pole system):

Integration in time is equivalent to multiplying by $\frac{1}{s}$ in Laplace (there we see the conversion to simple algebra!). This is mapped to the following in the z domain:
$$\frac{T}{1-z^{-1}}$$

Differentiation in time is equivalent to multiplying by $s$ in Laplace (along with subtracting the initial condition which I didn't include).  This is mapped to the following in the z-domain:
$$\frac{1-z^{-1}}{T}$$

where $T$ represents the time of one sample period.

Notice the factor of $T$ which is the sample period. The differentiation example is very easy to see what is going on: $1-z^{-1}$ is a simple difference of two successive samples:

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

Which is a (not the only) discrete approximation of a continuous time domain derivative. Consider taking the derivative of any time domain function: as we sample at higher and higher rates, the difference between two successive samples will get smaller, thus by dividing by the sample duration for the case of the derivative which is also getting smaller, we normalize the result to be independent of the sampling rate.

In your case specifically, this would occur in the mapping from s to z for the continuous time transfer function given by $H_b(s)$.


Confusion about result of FIR All-pass filter design I try to design a FIR All-pass filter with random phase in the frequency domain. I am a bit confused by my result and am not sure if the reason is a programming error or a misconception about convolution and impulse response of myself. I use Python and Numpy:

import numpy as np

mag = np.ones(129) # magnitude is unitiy
phase = np.zeros(129)
random_numbers = (np.random.rand(128)-0.5)*2.0
phase[1:] = random_numbers * np.pi # phase is random except for dc offset
spec = mag*np.exp(phase * 1j) # get spectrum from magnitude and phase
coefficients = np.fft.irfft(spec) # do inverse fft to get coefficents, result is of length 256

If I do a fft on my coefficients the magnitude looks nearly right. My first question: Why is the highest frequency suppressed?

And secondly, if I use my coefficients for convolving with an impulse, I thought that should reproduce the coefficients? But it does not. Why are there these "random variations" around 1.0?

impulse = np.zeros(256, dtype=np.float64) # create impulse
impulse[256/2] = 1.0
result = np.convolve(impulse, coefficients, "same") # convolve it

Any explanations or hints are much appreciated!

The dip at Nyquist has to do with the fact that you didn't assign the value $0$ to the desired phase at Nyquist. You correctly assigned $0$ to the phase at DC, but you have to do the same at Nyquist, because a filter with real-valued coefficients has a real-valued frequency response at DC and Nyquist. So you just need 127 random numbers, not 128.

But the plot of the magnitude response of the designed filter is deceptive, because by taking an FFT (of length 256) of the coefficients you basically just plot the desired magnitude (apart from the Nyquist point), because you just invert your design process. Try to zero-pad your coefficients and take an FFT of length 1024 or more. You will be surprised to see that your filter doesn't at all resemble a very good all-pass filter. You see some indication of this in your second plot, which is more realistic.

Note that a non-trivial FIR filter (i.e. one which is not just a single delta impulse) can never have an exactly constant magnitude response, because its transfer function is a polynomial which can't be constant (unless it is a trivial polynomial). Your all-pass filter has an especially poor magnitude response because it is very hard for the filter to approximate those random phase jumps.


What are the various techniques for detecting Walls in the Building(Architectural) floor plan Images? I have to process building floor plan images to extract Walls from the structure. It is trivial in case of binary images; but the images in our case are colored, and each have different color of walls. This restricts us from using any thresholding operation as we cannot assume the Walls will always be the darkest. How possibly can I come up with an algorithm that can work  on all sort of Images?

Try using Template matching by taking a small cutout of the wall whose length is long enough so that it is not confused with other elements. Keep the threshold value high(around 0.9 or higher) and then check the bounding boxes predictions.


Whether convolution followed by down sampling equals wavelet decomposition I know that wavelet decomposition equals convolution followed by down sampling by two.I tried to test this on mat lab on haar wavelet. I write the below code and compared with the default wavedec command in mat lab. The code I wrote is shown below.

clear all
clc
h=[0.7071    0.7071];% Haar scaling filter
g=[-0.7071    0.7071];%Haar wavelet filter
x=[1:10];% Input
h=fliplr(h);% h(-m)
g=fliplr(g);%g(-m)
s=conv(h,x);
s=downsample(s,2)
d=conv(g,x);
d=downsample(d,2)
[a L]=wavedec(x,1,'haar') % for comparison

After running this program I get s =  0.7071    3.5355    6.3639    9.1923   12.0207    7.0710.
d =  0.7071    0.7071    0.7071    0.7071    0.7071   -7.0710  and the wavedec out put as a = 2.1213    4.9497    7.7782   10.6066   13.4350   -0.7071   -0.7071   -0.7071   -0.7071   -0.7071. L =   5     5    10

I dont know why wavedec and my code gives different out puts, will any one help???

Simply put, no!

  I know that wavelet decomposition equals convolution followed by down sampling by two

No again. Some linear wavelet decompositions can be turned into a cascade of convolutions followed by downsampling. But wavelets are much more general. Their discrete implementations involve intricate properties in filters, downsampling, and levels.

If you just filter a signal  by some filter, and downsample, in general, you get almost nothing but aliasing.


Time domain distortion vs Harmonics I have a periodic time domain sensor signal sampled at $500 \text{ Hz}$ whose amplitude spectrum I compute and obtain. From the periodicity of the signal, I was expecting the fundamental frequency to be around $f = 1.2\text{ Hz}$, and this is what I observe with amplitude say $A_1$. But besides the minuscule diminishing harmonics at integer multiples of $f$;  what may be the cause of having the harmonic at $f_2$ around $2.4 \text{ Hz}$ that has an amplitude $A_2\approx1.25A_1$  ?

The signal is of $29 \text{ sec.} \Rightarrow N = 14500 \text{ samples}$. Down is a section with a full period of the signal, the blue one.

Below is a section of the amplitude spectrum. NFFT = 2^nextpow2(N);
I have rectified the question.

I think this is simply a property of your signal and you indeed have two significant frequency components. I don't really know the nature of your signal, but for me, it makes perfect sense, you can distinguish two periods in your signal, as shown on the plot below:

As you can see there is period $\color{green}{T_1}$ and $\color{brown}{T_2}$ between consecutive "peaks". So by doing very rough approximations from your plot:

Corresponding frequency for $\color{green}{T_1}$ can be estimated as: $\dfrac{1}{7.1-6.3}=\dfrac{1}{0.8}=1.25 \mathtt{Hz} $
In case of $\color{brown}{T_2}$ you get: $\dfrac{1}{6.7-6.3}=\dfrac{1}{0.4}=2.5 \mathtt{Hz} $

Some quick and dirty plot of such case with two sinusoids of a given frequencies and amplitudes:

So indeed you have two frequencies present in your signal. What's more, component around $2.5 \mathtt{Hz} $ can be stronger. Please consider that your signal also has very prominent DC component (also visible at your spectrum).


Improvement of Minimum description length (MDL) estimate I earnestly request apology if this question is inappropriate for the forum. The question has two parts one technical and the other is not technical. I would appreciate any response.

Let me consider a specific use of MDL: model order estimation. For radar/sonar it is equivalent to the estimation of the number of targets (location of target is not considered here). MDL has two parts: maximum likelihood (ML) estimate of the parameter and a penalty function. Many variations exist in the literature related to the penalty function but ML estimate is hardly touched as ML is considered optimal. The model, the ML estimate and the MDL estimate all fine with large sample size; but for small sample size or high noise it is not the case as the 'model' is not accurate. The model is valid only as an expected value which is an asymptotic value.

I have an algorithm which improves the estimate of the number of targets considerably at low number of sample size or higher noise. The algorithm uses concepts of quasi-maximum likelihood.

My questions are: 1) Can we prove that a asymptotically good model, whose parameter we want estimate, for small sample size can have better estimator than ML?

The non-technical question is:

2) Anyone interested in this problem or use MDL for number of targets?

Thanks a lot.

For lower number of samples the preferred algorithm is the Akaike Information Criteria (AIC). Both are not optimal in any way, though MDL is a consistent estimator of model order.  The difference between them is the penalty term and not another estimator (other than ML)

I refer you to Petre Stoica paper on model order selection which gives a great summary of algorithms.

I am also interested to hear more about your algorithm.


Objective Noise Classification for speech enhancement I'm working on my thesis and it is about an speech enhancement system to use proper noise reduction method/algorithm for different types of noise.

I want to test and compare some filters/algorithms for different noise types added to speech and find the best filter/algorithm for each noise type(s) so the designed system would use the best algorithm when it detects the noise type.

Are there any methods to automatically detect noise type so a speech enhancement system switch the filter/algorithm? Considering noise type would change during the speech?

Any advise on noise classes and proper filter for each class?

Environmental Sound Classification might be suitable for what you are asking for. Datasets include Urbansound8k, ESC-50 and Audio set. Pretrained classifier like VGGish are available, which cover the 527 classes in Audioset.


Simple way to replicate bloom (glow) effect in image using Matlab? I have some simulated images of light bulbs in Matlab. What I want to do is replicate the bloom (glow) effect on these images. The degree of the effect (i.e. how much the light bulbs bloom) should be controllable.

Is there any implemented function in Matlab Image Processing Toolbox help me to achieve this. Or is there any simple algorithm for this. I don't want a sophisticated algorithm since I don't require the bloom effect to be precisely replicated.

I would appreciate any suggestion.

Thank you.

Work in a linear gamma (intensity-linear) color space. Create a Gaussian blurred version of the image. I don't have MATLAB to try it, but it can be done using imgaussfilt. Then add some of the blurred image to the original image, or mix the two images to taste.

Figure 1. Bloom simulated in Gimp by adding to the image its Gaussian blurred copy at 31.5 % opacity.


Why does python's scipy.signal.dimpulse introduce delay in impulse response? Consider a simple linear, time invariant system of the form:

$y_k = cy_{k-1} + (1-c)x_k$

The impulse response of this system can be computed by either dimpulse or by applying lfilter to a vector composed of a one followed by zeros:

import scipy.signal as sp_signal
import numpy as np

Ts = 1
c = 0.9
A = [1, -c]
B = [1-c]
time, imp_resp1 = sp_signal.dimpulse((B, A, Ts))
x = np.zeros(100)
x[0] = 1
imp_resp2 = sp_signal.lfilter(B, A, x)
print(imp_resp1[0][:5,0])
print(imp_resp2[:5])

which yields:

array([ 0.    ,  0.1   ,  0.09  ,  0.081 ,  0.0729])

[ 0.1      0.09     0.081    0.0729   0.06561]

Why does dimpulse introduce a one-sample delay in the impulse response?

In lfilter the transfer function is described in decreasing powers of z, as shown below copied from the python doc for signal.lfilter:

                    -1              -M
         b[0] + b[1]z  + ... + b[M] z
Y(z) = -------------------------------- X(z)
                    -1              -N
         a[0] + a[1]z  + ... + a[N] z

While the parameters B and A in dimpulse are represented in increasing powers of z as described in the doc for signal.TransferFunction:

  If (numerator, denominator) is passed in for *system, coefficients for
  both the numerator and denominator should be specified in descending
  exponent order (e.g. s^2 + 3s + 5 or z^2 + 3z + 5 would be represented
  as [1, 3, 5])

Along with the following example given for:

$$\frac{z^2+3z+3}{z^2+2z+1}$$

Contruct the transfer function with a sampling time of 0.1 seconds:

H(z)=z2+3z+3z2+2z+1
H(z)=z2+3z+3z2+2z+1

>>>>signal.TransferFunction(num, den, dt=0.1)
TransferFunctionDiscrete(
array([ 1.,  3.,  3.]),
array([ 1.,  2.,  1.]),
dt: 0.1
)

So in your case, for dimpulse:

$$H_1(z)=\frac{1-0.9z^{-1}}{0.1}$$

And for lfilter:

$$H_2(z)=\frac{z-0.9}{0.1}$$

So

$$H_1(z)= H_2(z)z^{-1}$$

And mystery solved!


Sampling a signal with varying frequency Question
I'm trying to figure out the sampling rate for my ADC to sample essentially signal essentially of the form:
$$y(t) = \sin(\max(t, \omega_{max})\times t) + n$$
where $n$ is noise.
Context
This signal is from a motor, where the rotational speed linearly ramps up to some $\omega_{max}$. To figure out the required sampling rate, the only method I know of is to use the Nyquist rate of 2 times the maximum frequency content. The above signal can be re-written as
$y(t) = rectangularPulse(t \epsilon [0, \omega_{max}]) \times \sin(t^2) + heaviside(t-\omega_{max}) \times \sin(\omega_{max}t)$
This isn't bandlimited since a rectangular pulse is a $\text{sinc}$ in the frequency domain. What I want to do is low pass filter the signal first as described here. But I'm at a loss as to figure out what I should set my cut off frequency to, as the fourier transform of the signal doesn't die off at any frequency (unless I'm mistaken).
I feel like this should be a classical problem, but I can't find anything online. If you know of how to solve this problem, or you could point me to where to look, it would be greatly appreciated!

If your frequency increase linearly, then you have a chirp signal with linearly increasing frequency.

https://en.wikipedia.org/wiki/Chirp_spectrum

There is a closed-form solution for the Fourier transform. I'm not sure how that will help you. Why not pick a sampling frequency

$\omega_s \gt 2 \omega_{max}$

Since your highest frequency is $\omega_{max}$, no point in considering lower frequencies for your sampling rate. Is this for a control-loop algorithm? In that case, I would use a more rule-of-thumb criterion.

$\omega_s \gt 20 \omega_{max}$

This makes you less susceptible to sampling delays (among other things)


Calculation of autocorrelation - does calculation "loop"? With a delay $l$, autocorrelation is defined as:
$$r_{xx}(l) = \sum_{n=-\infty}^{\infty}x(n)x(n-l) = \sum_{n=-\infty}^{\infty}x(n)x(n+l).$$

I want to calculate the autocorrelation of a signal
$$x(n) = [1,2,3,4,5]$$
with delay $l=1$. The sum starts of as
$$1*2+2*3+3*4+4*5..., $$
but what happens when I reach $n=5$? Do I start of in the beginning again, i.e is the autocorrelation
$$1*2+2*3+3*4+4*5 + 5*1$$
? Am I done there?

You didn't specify for which values of $n$ your signal is not zero, so I'll assume that happens from $0$ to $4$. You have a finite sequence $x(n)$ such that

$$x(n) = 0 \qquad \forall n\in\mathbb{Z}:n<0 \ \mathrm{or} \ n>4$$

This means that $x(5)$, $x(6)$, etc. are equal to $0$. The same happens for $x(-1)$, $x(-2)$, and so on. Thus, products involving these quantities will be zero as well.


Detect static in audio with Python I have a program running in python. It samples the demodulated FM output from an SDR. When there is no transmission, the input to the program is completely static.
How can I go about squelching the audio in a way like CSQ on a radio?
Here is what I have so far:
                    if tx_chan == 100:
                        print(np.max(abs(amplitude)))

                    if debug == True:
                        print('%s:%s '%(tx_chan,sqelch_count))

                    #Figure out if a chunk is recordable
                    if CTCSS != 'CSQ':
                        CTCSS = int(float(CTCSS))
                        for peak in peakutils.indexes(abs(np.fft.rfft(amplitude)),thres=0.2):
                            freq = (np.fft.rfftfreq(len(amplitude))*26666.66)[peak]
                            if int(freq) in range(CTCSS-2,CTCSS+2):
                                is_good = True
                                curr_ctcss = freq
                            else:
                                if sqelch_count == 0 and HighCount > 10:
                                    is_good = True
                                    curr_ctcss = 'CSQ Overide; Squelch Count: %s'%sqelch_count
                                if sqelch_count == 0:
                                    is_good = True
                                    curr_ctcss = 'CSQ Overide; Squelch Count: %s'%sqelch_count

                                else:
                                    is_good = False
                                    curr_ctcss = 'CSQ & CTCSS Fail; Squelch Count: %s; CTCSS: %s'%(sqelch_count,freq)
                            if freq > 255:
                                break
                        if is_good == False and total_count == 0 and sqelch_count == 0:
                            is_good == True```

You'd avoid having this problem by detecting whether FM is present before demodulating it.
The reason is simple: FM demodulation throws away a lot of the info that's still in the RF signal – most prominently, the actual received power – and that info is what you'd want to use to squelch.
That can (unless there's other types of transmissions than audio FM going on on the same channels) simply be done via power detection. There's other, elegant, methods working on the radio signal, but many of them depend on the way you demodulate FM.
If you're really stuck with doing it on the audio side: The output of FM-demodulating radio noise only doesn't have the same statistical properties as audio. Sadly, it mostly depends on the type of demodulation you do, again, but what should work relatively OK would be something like estimating the variance of your signal in multiple subbands. If you've got a similar variance all over the bands, that would be unusual for speech. Music, however, might still be caught with that, so be a bit careful.
TL;DR: Do the presence detection on the radio signal, not the demodulated audio.


How to, efficiently, identify telephone number from the key pressed data audio file? A very common homework question. You are given an audio file(.wav) which has the recording of the keys pressed during a telephone call. You have to identify the telephone number from this audio file when you are given the sampling frequency & the duration for which one key is pressed.

From the given sampling frequency & duration of key press, I could manage to write a MATLAB code to identify the telephone number by taking (samp_freq * key_press_time) number of samples at a time, doing fft on it & finding the bin that gives me a peak. I used WAVREAD command to read the audio file.

One problem I noticed is that, if the duration of key press is different(1st key pressed for 1 sec, 2nd key for 2 sec, 3rd key for 0.5 sec & so on), my technique of separating the data & finding peak will fail.

How to get around this problem of random duration key press? I could think of 2 inefficient ways to do this.

During 1st iteration, take entire data & plot fft. If (no. of peaks > 1), take half data at a time plot fft. Keep on doing this until there is only 1 peak in the plot. But this technique will take a lot of time
WAVREAD the file, each key press data is separated by 0's. Thus take data between the 0's one at a time, plot fft & find peak & frequency. But this looks too simple to be used.

So is there any other effiecient method of finding the telephone number when the key press duration is random?

If you interested in computational efficiency then performing FFTs is not the best solution. Here's something you might try: apply the samples of your .wav file to a parallel bank of 8 Goertzel filters. (One filter for each of the 8 possible tone frequencies.) For example, below is the Goertzel filter for detecting the 770 Hz tone.

Monitor (measure) the output powers of the 8 filters every, say, 25 milliseconds. If the sample rate is 8000 Hz, then 25 msec means measure the 8 filter output power values every 200 input samples. So, in the figure below, apply the first 200 input samples to the filter and compute a single v(n) value. Then, and only then, compute a single complex-valued y(m) sample and single real-valued P770(m) sample. For the first 200 input samples do the same for the other 7 filters. At the end of the first 200 input samples see if two the eight Pxxx(m) power sample values exceeded some threshold value. If so, then you know a phone button key was pressed. Using some sort of logic, decode the two high-level Pxxx(m) power values to determine which key was pressed and store that single 'key number'. If no key was pushed during those first 200 input samples set the stored 'key number' value to Nan in Matlab.

Repeat the above processing for every block of 200 samples in your .wav file. If your .wav file had 16,000 samples then you should end up with a sequence of 80 'key number' values. Plot those 'key number' values versus the time index m.


Spectrogram - Am I missing something? I have the following signal:

I'm trying to compute a Spectrogram algorithm, but, don't think I'm doing it right..

I have computed the following:

1) STFT (size 256 with an overlap of 128)
 2) Computed the logs using: '10 * log10(sqrt(re * re + im + im)

This is the result that I get:

But when I use pylab in Python (for the same signal):
  x = pl.specgram(signal)

I get the following result:

Using the matplotlib I get the following:

Obviously, these are very different results.. I don't know why I'm getting these, I'm new to signal processing and spectrograms. Hope someone can help

EDIT:

This is the result I have when doing imshow in python:

This appears to be just a matter of projection. Try "imagesc" in matlab.


Pulse wave from sum of sinusoids double** table_new_pulse(ULONG length, ULONG nharms, double dcyc)
{
    ...

    ULONG harm = 1;
    double step = TWOPI / length;
    double dcyclerad = M_PI * dcyc;
    for (ULONG i = 0; i < nharms; i++)
    {
        double amp = sin(harm * dcyclerad) / harm;
        for (int j = 0; j < length; j++)
        {
            table[j] += amp * cos(step * harm * j);
        }
        harm += 1;
    }

    // normalize
    ...

    return table;
}

The preceding C function tries to approximate a pulse wave of duty-cycle dcyc by summing cosines of the first nharms each with an amplitude of sin(sin(harm * dcyc) / harm) and storing the values in a lookup table. If length is 1024 and nharms is 8, when dcyc is 50% the resulting waveform  has equal positive and negative peak amplitudes:

But as dcyc decreases from 50%, the negative peak amplitude also decreases. If dcyc is 33.3%, the negative peak amplitude is -0.540:

If dcyc is 10%, the negative peak amplitude decreases even further to -0.167, rendering the resulting waveform almost unipolar:

Can anyone explain why this happens in the simplest way possible? How can one generate a bandlimited pulse wave with varying duty-cycle and with equal positive and negative peak amplitudes?

Essentially, you are evaluating the partial sum of a Fourier Series of the rectangular periodic pulse. The pulse is symmetric around $t=0$, so the sine-parts can be omitted, leaving only the cosine-parts, as you have shown. However, you miss an essential thing: the frequency zero (i.e. the DC part), which is the overall mean of your waveform.

I dont know exactly, how you came up with the different values for amp, but it turns out the waveform with dcyc=0.5 has an amplitude of approximated 0.78. So, the mean of a waveform with a given $d_{cyc}$ is given by

$$a_0=0.78d_{cyc}+(1-d_{cyc})(-0.78)$$

You need to add this mean value to your Fourier Series to generate the equal-amplitude waveform. See below:

length = 1024
nharms = 8
dcyc = 0.2
step = 2*np.pi/length
harm = 1

A = 0.78
mean = (dcyc) * A + (1-dcyc)*(-A)

j = np.arange(4*length)

table = 0
for i in range(nharms):
    amp = np.sin(harm*dcyc*np.pi)/harm
    table = table + amp*np.cos(step*harm*j)
    harm = harm + 1
table = table + mean

plt.plot(j, table);
plt.grid(True)


Contour and area, raw (spatial) and central image moments I recently started using image moments for image processing of binary images. I read that the $0^{th}$ order contour moment is the perimeter and the $0^{th}$ order area moment is the area. These raw moments are both given by:

$M_{ij} = \sum_{x}\sum_{y}x^iy^j$.

This means If I have an image like this (but binary, foreground pixels shown in blue), the $0^{th}$ moment will correspond to the the perimeter, since it's an image of a contour :

On the other hand, if I have an image like this (foreground shown as while), I will get the area of the object as the $0^{th}$ moment:

Since I want to use the contours to get more properties, I also calculate the higher order ($1^{st}$, $2^{nd}$, $3^{rd}$ order) raw contour moment. I want to use these to get the central moments.

Formulas I am using to get the central moments are:

$\mu_{00} = M_{00}$

$\mu_{01} = 0$

$\mu_{10} = 0$

$\mu_{11} = \frac{M_{11}}{M_{00}} - x_c*y_c = \frac{M_{11}}{M_{00}} - (\frac{M_{10}}{M_{00}})*(\frac{M_{01}}{M_{00}})$

The formulas for calculating central moments are using raw moments. My question is: Which raw moments are used to calculate central moments, area or contour?. My guess is area moments, since the $0^{th}$ order central moment is also equal to the area, which is in fact the $0^{th}$ order area moment.

Additionally, can I calculate the central moments based on contour raw moments?

Actually, I was surprised how hard was it do deduct a proper definition of contour versus "normal", non-contour moments of an image. After reading a bunch of materials, here come my conclusions.

Firstly, in order to understand moments, and especially the difference and the usage of spatial (what the OP calls "raw"), central, and central normalized moments, I found two very good materials:

(manual) Johannes Kilian: "Simple Image Analysis By Moments"

Excellent manual with simple mathematics. Don't be scared by the integrals - you can read all of them as summations.

Also, it has a small overview on OpenCV functions used to operate with this moments. It's very old material (2001), so the OpenCV manual it is referring to is a bit old, but it still helps.

And than there's the wonderful third chapter, specifying which moment is used to describe which characteristic of a moment.
(image processing blog) Utkarsh: Image Moments

Simple, short and friendly. I found a lot of good material on this blog before.

Disclaimer AI Shack seemed to be offline at some point. Here is the homepage from the AI Shack author, where he talks about this project, so it still seems to be supported. I hope it comes back online soon, but if not maybe it can be tracked through the author's webpage.

Shortly, the spatial moments give information about the object in the image, i.e. related (dependent) on the object position.

The central moments are adjusted for translational invariance, by moving the origin of the "coordinate system" used for calculations to the centroid (center of gravity) of the object in question.

Finally, the central normalized moments are scaled by the area of the object, and are thus scale invariant in addition to translational invariance.

Now for the actual question part: what about contour moments?

The deductions from this part are mostly based on

OpenCV reference manual, v. 2.4.4.0, Chapter 3.5: Structural Analysis and Shape Descriptors
wikipedia page for Green's theorem - following the link from OpenCV manual, the method used when calculating contour moments (I'll explain in a few lines)
Gary Bradski, Adrian Kaehler: "Learning OpenCV: Computer Vision with the OpenCV Library"

(Link is to google books, but the relevant pages are accessible. The section that I linked to, plus the next 2-3 sections are relevant. That's about 3 pages total)

And the most important quotes from those sources:

  The moments of a contour are deﬁned in the same way but computed using the Green’s formula.

(OpenCV reference manual)

  In plane geometry, and in particular, area surveying, Green's theorem can be used to determine the area and centroid of plane figures solely by integrating over the perimeter.

(wiki for Green)

  Moreover, cvContourMoments is now just an alias for cvMoments.

(Bradski Kaehler book)

Based on that, I would deduce that contour moments does not refer to special measures of the object contours, but instead to a particular way to calculate image moments, only using the contour information (instead of pixel information for the whole image).

The difference, in the fundamental case, would be how both are calculated.

My guess would be that the direct implementation would work by pixel-by-pixel summation, directly implementing the formula. The object is expected to be filled.
My guess for the contour moments would be that the image contours are first determined (consult OpenCV manual) and then the Green theorem is applied on the contour data.

That would make the measurements slightly different for real images because the methods would differ in their: sensitivity to: noise, scaling, discretisation (pixel grid instead of continuous image). Also, the speed: calculating using contours is faster than using the direct approach. I would speculate that they would give perfectly equal results for an (idealized) continuous black and white image with no noise.

So, to answer your questions: the moments should be the same (differing because of noise etc). You can use spatial (raw) moments calculated by both methods to determine central moments (that will still describe the same thing).

Further support of this claims is the existence of this article (I only read the abstract, but should be very relevant, and even the abstract is informative) from 1994:

Luren Yang, Fritz Albregtsen: "Fast and Exact Computation of Moments Using Discrete Green's Theorem"

Note about getting the perimeter measure: I think, to get the "perimeter" which is actually just area of the contour, I would calculate the $0^{th}$ moment of the image of the contours of the objects, but treat the contours as a really thin object, instead of as "contours of an object".

All further measurements would of course differ if you used this moment further.


Getting an specific frequency component from sampled audio I have a half second audio data sampled at $44.1\textrm{ kHz}$. If I FFT it I get frequency components in $2\textrm{ Hz}$ increments. Now suppose I want the component of an arbitrary frequency, like $32.45\textrm{ Hz}$, how can I get it?

Looking at the DFT formula from Wikipedia:

$$
X_k = \sum_{n=0}^{N-1} x_n e^{-i 2 \pi k n / N}
$$

I wondered, what if I calculate the relevant $k$ and plug it into the equation? Like:

$$
{32.45 \times 22050 \over 44100} = 16.225\quad \textrm{then}\quad X_{16.225} = \sum_{n=0}^{N-1} x_n e^{-i 2\pi \times 16.225 n / N}
$$

would that get me what I want?

The commonly used Goertzel filter (or algorithm) used to detect a single frequency component is essentially a computational optimization of computing 1 bin of a DFT this way.

Note that a both a Goertzel filter and 1 bin of an un-windowed (or rectangular windowed) DFT don't just respond to a single frequency, but have a non-zero bandwidth filter response in the shape of a Sinc function, with most of the main lobe response on the order of Fs/N wide.  So the 32 Hz bin of your half-second-long FFT will contain most of a 32.45 Hz component, as well as vice versa with the Goertzel result.

Note that if the k is fractional, or not integer, which corresponds to a Goertzel filter length that is not an exact integer multiple of the period of the frequency in question, there will be some dependency of the magnitude result on the phase of the frequency component with respect to some point in the sampling window (say, the start).

And, unlike with a full FFT, 1 bin of a DFT won't tell you whether the measured component is significant compared to adjacent frequencies (is a peak), or to the rest of the signal (is above the noise floor).

But otherwise, yes.


Causal LTI system having exponential input I know that for an LTI system having complex exponential input,
i.e, $x(t)=\exp(j w_o t)$ & $h(t) \to $ LTI System ;
then, its output { $y(t) \} =M \exp(j w_o t + \phi)$ , where $M= |H(j w)|_{|w= w_o}$ & $\phi = \angle \{ H(j w)\} _{|w= w_o}$

But for a causal LTI system having exponential input,
i.e, $x(t)=\exp(a t)$ & $h(t) \to $ Causal LTI System ;
Is it true to say its output $\{ y(t)\}$ is given by:
 $y(t)  =K \exp(a t )$ , where $K= H(s)_{|s=a} \quad$  ?

Any LTI system with transfer function $H(s)$ reacts to an input signal $x(t)=e^{s_0t}$ with an output signal $y(t)=H(s_0)e^{s_0t}$, if $s_0$ is inside the region of convergence (ROC) of $H(s)$. This is equivalent to saying that the functions $e^{s_0t}$ are eigenfunctions of LTI systems, and $H(s_0)$ is the corresponding eigenvalue.

Note that $s_0$ is an arbitrary complex constant inside the ROC of $H(s)$, so the relation mentioned above is also valid for $s_0\in\mathbb{R}$. Also note that causality is not required for the above relation to hold.

However, if your question is about input signals of the form $x(t)=e^{at}u(t)$, where $u(t)$ is the unit step function, i.e., for input signals that are zero for $t<0$, then the output is generally not equal to $H(a)e^{at}u(t)$, because the functions $e^{at}u(t)$ are no eigenfunctions of LTI systems.


Understanding FFT bin resolution For a signal sampled at $f_s$, the frequency resolution (or "bin width") for an $N$ point FFT is $f_s/N$. Does this mean that the $k$th bin will contain energy from sinusoids within $0.5f_s/N$ on either side of the center frequency of that bin (ignoring spectral leakage for now)? Is this energy averaged? Is it the same as taking a $2N$ point FFT and then averaging every two bins?

The energy bin "width" of an FFT is the width of the transform of the window used (a rectangular window by default).

The transform of a rectangle the width of the FFT is a Sinc function (sin(x)/x) with a main lobe width equal to twice the DFT/FFT bin spacing, plus smaller ripples extending the full width of the FFT result (and actually circling around forever making it technically a Dirichlet function).  Note that the frequency response of each bin filter overlaps with others except at exact bin centers (representing exactly periodic-in-aperture pure sinusoids).

For the default rectangular window, only about 87% of a bin's passband is within 0.5 of a bin width from the bin center frequency, which is why ignoring so-called leakage might not be a good idea.  A longer FFT has a narrower main lobe bin response, with the ripples ("leakage") dying out faster.  But 2 half width Sinc functions do not sum to a single Sinc, so averaging a longer FFT will not produce identical results (except for sinusoids that are exactly periodic in both FFT lengths.)

A window other than rectangular (von Hann, etc.) will have a main lobe wider that 2 bin spacings, but, in return, far smaller ripples away from the main lobe.

The DFT/FFT resolution depends on whether you are trying to resolve closely spaced peaks as clearly separate peaks (requiring around double or more the bin spacing), or estimating the frequency of a single isolated peak with a relatively high local S/N ratio (sometimes in low noise allowing a fraction of a bin spacing in resolution after high quality interpolation).


Averaging multiple FFTs in matlab I have multiple, real, audio signals which I am taking the FFTs. Although they don't come from the same recording, those signals represent the same thing (biologically). I was told that I could combine those ffts and average them together to get a good representation of the averaged signal. I am thinking of using the matlab function 'pwelch' for this purpose but I have lots of questions and I could not find answers by googling.

I have already averaged the ffts themselves with the matlab function 'pwelch' to decrease the noise. When combining them together, should I just use the output of a normal FFT rather than the output of pwelch (P)?
My signals are not exactly of the same length, which will make them not usable to combine in a matrix. Should I cut them all to the same length before averaging? Or is there a way round that?
How do I actually combine them? I was thinking something like this:

:

avg = (abs(P1) + abs(P2) + abs(P3)) / 3;

and then use the function pwelch like this:

[Pavg, Favg] = welch(avg,ones(SegmentLength,1),0,NFFTavg,fs,'power');

Any insights or help would be really appreciated. Thank you!

Welch’s method is a scheme to reduce the fluctuations (noise) in spectral power estimations when you have only a single time sequence with which to work. Because you have multiple time sequences, I’m not convinced that Welch’s method is what you want. If you have time sequences $x1(n)$, $x2(n)$, and $x3(n)$ of different lengths, I suggest you zero-pad the shorter sequences so they’re equal in length to the longest of the original $x1(n)$, $x2(n)$, and $x3(n)$ sequences. (It’s my guess that no time-domain Hanning or Hamming windowing is necessary.) Perform three FFTs and compute the three corresponding spectral power sequences $P1(m)$, $P2(m)$, and $P3(m)$. Finally, average the three spectral power sequences, element-by-element, as:

Spec_Pwr_Ave = P1 + P2 + P3;

Then plot your ‘Spec_Pwr_Ave’ power-spectral sequence.


How to recover $f(t)$ from Fourier Transform of its absolute value $\mathcal{F}|f(t)|$? Let the Fourier Transform of a real signal, $f(t)$, be $\mathcal{F}(\omega)$. And the FT of the absolute value of the same signal, $|f(t)|$, be $\mathcal{F}(u)$.

Can $\mathcal{F}(w)$ be recovered from $\mathcal{F}(u)$?

For instance, the FT of $a \cdot \cos(ft)$ returns a spectrum in which the frequency $f$ has amplitude $a$.

Can $f$ and $a$ be recovered from the FT of $a \cdot \cos(ft)$?

Matt L.'s comment said it. Once you take the absolute value of the function, you lose its sign. Whatever you do to the result, like taking the Fourier transform, doesn't help. The information is lost.

It is not a part of your question, but the function could have known properties that enable recovery. If it is band-limited, like your cosine, then recovery is possible up to a sign ambiguity $f(t)$ vs. $-f(t)$. In other words, if you know the absolute value of a band-limited function at all times $t \in \mathbb{R}$, and if you know the sign of the function at time $t_1$, then you also know the sign of the function at time $t_2$. For an algorithm, see Gaurav Thakur (2010) "Reconstruction of Bandlimited Functions from Unsigned Samples".


Finding audio data in recording I am trying to find a specific block of data that should hold valid audio data that I want to process. For instance, I record at 44.1kHz, but I only transmit a .5s of content. Since the buffer I store the recordings in contains both the recording as well as background noise, is there a way to isolate the sample I desire? Is this a simple filtering problem? I know the range of frequencies I need in the sample, so would a bandpass filter solve this issue? My ultimate goal is to run the time series data through an FFT to see which of the desired frequencies were transmitted. Note that this is being done in software. Thanks for any help!

Knowing the range of frequencies is a good start and yes a bandwidth filter will certainly help you isolate them but you must consider the overall level of background noise vs signal noise.

Performing an FFT on filtered (bandwidth or otherwise) data will only yield results on the frequencies that are left over! This will either be of use to you or not.

If the number of frequencies that you are after are not that numerous, a collection of IIR (infinite impulse response) filters might be a better choice.

Consider DTMF tones used by telephones and cellphones, typically IIR filters have found their way in software based detection of these dual-tone based signals. There are online resources for calculating each of the  filter's co-efficients. You must use one filter per frequency that you are after.

Here is an iPhone version of source code for DTMF decoding using IIR filters

The file DTMFDecoder.m contains the bulk of the code with the co-efficients defined near the top.

Without knowing any further detail, this is all I am able to help you with.


DC blocking IIR filter clipping? I'm processing audio data for voice input from a mic. The data arrives in 32 bit floats [-1 ~ +1].

My first filter is to remove DC:

// x = new input value, y = filtered output value
m_x += ( 0.01 * ( x - m_x ) );
y = ( x - m_x );

When I feed it audio that is close to clipping (but not actually hitting -1 or +1), I'll get values back that actually do go above the [-1 ~ +1] limit - sometimes way above. I find this behavior curious.

Can anyone explain why this happens?

Also, what's the best way to "fix" this? Do a simple clamp for the returned y value? Pre-scale the input down via (x * 0.7071) first?

Thanks!

Your DC blocker transfer function is

$H(z) = (1-a)\frac{1-z^{-1}}{1-(1-a)z^{-1}}$

and an alternative (equivalent) difference equation is

$y_n = (1-a)(x_n - x_{n-1}+y_{n-1})$

Although this filter provides no gain it is possible to find inputs that exceed +1. For instance this input [1 -1 -1 -1 1].

Because you seem to work with a floating point realization you should not need to scale down before filtering.

Edit:
The output of your filter is bounded this way: $|y_n|<2(1-a)$. So an input gain of 0.5 guarantees that the output of the filter can't exceed +1 if the input is bounded between -1 and +1. Can you come up with an input that produces an output value that exceeds +1.5 or an input that produces an output values that is less than -1.5?

Edit2:
Filters can provide amplification as well as attenuation, otherwise they wouldn't be that interesting. Your filter does not provide amplification in the sense that a sinusoid passing through your filter will not have an increased amplitude. You can see that by plotting the amplitude/frequency response of your filter. However, for more complex inputs the non-linear phase response of your filter can cause 'overshoots' in the output. An allpass filter for instance, even though it has a completely flat amplitude response can also provide 'overshoots' for fullscale inputs.

How to deal with the spikes depends very much on what the other modules in your signal path are doing. Maybe some of the modules create headroom maybe some of them consume headroom. I don't know. Considering your DC blocker in isolation then you can apply a headroom gain of say 0.7 or 0.8 and then saturate your output. Although this setting is likely to perform some saturations my guess is that they will be completely inaudible. You will have to confirm this by experiment.


Calculating Stability of Phaser Feedback Loop I have been doing some experiments with phaser programs. I really like setting the phaser feedback to really high levels, such that the effect is really metallic and resonant, almost like a physical modeling instrument. The problem with this is that it's easy to go overboard and get runaway feedback which gets harsh extremely quickly.

I am wondering if there is a way to calculate the point at which self oscillation will occur, given the number of phaser poles, the cutoff frequency of the allpass filters, and the feedback level.

The phasor I am using uses a variable number of 1 pole allpass filters in series, each defined by the equation:

$$y[n] = A \cdot ( x[n] + y[n-1] ) - x[n-1]$$

where:

$$A = \frac{1 - \pi \cdot \frac{cutoff}{SampleRate}}{1 + \pi \cdot \frac{cutoff}{SampleRate}}$$

Here is a signal flow diagram:

Ultimately what I am looking for is an equation I can use to calculate the following:

Given a phaser with x poles with y cutoff, self oscillation will being at feedback level z

Can anyone recommend a good method or some resources I could use to figure this out?

A feedback loop like this (with enough phase diversity) will be stable if the magnitude of the loop gain is less than unity. In this specific example, it's really simple: make sure your feedback coefficient needs to be less than 1.

So it's $|z|<1$, x and y don't really matter.


What is convolution of two sine waves (tones)? Convolution of two sine waves (or tones as called in audio) is theoretically not defined as the integral is infinite. Taking finite duration windowed sine waves and doing there convolution computationally always contains a fundamental frequency equal to that of the lower frequency sine wave. I am not getting an intuitive understanding of this. Anybody who can share more understanding in this?

Convolution in the time domain is equivalent to multiplication in the frequency domain.

If you window two sinusoids in the time domain to get finite length waveforms, and the two sinusoids are exactly integer periodic in the window width, then the DFT will be impulses.  If the frequencies are different, the impulses in the two DFT will be disjoint, and the multiplication of the two spectrum will result in zero output.

If one or both of the two sinusoids are not exactly integer periodic in your window width, then the FT will result in a Sinc function.  Then the convolution will be equivalent to multiplying a Sinc function against an impulse or another Sinc function, which will produce a non-zero result.


Additional peak after resampling I prepared a Matlab script, which re-samples spectrum from 24kHz to 20kHz. When I used up-sampling, a fir filter, and then up-sampling again at the end, I get an additional peak at 4kHz which is not expected. When I checked via the resample function from Matlab paltette, there is no additional peak. What did I do wrong? What do I have to add to my script to avoid this kind of situation. Thank you in advance

clc;
close all;
clear variables;
Fs1 = 24000;
TimeStep1 = 1/Fs1;
Lenght = 10000;
time1 = (0:Lenght-1)*TimeStep1;
InputSignal = sin(2*pi*8000*time1);
NFFT1 = 2^nextpow2(Lenght);
Y1 = fft(InputSignal,NFFT1)/Lenght;
f1 = Fs1/2*linspace(0,1,NFFT1/2+1);
figure
plot(f1,2*abs(Y1(1:NFFT1/2+1)))
title('Original Spectrum')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
Fs2 = 20000;
T2 = 1/Fs2;
nGcd = gcd(Fs1, Fs2);
upY = upsample(InputSignal,Fs2/nGcd);
firCoeff = fir1(255, 0.9, 'low');
filteredY = filter(firCoeff,1,upY);
dataOut = resample(InputSignal,Fs2,ceil(Fs1));
L2 = length(dataOut);
t2 = (0:L2-1)*T2;
N=2;
lenght = (L2/N);
NFFT = 2^nextpow2(lenght);
Y = fft(dataOut,NFFT)/lenght;
f = Fs2/2*linspace(0,1,NFFT/2+1);
figure
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Resampled Spectrum')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
figure
dataOut2 = downsample(filteredY,Fs1/nGcd);
Y2 = fft(dataOut2,NFFT)/lenght;
plot(f,2*abs(Y2(1:NFFT/2+1)))
title('UpFIRDwn Spectrum')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

You are upsampling by a factor 5, which creates mirrors of the original spectrum (see the spectrum of $Yup$). With the filter you intent to remove these mirrors such that you can safely downsample the signal without aliasing.

To prevent aliasing, the upper limit for your cutoff frequency is half the lowest sampling rate, which is $0.5*Fs2 = 10 kHz$. However, the finite steepness of your filter creates a so-called transition band, which leads to some aliasing for higher frequencies for this maximum cutoff frequency. Therefore, a cutoff frequency a little below the absolute maximum is desired, for example 9 kHz. Choosing the transition band depends on your filter length and singal content.

For a cutoff frequency of 9 kHz at 5 * 24 / 2 khz, the normalized cut-off frequency is 0.15. You can use $Fs2/nGcd * Fs2$ to calculate this normalized cutoff frequency automatically.

Some guessing from my side:
In your code you use a normalized cutoff frequency of 0.9, which might be related to the 9 kHz and the output sampling rate of 20 kHz. However, the filter is applied to the $Yup$ signal, which has a sampling rate of 5 * 24 kHz.

In short: changing the normalized cutoff frequency of your filter from 0.9 to 0.15 should solve your problem.

EDIT: Follow up on using fir2

With the fir2(N, F, A) method of Matlab you can design an N points FIR filter by prescribing the filter magnitude at certain frequencies. In your case, you want to design a low pass filter, which means that the frequencies and magnitudes should be something like the following:

A = [1 1 0 0];
F = [0 wpass wstop 1];

where wpass and wstop are the relative frequencies up to where you want the filter to let the signal pass or stop, respectively.

In your case you can set these for example as follows:

wstop = Fs2/(Fs2/nGcd*Fs1);
wpass = 0.9 * wstop;

Selecting the filter length N is a tradeoff between complexity and the ripple of the filter with respect to an ideal low pass filter.

Note: In your question you have m=[0 0 1 1], which ticks on a high pass filter.


Discrete Fourier Transform: Prove that only N distinct values of {X(k)} can be computed Please can anyone answer this question:

The Discrete Fourier Transform (DFT) is a sequence of $X(k)$ in frequency domain of the  time sequence {$x(n)$} of $N$ terms. The periodic property of the DFT is often used to reduce its computational time. Prove that only $N$ distinct values of {$X(k)$} can be computed.

It's not quite clear what you are supposed to consider as starting point for your proof. You have to make some assumptions of what is known before you can prove something else. So unless you give a more detailed description of what you want to prove from what, I will give you a proof that uses a sensible set of assumptions.

Let's characterise the discrete Fourier transform as a map between two complex vector spaces $$\mathcal{F}:\mathbb{C}^N\to\mathbb{C}^M$$ of dimension $N$ and $M$ respectively. The map is supposed to preserve the structure of the vector space, so it must be linear.

A linear map between finite dimensional vector spaces can be completely described using a matrix, which in this case is $M\times N$. The rank of this matrix $F_+$ is less or equal to the smaller dimension, i.e. $\mathrm{rank}(F_+)\le\min(N,M)$.

We also require the map to be a transform, meaning that it can be inverted. So $$\mathcal{F}^{-1}:\mathbb{C}^M\to\mathbb{C}^N$$ must exist and $$\mathcal{F}^{-1}\circ\mathcal{F}=\mathrm{id}_N$$ as well as $$\mathcal{F}\circ\mathcal{F}^{-1}=\mathrm{id}_M$$
and the matrix $F_-$ of the inverse transform is also less or equal to the smaller dimension, i.e. $\mathrm{rank}(F_-)\le\min(N,M)$.

We can write the inversion conditions as a matrix equation
$$F^{-1}\cdot F = E_N$$
$$F\cdot F^{-1} = E_M$$
where $E$ is the identity matrix of the size specified in the subscript. The matrix multiplication can be translated to a rank inequality:

$$\min\left(\mathrm{rank}(F_-),\mathrm{rank}(F_+)\right)\geq\mathrm{rank}(E_N)=N$$
and
$$\min\left(\mathrm{rank}(F_+),\mathrm{rank}(F_-)\right)\geq\mathrm{rank}(E_M)=M$$
because the product of two matrices has at most the rank of the matrix with the lesser rank.
From above we also had
$$\mathrm{rank}(F_-)\le\min(N,M)$$
and
$$\mathrm{rank}(F_+)\le\min(N,M)$$

These four inequalities only allow for one solution, which is $M=N$. To see this, assume that for example $N<M$, then it follows that $\mathrm{rank}(F_-)=\mathrm{rank}(F_+)\leq N$ but also $\mathrm{rank}(F_+)\geq M$, from which $N\geq M$ would follow, making the assumption inconsistent. The same goes for assuming $N>M$.

That means the discrete Fourier transform must be between complex vector spaces of equal dimensions. Or in other words, the number of frequency bins must be equal to the number of time samples.

The proof above can of course be a single line if you already know that vector space isomorphisms require matching dimensions. Of course, some linear algebra is still necessary, but the rank formulas should be rather intuitive.


UWB Indoor Positioning I am trying to understand the UWB technology and how it can be used for localization.

Are there any resources where I can study about the algorithms it uses to calculate the direction?

Also why does it able to provide a centimetre accuracy?

For my understanding, it sends pulses in different frequencies within the frequency band. I am not quite sure how the communication is established. Is there a specific protocol used similar to Bluetooth?

Thanks a lot!!!

Are there any resources where I can study about the algorithms it uses to calculate the direction?

"positioning" and direction finding are typically two different things; you can use direction finding with enough observers for positioning, but it's not the only way, and probably not used here. UWB and angular resolution is hard, because the usual "beautiful" algorithms for direction of arrival assume something narrowband.
There's nothing inherently fundamentally different about UWB pulse radar and pulse radar in general. And there's about 116 years of pulse radar history; I'm going to go ahead and assume you're a student; maybe ask your supervisor or professor in a radar or high-frequency lecture what they can recommend; it's probable they have a good idea of what would embed well into your educational background.

Also why does it able to provide a centimetre accuracy?

Radar systems are limited by noise and by bandwidth. Remember what "UWB" stands for!

For my understanding, it sends pulses in different frequencies within the frequency band.

A pulse is a wideband thing, it doesn't have a single frequency; so there's something a bit questionable about that wording.
Now, what could very well be is that it's a stepped-frequency wideband radar.
However, we don't really know which kind of UWB radar you're referring to, so I'll refrain from speculating here – might only confuse the both of us!

I am not quite sure how the communication is established.

Localization or radar are not communication methods.

Is there a specific protocol used similar to Bluetooth?

Almost certainly not.


Why does FEC seem to make BER worse at low Eb/N0? I've been plotting BER curves and am wondering why FEC seems to make a BER curve worse at low Eb/N0 than the un-coded theoretical BER curve?

My example plot below is for Differential BPSK using a Rate 1/2 k=7 convolutional code on the AWGN channel.

This is a behaviour that is commonly outlined in some textbooks and tutorials on FEC, but usually in the form of an observation. For instance, Turbo Codes are sometimes characterised by their Threshold $E_b/N_o$, which can be considered as the minimum SNR per bit from where the BER starts to decrease. I have not seen a formal explanation for this property, but it can be deduced at least for some types of codes I assume. In this Wikipedia article, for instance, it is explained as a consequence of the risk-polling effect that comes from the noise averaging capability of FEC. Let me complement this, also by an intuitive approach.

Consider for instance the use of a $(n,k)$ binary block code with hard decision decoding. The maximum number of correctable bit errors for this code is commonly denoted by $t$.
If the number of bit errors within each transmitted block is higher than $t$, then the receiver will not be able to recover the original codeword. Moreover, the received block will be at roughly the same distance of multiple codewords within the set of size $2^k$. The decoder then picks one of those codewords with high uncertainty, which results in a BER close to 0.5.

Simply put, a few bit errors - though above the error correction capability $t$ of the code - may result in a complete block error, which explains why FEC underperforms the uncoded system for a certain SNR range.

This is why in systems that operate at low SNR or when errors appear in bursts, it is a good idea to use interleaving in order to spread the bit errors along several blocks and keep the number of bit errors per block below $t$.


How to remove certain background from image? This was the image obtained after Gabor filtering... Is there some way of removing the lines inside the image except the bright white rounded defect.

Did try another orientation of gabor filter go the following result:

After thresholding :Which is ok enough but not great

After Median filtering:

Thanks in advance

If the output of your Gábor-filter is reliable and the variation in your image-data is not so high that the outcome looks completely different, then you can use the following approach (parts of it where already mentioned):

Binarize your second image with any automatic thresholding algorithm. The range of  thresholds which will work is large as you will see.
Use a component labeling algorithm to label each connected region of pixel with a unique value.
Calculate for every component of your image a property which describes, how close your object is to a filled circle. For this, you can use for instance the compactness. I used below the fraction of pixels within the equivalent disk radius. This radius is the radius a disk with the same area as your object would have.

Area is simply counting the number of pixels of your object
Calculating the radius of a circle given its area should be manageable too
For the pixel inside this radius, you need the center of each image object, but this is just the mean of all positions of the object-pixels. If you are a physicist then you know the center of mass of several point-masses. This is equivalent.
Now you calculate for every pixel of each object the distance from its center and check whether or not it is smaller then the circle radius. Dividing the two numbers and you get the fraction of inside and outside pixel.

Take the object with the highest fraction. Note: objects with only one pixel will get a value of 1. So you should set a size-threshold and take only objects larger, say 10 pixel.

To see, that an automatic binarization should work, here are the results for a very low and a very high threshold:

Update Component Labeling

The choice of the component labeling is not critical for your application. I would suggest, if you have to implement it by yourself, to use a very simple method. The one pass version from the Wikipedia site is very easy. Basically you iterate through your binary image and when you meet a pixel which is white and is not already labeled you use a new label to for this object and start with this pixel.

The process of labeling this object with a label is basically similar to a flood fill. This is on the Wikipedia-site the inner steps 1-4 in the algorithm. You start with this labeled pixel and put all its neighbors on a stack (they used a vector). For a pixel on the stack you check whether it is foreground and not already labeled. If you have to label it, you put again all its neighbors on the stack. Do this until your stack is empty.

Then you continue your scan through the image. Unlike the description on the Wiki-site, you don't have to remove a pixel from your original image, you just skip when you have a value different from 0 in your label-image.


Experimental determination of impulse response function I am doing an experiment (I major in physics, not familiar with signal processing so there may be some naive mistakes, please be kind) which requires acting arbitrary pulses on a system.

I use an AWG to generate the pulse (current), bring it to the linear amplifier, and eventually, to a coil (than the current is converted to the magnetic field and act on the system). The pulse will be distorted by the coil and the amplifier, and U want to know the impulse response function (f), so as to modify my input pulse.

Is there a method which can determine the impulse response function experimentally?

I know a little about the method which says: generate a number of input pulses, use an oscilloscope to get the distorted pulse, and f is treated as a matrix which can be determined by multiply the "inverse of input pulse matrix" and the output pulse matrix. Can this really work? Is there a standard method to determine f (I suppose there is)?

You can determine the impulse response of a system using Gaussian random noise as a test signal, let's name it $x$. Then, you have to measure the system response $y$ to that noise. With a considerable amount of digital samples from both signals then you can compute the crosscorrelation between $x$ and $y$ as follows:

$$
R_{xy}(n\Delta t) = \frac{1}{N-n} \sum_{k=1}^{N-n-1} x[k] \cdot y(k-n\Delta t)
$$

If the white noise has a bandwidth that is much greater than the bandwidth of the system (for instance, 100 times) then the crosscorrelation is equal to the impulse response times the power spectrum density of the noise (which is constant) divided by 2. Of course, there is some uncertainty in the estimation, which decreases with large values for $N$.

Source: Oliver-Cage: Electronic Measurements And Instrumentation


Is it possible to lose half the samples, and yet see the full spectrogram as the original signal? Consider an speech audio signal sampled at 16,000 samples per second. Plot its spctrogram.

If we delete every other sample from input signal, we get y[n] = x[2n] for all n.

It seems that if we plot the spectrogram of y[n], we would only see frequencies upto 4kHz on the Y axis (all frequencies above 4k would be zero). So this scheme of resampling is not useful.

Is it possible to do some operation on x[n], and get a signal (y[n]) with only 8000 samples per second, and yet see the complete frequency range 0Hz - 8kHz?

Of course, this will not be as accurate as the original spectrogram, but still I'm interested whether the information in original "4kHz - 8 kHz" frequency range is still somehow partially preserved after some resampling operation.

If yes, what's the mathematical formula for this new resampling scheme?

To rephrase your question: Is there a process that halves the number of samples of a full-bandwidth signal without significantly affecting the appearance of the spectrogram, except for scaling of time or frequency axes?

Resampling does not do that. The only such process I know of is time stretching / pitch scaling followed by decimation. Actually, one definition of time stretching / pitch scaling would be to "modify the signal so that the spectrogram is stretched".

Time stretching and pitch scaling are equivalent from the point of view that the number of samples needed to represent the processed signal changes. Time stretching to half the original length requires only half the number of samples, and pitch scaling to half the frequency allows to discard every second sample without causing aliasing, assuming that the empty upper part of the frequency spectrum is not filled by something made-up by the algorithm used.

Time stretching / pitch scaling does not have a unique definition or a unique mathematical formula. Well-done time stretching should match one's expectation of, in your case, "speeding up" the process that generated the audio. Speeding up natural processes is not uniquely defined. For example, should the pianist be made to move their fingers faster, resulting in changes in note attacks in the least, or should we be more agnostic of the actual process taking place and just isolate the notes in the audio and move them around into a more compact pattern? And should the note decay times be made shorter, too, because we want to make everything shorter? Depends on one's expectations or needs.

There are various algorithms though. Here are before and after spectrograms of the phrase "hello" (freesound), processed by rubberband -t0.5 speech16kHz.wav out.wav with Rubber Band v1.8.2 command line and analyzed in Adobe Audition 3.0 with 512 bands and Blackman–Harris window:

Figure 1. Top: original speech spectrogram, bottom: spectrogram of the speech time-stretched to half the original length and half the number of samples. The spectrograms are zoomed out to full view for easier comparison.

The spectrogram of the time-stretched signal looks essentially the same, but has reduced detail and resolution.


How to amplify sound wave I'm developing a singing game and I want to output player's voice through speakers. I'm using portaudio for that, I acquire the input signal and just write it to the output buffer. The problem is the voice is bare audible. I tried to write to output buffer multiplying by a constant but that distorted the voice.

How would I do that?

Thanks.

The first and best solution is to increase the input gain which should be possible on most platforms. If not, then read on...

If multiplying by a constant caused distortion then either:

it is too loud and distortion is inevitable (barring what I'll just call "fancy" signal processing)
you are doing something wrong.

I'll assume #2 since you started with the case of "barely audible", and you should learn how to apply gain before doing anything "fancy" anyway. To apply gain, make sure you multiply each sample individually by the same constant. Remember that the sample count = frames * channels. Loop through each each sample in the buffer and multiply it by the gain.

For this case, you probably want to use the float datatype in PortAudio, but any type will work.


When do we normalize data in audio processing? Lets assume we have a simple audio processing flow with preprocessing, feature-extraction and classification.

Where and why do we normalize the data?

I know the purpose of normalizing the data in the preprocession-step. But as I implemented some feature-extracrion methods like spectral flux, I have read about normalization of the extracted features.

Do we need to normalize the extracted features, when we normalized the data in the preprocessing step and vice versa?

There is no definite answer to that question since different workflows are applied to different tasks. For example in some applications you might want to normalize your input data samples or perform some kind of Automatic Gain Control (i.e. you work in fixed point and even for some low signal values you want your FFT to have good resolution). On the contrary you might be using some energy features, such as very simple energy based Activity Detection, thus normalisation of your input samples will mess your results. Therefore you must use your common sense and think (or test) whether your features are dependent on input amplitude.

Other thing that you probably had in mind, is to normalize the features itself. This is called a whitening and is used for example in speech recognition for removal of channel effects. It is done by removing mean of your features and (if you wish) normalizing by the variance. This make your features to be in the same space.

To wrap up, it's up to you. Train/test both systems with and without any type of normalization and observe the performance.


2D version of Shift Theorem What is the 2D version of Shift Theorem?

The 2D version is just the same as the 1D version:
$$
x(n,m) \leftrightarrow X(k,l) = \sum_{n=0}^{N-1} \sum_{m=0}^{M-1} x(n,m) e^{-j2\pi nk/N - j2\pi ml/M}\\
x(n-a,m-b) \leftrightarrow
\begin{array}{ccl}
&&\sum_{n=0}^{N-1} \sum_{m=0}^{M-1} x(n-a,m-b) e^{-j2\pi nk/N - j2\pi ml/M} \\
&=& \sum_{n'=0}^{N-1} \sum_{m'=0}^{M-1} x(n',m') e^{-j2\pi (n'+a)k/N - j2\pi (m'+b)l/M}\\
 &=& X(k,l) e^{-j (2\pi ka/N + 2\pi lb/M)}
\end{array}
$$


Quantization level versus bandwidth Why is it that increasing the number of quantization levels $L$ results in an increase in the signal's bandwidth?

Say you quantize a signal with $L=2^n$ levels. Then, every quantized sample is represented by $n$ bits. So, as you increase $L$, you need more and more bits to represent the samples. This means that transmitting the signal will require transmitting more bits.

Also, keep in mind that, in order to convey bits at a rate of $R_b$ bits per second, you need a signal with bandwidth at least $B=R_b/2$ hertz (assuming binary communications).

Now (this is important) assume that you want to transmit the signal in real time, like in a voice application, or in streaming. If you increase $L$, you will have to transmit more bits per second, which means a larger $R_b$, which in turn implies a larger $B$.

As an example, say you're transmitting telephone-quality voice with $L=256$, $n=8$, and with 8000 samples per second. Then, $R_b = 8 \times 8000 = 64000$ bits per second, which require a signal of bandwidth $B > 32000\text{ Hz}$.

Say you want to improve the signal quality and go with $L=1024$ ($n=10$). Now, you need $R_b = 80000$ and $B>40000\text{ Hz}$. An increase in $L$ produced an increase in $B$.

Note that, if there is no real-time requirement, you can instead keep a signal with the same bandwidth (so, same bit rate) and instead take longer to transmit. So, in general, increasing $L$ will cause an increase in bandwidth and/or an increase in transmission time.


OFDM demodulation, how to set the sampling rate? I have a signal for which I know the OFDM burst length (from the autocorrelation function of the signal) T_burst and the FFT size 1024 (by measuring the channel spacing) and the guard interval (128).

How should sample such a signal for OFDM demodulation ?
What is the ideal sample rate for that ?
Is it (1024+128)/T_burst?

I plan to use GNU radio for demodulation, like this (this image is from IRC channel - credit to drmpeg):

Your receiver will consist of, apart from synchronization blocks, an FFT. You want that FFT to generate one output bin for each of the subcarriers. Therefore, you should make use of the relationship between the input sample rate of a DFT to its output bin spacing:

$$
\Delta f = \frac{f_s}{N_{FFT}}
$$

where $\Delta f$ is the subcarrier spacing and $f_s$ is the sample rate (both in Hz). Assuming you know the subcarrier spacing and the FFT size, then it's straightforward to choose the sample rate:

$$
f_s = N_{FFT} \Delta f
$$


Howto calculate SNR for EEG data? This question is about SNR in the context of EEG. (a related question  about SNR)

I am interested in calculating The SNR of an ERP. My motivation is: To calculate the "signal-to-noise ratio" and from this to extract the information (how?) on how many trials are needed to construct a stable ERP.

my Main question is: how do I calculate and SNR for such scenario  ?

some more explation about EEG and ERP are enclosed:
EEG:

  The on-going electrical activity of the brain measured from scalp electrodes is called the electroencephalogram or EEG.  within the EEG is a signal which is more revealing about information processing in the brain. This signal can be obtained by time-locking the recording of the EEG to the onset of events such as a person reading a word on a computer screen, listening to a musical note played on an instrument, or viewing a picture in a magazine. The resulting activity is called an "event-related potential" (ERP)

SNR:

  The number of trials necessary to obtain an ERP depends on a number of factors, the most important being the "signal-to-noise ratio", that is, the relative size of the signal (the ERP) relative to the size of the noise (the background EEG). In cognitive experiments 30 to 50 stimulus presentations are typically required to obtain a good clean average ERP.

from: http://neurocog.psy.tufts.edu/images/erp.htm

First you have to clearly define what you mean by "noise". Is simply the electronic noise in your measurement system or is it also brain activity that's not related to your stimulus.

One easy way to determine SNR is to compare the signal level during a stimulus with the signal level without stimulus. If noise and signal are uncorrelated you can determine the signal level through energy subtraction and thus determine the SNR.

Depending on your noise definition, this may not be practical cause you can't easily shut off the brain of the subject :-).


FFT of Square Wave and Sin Wave When we have FFT of Square wave of a pure real signal with high frequency(close to f(max), less than half of sampling frequency) we see only 2 symmetric high peaks but when we decrease the frequency we see the two peaks and peaks of harmonics and they are getting smaller and smaller, why? Also the highest peaks are getting away from each other, movig <- and -> , vice versa -> <- when we increase the frequency. why?

FFT of Sin wave has only 2 peaks no matter what and they also get away from each other when the frequency is low and get more and more closer when the frequency gets higher. Why? I would be glad if you can enlighten me!

Lets say our sample count is 64 and sampling frequency is 8kHz and signal frequency 2kHz. We should see 2 peaks one at 16 and one at 48.

The picture above is only the left half of my graph, there is also a symmetrical part on the right. My signal is pure real. These small harmonics dont appear at my above settings but when I decrease the freuency to 600 Hz for example I see harmonics like this for square wave and none for Sine Wave. So on high frequencies they disappear and low frequencies they appear, the lower the more of these harmonics. Also the peaks are  getting away from each other when I lower the frequency and getting closer when I increase the frequency. At 4 kHz I expect to see only 1 Peak, two of them will unite. Why? Also I couldnt test what happens if I run with 5kHz :( and when I run the code with lots of samples like 2048, the graph displayed makes no sense at all. Why?

I know I'm asking a lot questions but I would be glad if you have some time and help me!

The odd harmonics are what makes a square wave different from a sine wave in the frequency domain.  Otherwise the FFT would lose information about the wave shape of the time domain input.

The higher the fundamental frequency, the farther apart are any multiples of that frequency, such as overtones or harmonics.