Understanding spectra in polyphase filter banks I'm trying to fully understand the spectral effects in a polyphase channelizer / polyphase analysis filter bank. I don't really have a signal processing background, so please excuse some incorrect terminology ;-).

So let's consider a very simple case with two output channels. I'll be referring to Polyphase Channelizer Demystified
[Lecture Notes] in this question. Let's start with this block diagram:

and this input spectrum:

I want to understand what spectra at A1, A2, B1 and B2 look like. So obviously as the input signal is sampled we have a repeating spectrum somewhat depicted in the graph above (showing one repeated half of the grey channel on each side) the paper explains downsampling by factor 2 will scale the signals by factor 1/2 and linearly add one shifted version.

At A1 we have such a simple down-sampling operation, so we should get this:

the signal is downsampled by factor 2, so we have a new sample frequency. Still, the linearly combined signal is plausible, as the spectrum repeats according to the new output sample rate.

Equation (6) and (7) in the paper describe the spectrum of a downsampled signal with one unit delay. This should be the signal B1 and it should look like this:

However, this is not plausible, as the signal spectrum does not repeat properly according to the new sample frequency. What did I miss here?

The next question is: What is the effect of the downsampled filter phases? The filter-downsampling will lead to a 'repeated spectrum' filter transfer function as well, correct? So is this correct for A2?:

And what's the effect of the offset/shift of the filter coefficients in p1 caused by the polyphase partition?

Now if we only add up the Nyquist band of the A1 and B1 pictures we end up with the orange spectrum representing channel 0. If we do the same but with a phase shift A1 and -B1 gives the grey spectrum, channel 1. So the end result is plausible, but the intermediate signal at B1 is not.

EDIT:
To add some equations to this topic: According to R. E. Crochiere and L. R. Rabiner. Multirate Digital Signal Processing, pp. 84–88

Filter partitions:
$p_\rho(n)  =  h(nM + \rho)$

Equation for polyphase branches:
$x_\rho(n)  =  x(nM - \rho)$

Since you define the input to the filter $p_1[n]$ as $x_1[n]=x[2n-1]$ and the channel impulse response as $p_1[n] = p[2n+1]$ the spectrum at the output of that channel will be as follows:

First express $X_1(e^{j\omega})$, the DTFT of the input to the channel no 1.

$$ x[n] \longrightarrow \boxed{z^{-1}} \longrightarrow v[n]=x[n-1] \longrightarrow \boxed{\downarrow2} \longrightarrow w[n]=x[2n-1] $$

Now as can be seen, $V(e^{j\omega}) = X(e^{j\omega}) e^{-j\omega} $ and also $W(e^{j\omega}) = \frac{1}{2} \sum_{k=0}^{1} V( e^{j( \frac{\omega-2\pi k}{2} )  })$ When combined it yields;

$$ X_1(e^{j\omega}) = \frac{1}{2} \sum_{k=0}^{1} X( e^{j( \frac{\omega-2\pi k}{2} )  }) e^{-j( \frac{\omega-2\pi k}{2} )} $$

In this case

$$ X_1(e^{j\omega}) = \frac{1}{2} e^{-j\frac{\omega}{2}} \cdot(X(e^{j\frac{\omega}{2}}) - X(e^{j(\frac{\omega}{2} - \pi)})) $$

Note: This is where the linked paper an therefore image four in the question is confusing. It does not include the linear phase part $e^{-j\frac{\omega}{2}}$ in equation (7). Without this linear phase part, as mentioned in the question, the signal is not a valid factor two down sampled signal. The phase term 'fixes' the equation to make the phase 'repeat' properly according to the new sample frequency.

That's the signal spectrum, and the filter spectrum similarly is:

$$ h[n] \longrightarrow \boxed{z^{1}} \longrightarrow v[n]=h[n+1] \longrightarrow \boxed{\downarrow2} \longrightarrow w[n]=h[2n+1] $$

$$ H_1(e^{j\omega}) = \frac{1}{2} \sum_{k=0}^{1} H( e^{j( \frac{\omega-2\pi k}{2} )  }) e^{j( \frac{\omega-2\pi k}{2} )} $$

In this case $$ H_1(e^{j\omega}) = \frac{1}{2} e^{j\frac{\omega}{2}} \cdot(H(e^{j\frac{\omega}{2}}) - H(e^{j(\frac{\omega}{2} - \pi)})) $$

Note: The linear phase term has opposite sign to the phase term in $X_1$ so these terms will cancel each other when we multiply $X_1$ and $H_1$.

Hence the spectrum at $B_2$ is:

$$B_2(e^{j\omega}) = H_1(e^{j\omega}) X_1(e^{j\omega}) $$
$$ B_2(e^{j\omega}) = \left( \frac{1}{2} \sum_{k=0}^{1} H( e^{j( \frac{\omega-2\pi k}{2} )  }) e^{j( \frac{\omega-2\pi k}{2} )} \right)  \left( \frac{1}{2} \sum_{k=0}^{1} X( e^{j( \frac{\omega-2\pi k}{2} )  }) e^{-j( \frac{\omega-2\pi k}{2} )} \right) $$

in this case:

$$ B_2(e^{j\omega}) = \frac{1}{4} \cdot (X(e^{j\frac{\omega}{2}})H(e^{j\frac{\omega}{2}}) - X(e^{j\frac{\omega}{2}})H(e^{j(\frac{\omega}{2} - \pi)}) -X(e^{j(\frac{\omega}{2} - \pi)})H(e^{j\frac{\omega}{2}}) + X(e^{j(\frac{\omega}{2} - \pi)})H(e^{j(\frac{\omega}{2} - \pi)})) $$

For $A_2$:

$$ A_2(e^{j\omega}) = \frac{1}{4} \cdot (X(e^{j\frac{\omega}{2}})H(e^{j\frac{\omega}{2}}) + X(e^{j\frac{\omega}{2}})H(e^{j(\frac{\omega}{2} - \pi)}) +X(e^{j(\frac{\omega}{2} - \pi)})H(e^{j\frac{\omega}{2}}) + X(e^{j(\frac{\omega}{2} - \pi)})H(e^{j(\frac{\omega}{2} - \pi)})) $$

Note that the aliasing cancellation happens at the reconstruction part of the system: As can now be seen $A_2 + B_2$ results in one channel, whereas $A_2 - B_2$ will produce the other channel.


Limitation on the shift theorem of DFT due to frequency resolution? The shift theorem states that shifting a sine wave in time domain by t is equivalent to multiplying the corresponding DFT coefficient of the signal by a complex exponential e^(-jwt). Described by the following operation:

However, I realized that in order to apply the time shift precisely, the frequency of interest has to be multiples of the frequency bin size.

Here is a quick example I did in MATLAB - I simulated two signals (blue & red), red being the delayed version of the blue (with arbitrary Fs and signal length). My goal is to time shift the red signal by the delayed amount of time so that it overlaps with the blue (forget about the discontinuity caused by the circular shift property for now!!) :

In this image, the frequency of my signal (104.65Hz) is the exact multiples of my frequency bin size, therefore the shifted signal matches the original perfectly.

In this cases, the frequency of my signal (150Hz) lies in between my frequency bins. Despite the fact that 150Hz is very close to 104.65Hz, the shift is off from my desired delay amount.

Now imagine if I have a slightly more complicated signal with multiple sinusoidal components, and I would like to shift each of them by a different t. It is virtually impossible to choose my frequency bin such that it can be divided nicely by all these frequencies. What should I do in that case?

I guess I was wondering:

If someone can explain to me where exactly does this kind of shift offset come from.
If I have a more complicated signal described above and I want to shift different sinusoids by different amount, how can I get around this problem?

Many thanks!

But what if I want to shift different frequency components by a different number of samples? In that case, you can't find FFT of the whole data set and multiply by the same exponential. You would have to treat different frequencies separately right?

What you're describing is a system with a non-constant (i.e. frequency-dependent) group delay.

Group delay is the derivative of phase, hence,  $\frac{\mathrm d\phantom f}{\mathrm df}\phi(f)$ of a system.

In filter design, you (most of the time) want a constant group delay, so you design your filter with linear phase (which means you need to use FIR filters, and you need to make them symmetric); the derivative of a linear function is a constant!

Now, you want a given "group delay profile". Sooo, what about actually implementing the following idea:

Essentially, you approach this via the window design method that one typically uses to design linear-phase filters of a wanted amplitude response. You, instead, want a to define a group delay, and hence, a phase response over frequency.

Convert a definition of your group delay to a phase definition by numerically integrating that profile (i.e. do a cumulative sum). (the integration is the inverse operation to the differentiation that delivers the group delay from the phase response; so, the result should be that each value is a $e^{j2\varphi}$, with $2\pi\varphi$ being the relative delay you want). Do that definition at a high "oversampling", i.e. if you just need an actual "resolution" for your frequency-dependent delay of sampling rate / 10, define that function at 200 points. Note down that you have a "frequency-sampling-rate" of $r_\text{freq}=\frac{f_\text{sample}}{200}$.
Convolve the result with a low-pass linear filter impulse response. Why will become clearer in the next step.
Convert the result to time domain by doing an IDFT (IFFT). The effect of the low-pass filter above is now, a lot of the resulting values are (close enough to) zero. So, you can omit these.

If I'm not mistaken, you now have a filter of desired group delay profile plus a constant (the constant stems from the low-pass filtering in frequency domain).

Another route if you only care about non-overlapping bands and these should have very different group delays, making it hard to properly build that filter via the filter design method:

Find the "hardest" (i.e. narrowest) bands; design a FIR band pass filter for these using your favorite math toolkit's filter design methodology (I propose designing an equivalent-width low pass and shifting it in frequency). The length of that filter will be the longest you'll need for all the band pass filters. Sufficiently over-sample, i.e. re-design that filter, but use let's say 20 times the length. Use that length for all other filters, too.
The group delay of a linear-phase phase filter is half its length. Now, all your bandpasses have the same group delay $\tau_\text{offset}$.
Adjust the group delay of the individual filters by transforming (FFT), multiplying by a complex sinusoid of appropriate frequency, and transforming back (IFFT). If your frequency isn't (close enough) to $f_s / \text{length}(filter)$, you'll get leakage / Gibb's phenomenon. You can trade accuracy for oscillations by applying a window. If your filters are long, it's probably not even necessary to window.
Sum up all the filter impulse responses.
(optional) if it turns out that a lot of the resulting taps are close to zero, truncate.


Frequency domain filtering: must filter coefficients always be conjugate symetric? While trying to do frequency domain filtering of audio signals using windowed overlap-add methodology, I came across some pre-calculated frequency domain coefficients. I noticed that the coefficients, complex numbers, are not conjugate symetric. Can that be? What does it tell me about the filter?
If the input filter is FIR, described by y(n)=b(1)x(n)+b(2)x(n−1)+⋯+b(nb+1)x(n−nb), then the discrete Fourier transform of this real filter coefficients should be conjugate symetric as for any other real input. Furthermore, after the input signal is transformed to frequency domain, multiplied with the filter coefficients, and transformed back to time domain using inverse DFT, the signal (audio_filtered) should again be real?
A = fft(audio);
H = fft(filter);
Afiltered = A .* H;
audio_filtered = ifft(Afiltered)

If the input filter is IIR, are things different with respect to conjugate symetry?

I noticed that the coefficients, complex numbers, are not conjugate symetric.

Check what is the frequency associated with each coefficient, usually the coefficients are defined such that $a_i$ corresponds to the frequency $f_s i / N$, where $f_s$ is the sampling frequency, and $N$ is the number of coefficients. So what you have is $a_i = conj(a_{N-i})$, in some cases they may be circularly shifted, and so on.

the signal (audio_filtered) should again be real?

Yes, if you think in the time domain, the convolution involves only real values, so the result is a real signal. If you think in the frequency domain, both the filter and the signal have Hermitian symmetry, so the product of them will have hermitian symmetry as well.
Proof, let the signal with FFT $X$, and the filter with FFT $H$, the FFT of the filtered signal $Y$ is
$Y_{i} = X_i H_i = conj(X_{N-i})conj(H_{N-i}) = conj(X_{N-i}H_{N-i}) = Y_{N-i}

If the input filter is IIR, are things different with respect to conjugate symetry?

No, if you look at the fourier transform of an infinite signal
$$X(\omega) = \sum_{k=0}^{\infty} u_k exp(-2j\pi k \omega)$$
Assume $u_k$ is real, and will will notice that since $e^{-j x}$ is the complex conjugate of $e^{+jx}$, each term will be conjugated and so, $X(\omega)$ is the complex conjugate of $X(-\omega)$


Convolution in Image Processing Does anyone know how to write arg min code in matlab with explanation?

I am referring to a paper for my research. I am not sure whether I can paste the paper link here. At first, I try to do feature extraction of finger using Gabor Filter. According to the paper, after feature extraction level, I need to do the following using formula to get my desired image:

OriCode(x, y) = arg j min {I(x, y) ∗ G(x, y, θj)}

∗ stands for convolution,
I(x,y) stands for input image,
G(x,y,θj) stands for gabor filtered image

There is actually a small j between arg and min. OriCode ranges from 0 to J-1
whereby j ranges from 0 to J. Here,J is set to 6.

I try to use conv2 to convolute my original input image and gabor filtered image. End up I am getting a blank white image. And what is the arg min used for? Please help me. Thanks.

I'm going to use a grayscale image for simplicity. Suppose you have a 3D matrix and each slice is the result of filtering the image with the Gabor filter for that angle.

I = randn(256, 256); % Just some random image
G = gabor_filter; % Assume G is 256x256xNum_Angles, e.g. 256x256x6 if you had 6 angles
Y = zeros(256, 256, size(G, 3)); % Preallocate a temporary matrix

for idx = 1:size(G, 3)
    Y(:, :, idx) = conv2(I, G(:, :, idx), 'same');
end

Ymin = min(Y, [], 3); % This takes the minimum along the 3rd dimension


Spectral leakage from mathematical point of view Does anyone can explain or propose a reference that mathematically and quantitively investigates the effect of leakage on the magnitude and phase of a signal?

The reference is simple convolution.  So-called spectral leakage is due to multiplication with a window function.  Feed any longer signal to a finite length DFT or FFT, and a default rectangular window (if not some other shape of window) is usually applied to the longer signal.  That multiplication in the one domain results in (circular) convolution of the transforms in the other domain.  The results of this convolution are sometimes called "leakage".


How does averaging increase the Signal to Noise ratio? When attempting to spot activities on an EEG waveform given a stimuli, it is recommended that several trials be taken and their average is done to increase the SNR.

Intuitively, since noise is not stationary or localized, averaging different trials will have a smearing effect on the noise.

But how exactly does it increase the SNR since the noise is still present (but averaged out)? Can we more rigorously define the effect of averaging on the SNR?

Intuitively this is true, because averaging a zero mean noise processes approximates its expectation value - which is zero. More rigorously:

If the signal $x$ that you want to observe (estimate, actually) is constant for all observations $y$ we can write the $n\text{th}$ observation as
$$
y_n = x + r_n
$$
where $r_n$ is the noise which is different for every observation, of course. The SNR $\gamma_1$ of $y_n$ is
$$
\gamma_1=\frac{\operatorname{E}[|x|^2]}{\operatorname{E}[|r_n|^2]}
$$
where $\operatorname{E}$ is the expectation operator. By averaging $N$ observations of $y_n$ we get
$$
\hat x = \frac{1}{N}\sum_{n=1}^N( x+r_n) = x + \frac{1}{N}\sum_{n=1}^N r_n
$$
The SNR $\gamma_2$ of the averaged observations is
$$
\gamma_2=\frac{\operatorname{E}[|x|^2]}{\frac{1}{N^2}\operatorname{E}[|\sum_{n=1}^N r_n|^2]}
$$
To simplify that expression, we assume that the noise sequences $r_n$ are uncorrelated (1). In that case, the expectation of the squared sum of $r_n$ is equal to the sum of the expectations of the squared $r_n$, i.e.
$$
\operatorname{E}\left[|\sum_{n=1}^N r_n|^2\right] = \sum_{n=1}^N \operatorname{E}\left[|r_n|^2\right]=N\operatorname{E}\left[|r_n|^2\right]
$$
Inserting this in the expression for $\gamma_2$ yields
$$
\gamma_2=N \frac{\operatorname{E}[|x|^2]}{\operatorname{E}[|r_n|^2]}
$$
Comparing the SNR $\gamma_1$ of the single observation and the SNR $\gamma_2$ of the averaged observation we see that averaging increases the SNR by a factor of $N$.

(1) This holds for a Gaussian noise process, for example. Right now I'm not sure for what classes of stochastic processes it holds. Comments on this are welcome.


Examples of linear/logistic regression applications I was wondering about super simple applications of linear or logistic regression to signal processing (image processing specifically)?

I'm a math guy that's been asked to explain some mathematical methods to signal/image processing people, and I'd like to be able to use a couple simple examples that are related to their field, but I know little about image processing. So, I'm not looking for anything complicated (I don't know the subject, so the simpler the better; I'm just hoping to find something to demonstrate the applicability of these algorithms).

Thank you for your help!

In image processing all I can think of is using linear regression to fit a plane to the intensity surface of an image, and using logistic regression for some  kind of supervised thresholding...

You might be better off using an example from computer vision instead. For example texture classification or object recognition.


Understanding "voltage" and "power" in SNR of a sampled signal I'm facing a gap in my understanding of when a signal is a "voltage" signal and when it is a "power" signal that I've always managed to avoid resolving until now...  First what I think I understand, then my question after the horizontal rule:

I understand SNR to be the ratio of average signal power to RMS noise power:
$$
SNR=\frac{S_P}{N_P}
$$

or in dB:
$$
SNR_{dB}=10\cdot\log(S_P)-10\cdot\log(N_P)
$$

When calculating SNR based on voltages, power is proportional to voltage squared ($P=V^2/R$) and the Rs drop out in the ratio so we get:
$$
SNR=\frac{{S_V}^2/R}{{N_V}^2/R}=\left(\frac{{S_V}}{{N_V}}\right)^2
$$$$
SNR_{dB}=20\cdot\log(S_V)-20\cdot\log(N_V)
$$

With the familiar 10log for power, 20log for voltage relation doing the squaring for voltage.

I have a radio signal that I receive at an antenna.  For simplicity, I assume I have a perfect receiver and noiseless amplifier so that my received signal power is $P_{Rx}$ and noise is only Johnson thermal AWGN: $N_T=k_B\cdot T\cdot B$, where $k_B$ is Boltzman's constant, T is temperature and B is bandwidth-- none of which really figure in after this.

I believe my SNR at this point to be:
$$
SNR_{RF}=\frac{P_{Rx}}{N_T}
$$

Now I sample, and this is where I start to have questions.

Ignoring quantization noise and such, I believe the analog-digital conversion process results in voltage signals-- that is, the sample values are measurements of the voltage across a resistor ladder or some such meaning the SNR of the sampled signal follows the voltage law:
$$
SNR_{ADC}=\left(\frac{S_{ADC}}{N_{ADC}}\right)^2
$$

Setting aside practical implementation losses, I believe $SNR_{ADC}=SNR_{RF}$.

Where I start to be less sure is when I start operating on the sampled signals.  Let's say I multiply the sampled signal by a delayed version of itself.  Ok, the noise terms get more complicated because I'm taking the product of two independent random variables with non-zero mean but, more fundamentally, is the result a "voltage" or a "power"?  Is there a physical explanation that will help me understand this?

That is: in order to maintain consistency among my SNR estimates, is this a 10log or 20log calculation?

By a pure units analysis, I should have voltage-squared which implies power-- but these are still ADC levels.  It would also seem odd to say that my signal is voltage-cubed if I multiply by two delayed copies...

Handling units in signal processing is tricky (what are the units of $e^{s(t)}$?). I find it useful to think about operations on signals physically, instead of trying to make sense of them mathematically.

For instance, a mixer multiplies two signals. Let's assume the signals are measured in volts. The mixer output is, obviously, also a voltage. This is true even when one of the mixer inputs is its own delayed output.

When calculating an SNR, you square the signal sample (which is a voltage value). In this case, you get a result in watts, because it turns out that the signal power is equal to its square. Physically, you could connect a watt-meter to get a reading of the signal power; you would find that, at every instant, the power turns out to be equal (or proportional) to the square of each sample.

So why, when multiplying two voltages, you get a voltage in one case (mixer) and a power in another (SNR)? I think the reason is that you interpret the result differently in each case. In the case of SNR calculation, the square of a sample physically matches its instantaneous power, so you can assign that meaning to the mathematical operation.


Help with literature pointers for improving estimation of underlying mean for signal with periodic noise I am in the starting phase of a master thesis and are having some trouble finding a good starting point in the literature for my thesis.
The aim of my thesis is to provide as good of an estimate as possible to an underlying constant signal, given measured data with periodic noise.
My degree is in informatics, with a specialization in machine learning, so I don’t have a strong background in signal processing.
The naïve approach is to just take the average of the measured data, it’s clear that the main issue with this approach, is the effect of incomplete wavelengths.
Of course if one could measure long series the periodic noise would average out, the goal of my thesis is rather to see how short it is possible to measure and still accurately estimate the underlying mean value.
A more sophisticated approach is to apply a filter before averaging the signal, but most filters are designed for this purpose so it is not clear which filters would be best suited.
Many filters does not alter the average at all, but what we need is to flatten out the waves such that the incomplete waves at the ends does not affect the signal.
One approach is tapering the edges of the signal, but we haven’t found any literature on that approach.
we will look into how to model the periodic noise as sine waves and then subtracting the periodic noise from the measured data. This approach works very well if you know the frequencies of the periodic noise but seems harder to do if little is known about the periodic noise.
I have looked a bit into using the autocovariance to sample the signal such that each wave top is accompanied by a wave bottom.
I had a bit of an unfounded hunch about Fourier analysis, but it seems our signal might not be of sufficient quality for this. In any case, I just need to find some literature tackling this problem, so I have some background on the solutions others have made before.
I would think this is a problem occurring with many different measurements, and that it would have been tackled by someone previously, but I seem to have a bit of trouble navigating the literature.

You're looking at a cyclostationary random process. That's a concept you'll really need to understand! This will require a relatively solid basis on Fourier analysis, quite likely, but you might probably just take excourses into math basics literature while you read up on analyzing random processes and then specifically cyclostationary processes (really: do it in that order!).
So, Steven M. Kay: Fundamentals of Statistical Signal Processing: Estimation Theory is definitely a book you'd want to read, at least the first half, if this is the topic of your master thesis. That'd be the foundation of your "state of the art" chapter in your thesis, and it's really a good standard book that explains things in a very canonical way. Chances is your university library has it, or you can find it online for cheap (ask the library, and if that doesn't work, ask your professor to buy that book, by the way, if the university library doesn't have it. It's really, as far as I can tell, the standard book on the topic.)

with a specialization in machine learning

That means that you probably have a very solid foundation in stochastics, and my experience is that often, the relatively difficult jump from single (or families of) scalar random variables to stochastic processes has already been done, but often under some other name.

The naïve approach is to just take the average of the measured data, it’s clear that the main issue with this approach, is the effect of incomplete wavelengths.

Damn straight!

Of course if one could measure long series the periodic noise would average out, the goal of my thesis is rather to see how short it is possible to measure and still accurately estimate the underlying mean value

Good news: Well, I must admit this sounds either very unspecific, or very basic! In other words, soon as you refresh your stochastics knowledge and maybe extend it to a theoretical understanding of stochastic processes, you might be making relatively big jumps forward.
For example, assuming your noise $N(t)$ is additive, i.e. your observation $Y(t) = X(t)+N(t)$, and seeing that integration is a linear operation, you'll be pretty quickly be able to calculate the variance of your observation if your calculate the integral over some window $w(t)$, $\int_a^b Y(t)w(t)\,\mathrm dt=\int_a^b X(t)w(t)\,\mathrm dt+\int_a^b N(t)w(t)\,\mathrm dt$. Then, you'd throw a bit of the definition of cyclostationarity, and maybe Fourier analysis at the problem :)

One approach is tapering the edges of the signal, but we haven’t found any literature on that approach.

You must be reading the wrong literature! Thomson's Multitaper Method is indeed one of the classics in the analysis of random, periodic signals.

we will look into how to model the periodic noise as sine waves and then subtracting the periodic noise from the measured data.

You can never subtract noise unless you know it, that's the thing about noise, so you'll need to estimate it first, and that's a spectral estimation problem.


Frequency estimation on chaotic waveform The target signal is generated by the following analytically solvable hybrid chaotic oscillator

\begin{equation}\label{eq:1}
\ddot{u}  - 2\beta \dot{u}  +(\omega^2+\beta^2)(u-s)=0
\end{equation}

where

$\omega$ determines the symbol period as $T = \frac{2\pi}{\omega}$.
$\beta \in (0,T^{-1}\ln 2)$ affects the fluctuating amplitude.
$s \in{\pm 1}$ is the random bipolar symbol sequence driving the system.

Its typical continuous waveform $u(t)$ for $f=100 \, \mathrm{Hz}$ is as follows:

I am interested in estimating the basis frequency $f$ of the signal. However, the amplitude spectrum of $u(t)$ has no obvious feature that can be used to estimate the frequency. Hence, I calculate

$$z(t) = |u(t)| - 1$$

The amplitude spectrum of $z(t)$ is as follows

Obviously, we can see a spike at $f = 100 \, \mathrm{Hz}$. But it is hard to calculate the accurate expression of $z(f)$ since $u(t)$ is a kind of non-stationary signal.

The STFT time-frequency spectrum is shown below

As shown in the STFT results, in some time intervals, no obvious frequency component can be found. However, in the part of approximately periodic waveform, the frequency component at $f = 100 \, \mathrm{Hz}$ is strong.

Is there any theory that can support my method of frequency estimation from the FFT of the whole $z(t)$?

Here is the MATLAB code for generating the chaotic waveform from an analytical expression (not provided above) and my current experiments on frequency estimation.

clear;close all;

finit = 100;%set frequency
% finit =2000;
% finit =2;
om=2*pi*finit;
be=om/(2*pi)*log(2)-0.01*rand(1,1);
T = 2*pi/om; %Period
N=100;% number of symbol
tall = N*T;%total time
h=T*0.01;%time step
t=0:h:tall-h;
ut=zeros(1,length(t));
u_in = 2*rand(1)-1;
% u_in =0.99;
for i=1:N
    if u_in>=0
        s_in = 1;
    else
        s_in = -1;
    end
    t_f = (i-1)*T:h:i*T-h;  %Time of each frame
    n_f = (i-1)*length(t_f)+1:i*length(t_f);
    ut(n_f)=s_in+(u_in-s_in)*exp(be*(t_f+(1-i)*T)).*...
           (cos(om*t_f)-be/om*sin(om*t_f));
    u_in=s_in+(u_in-s_in)*exp(be*T)*cos(om*i*T);
end

figure
plot(t, ut)
xlabel('Time t/s');
title('u(t)')

y = abs(ut) - 1; % Observed signal

%% FFT spectrum
L = length(ut);
Ts_v = t./(0:L-1);
Fs = ceil(1/Ts_v(2));
deltaf= 0.1;
NFFT = floor(Fs/deltaf);
fs = Fs;
f = Fs/2*linspace(0,1,floor(NFFT/2+1));
Y = fft(y,NFFT)/L;
Yabs = 20*log10(abs(Y(1:floor(NFFT/2+1))));
figure;
plot(f,Yabs);
xlim([0,max(f)/4]);
ylim([-100,0]);
xlabel('Frequency f/Hz');
title('FFT Spectrum');
%%  Wavelet Transfrom Time-Frequency spectrum
wavename='cmor3-3';
totalscal=256*2;
Fc=centfrq(wavename);
c=2*Fc*totalscal;
scals=c./(1:totalscal);
f=scal2frq(scals,wavename,1/fs);
coefs=cwt(y ,scals,wavename);

figure;
subplot(211)
plot(t,y,'linewidth',1)
xlim([0 1])
ylim([-1.5 1.5])
subplot(212)
imagesc(t,f,abs(coefs));
set(gca,'YDir','normal')
ylim([0 200]);
% colorbar;
xlabel('Time t/s');
ylabel('Frequency f/Hz');
title('Wavelet Time-Frequency Spectrum');
%% STFT Time-Frequency spectrum
window_len = 10*2*pi/om/h;
noverlap = 50;
nfft = NFFT;
[~,F,T,P]=spectrogram(y,hamming(window_len),noverlap,nfft,fs);

figure;
subplot(211)
plot(t,y,'linewidth',1)
xlabel('Time t/s');
title('|u(t)|-1')
xlim([0 1])
ylim([-1.5 1.5])
subplot(212)
surf(T,F,abs(P),'edgecolor','none')
axis tight;
view(0,90)
ylim([0 200]);
xlabel('Time t/s');
ylabel('Frequency f/Hz');
title('STFT Time-Frequency Spectrum');
```

Here is a stab.

Start with:

$$ \ddot{u}  - 2\beta \dot{u}  +(\omega^2+\beta^2)(u-s)=0 $$

Take the derivative:

$$ \dddot{u}  - 2\beta \ddot{u}  +(\omega^2+\beta^2)(\dot{u}-\dot{s})=0 $$

Now, $ \dot{s} = 0 $ Wherever it is defined.  So:

$$ \dddot{u}  - 2\beta \ddot{u}  +(\omega^2+\beta^2)(\dot{u})=0 $$

Let $ v = \dot{u} $, and you get

$$ \ddot{v}  - 2\beta \dot{v}  +(\omega^2+\beta^2)v=0 $$

Which is an ordinary second degree linear differential equation.  v will have the same frequency as u.

The most basic (good) differential approximation is:

$$ d[n] = \frac{u[n+1]-u[n-1]}{2} $$

So run that over your data and see what you get and please add a chart to your question with the results.

I think you may best solve this one by analyzing the zero crossing of those results.  The transitions are going to cause weirdness as your system "falls" towards the new equilibrium.

Followup:

A more traditional approach perhaps.

$$ \ddot{u}  - 2\beta \dot{u}  +(\omega^2+\beta^2)(u-s)=0 $$

Your homogenuous equation is:

$$ \ddot{u}  - 2\beta \dot{u}  +(\omega^2+\beta^2) u=0 $$

with solutions:

$$ u_h = \left[ C_1 \cos( \omega t ) + C_2 \sin( \omega t ) \right] e ^ { \beta t } $$

And your particular equation is:

$$ \ddot{u}  - 2\beta \dot{u}  +(\omega^2+\beta^2) s=0 $$

Assume:

$$ u = a t^2 + b t + c $$

$$ \dot{u} =  2 a t + b $$

$$ \ddot{u} =  2 a  $$

$$ (2a)  - 2\beta ( 2 a t + b )  + (\omega^2+\beta^2) s=0 $$

$$ (2 - 4 \beta  t)  a  - 2\beta  b   + (\omega^2+\beta^2) s=0 $$

$$ a = 0 $$

$$ b = \frac{\omega^2+\beta^2}{2\beta} \cdot s $$

So :

$$ u_p = \left( \frac{\omega^2+\beta^2}{2\beta} \cdot s \right) t $$

For any range of consistent $s$ you now have your general and particular solutions.

Therefore you should only need five points to solve for your five unknowns:  $ \omega, \beta, C_1, C_2,$ and $ c $.

With noise, that goes out the window.  If you take the second derivative, your particular solution will be squashed completely.  Then you should only have to worry about the $s$ transition points.

Line crossings will give you $\omega$, and the envelope will give you $\beta$.


Good test for periodicity between signals I have two timeseries signals. They look like this:

Each signal started out from the same array, but each received different preprocessing treatments. Ultimately, each signal represents the breathing rate of a subject. Each peak constitutes the peak of an inhale, so the number of peaks/time yields the respiratory rate. I want to determine which preprocessing treatments yield a more periodic signal so that I can conduct the most optimal peak detection.
As you can see, both signals are noisy, but it looks like the orange timeseries has more consistent periodicity, making it more suitable for my task. What is an optimal test to determine how periodic each signal is so that I can determine the signal without having to look at them?
Ideally, a solution implemented in python or c++ would be great, but not necessary.
Edit
Here is a way to compute average frequency per signal.

Am I correct to assume the higher value is more periodic, or is that too simplistic?
Edit 2
More "periodic" probably isn't the right term. This is what I mean:

The first plot has less noise, so it is cleaner/easier to determine the periods between peaks. Again, I am ultimately trying to solve a peak detection problem.

Something you could do is calculate the FFT of both signals and define some criterion for periodicity. For example, you pick the highest component of each FFT and compare the two (one from each signal) and thus determine which is more "periodic". You could consider an average of all the FFT components and compare. This would be a rather robust solution, since the FFT calculates the frequency decomposition of a signal.


How do I normalize a predesigned IIR filter? Obviously  my question here is very similar to  this previously asked question: How to normalize predesigned IIR filter coefficients to have unity gain?

Unfortunately, in the responses available to that question, it isn't super clear how the answers are obtained.
I have a pre-designed IIR filter that is applied in 2 steps:

Step 1:

    Yn =  Xn + a1*Yn-1 - b1*Yn-2

Step 2:

    Zn =  Yn + a2*Zn-1 - b2*Zn-2

with Xn being my initial input and Zn my final output.
Given these 2 pairs of coefficients (a1,b1) and (a2,b2), how do I go about finding the gain of this IIR filter to ensure it has unity gain.

Thanks for the help

As mentioned in the comments, you have to choose a frequency $\omega_0$ at which you want to normalize the gain. This could be DC (i.e., $\omega_0=0$) or any other frequency, depending on the filter's characteristic.

You have two filters in series, each with transfer function

$$H_i(z)=\frac{1}{1-a_iz^{-1}+b_iz^{-2}},\qquad i=1,2\tag{1}$$

The total transfer function is just the product of the two individual transfer functions:

$$H(z)=H_1(z)\cdot H_2(z)\tag{2}$$

The gain at frequency $\omega_0$ equals the magnitude of the total transfer function $(2)$ evaluated at $z=e^{j\omega_0}$. So if you scale one of the two transfer functions $H_i(z)$ by the inverse of that gain, your total system will have unity gain at $\omega_0$.

Scaling in your time domain equations simply means scaling the input variable.

Note that $\omega_0$ is a normalized frequency in radians:

$$\omega_0=2\pi\frac{f_0}{f_s}\tag{3}$$

where $f_s$ is the sampling frequency.


Formulating a function on Matlab for the Shannon interpolation formula I am trying to formulate an algorithm for applying to Shannon interpolation formula to the discrete signal

$$x[n]=\frac{c^2}{4}\int\limits_{0}^{n T}y\left(\tfrac{c}{2}s\right)\,ds,$$

where $c$ is constant. Now, if $\frac{1}{T}=f_{s}>2B$, where $B$ is the band-limit of $X(f)$, then we can reconstruct the continuous time signal as:

$$x_{\text{const}}(t)=\sum\limits_{n=-\infty}^{\infty}\frac{c^2}{4}{\int\limits_{0}^{n T}y\left(\tfrac{c}{2}s\right)\,ds}\cdot\text{sinc}\left(\frac{t-n T}{T}\right)$$

function x=interpolate(y,N,T,c)
     x=0;
     for n=1:N
         x=x+((c^2)/2)*integral(y,0,n.*T).*sinc((t-n.*T)./T);
     end

However, I am having trouble with incorporating the t variable here. Setting it as global and running interpolate(y,20,20,1) also produces nothing to note.

Note that the continuous time signal in question is given by

$$x_{\text{true}}(t)=\frac{c^2}{4}\int_{0}^{t}y\left(\tfrac{c}{2}s\right)\,ds.$$

Given your samples $x[n]$, here's a code to perform the interpolation by the sinc-functions:

F = 1000; % the sampling frequency used for drawing the continuous function
Fs = 50; % the sampling frequency
T = 1/Fs; % the sampling time

t = 0:1/F:1; % the time-samples for "continuous" plotting
ts = 0:T:1;  % the time at the sample points

% the original function. In your case, this would be integral expression
x = @(t) sin(2*pi*10*t) + cos(2*pi*11*t);
xn = x(ts);  % the discrete sequence

% perform interpolation from xn to x:
interpolated = 0;
for n=1:length(ts)
    % n-1 due to the 1-based indexing of matlab
    interpolated = interpolated + xn(n) * sinc((t-(n-1)*T)/T);
end

hold off;
plot(t, x(t), 'r', 'linewidth', 3); hold on;  % plot continuous
stem(ts, xn);  hold on; % plot sampled
plot(t, interpolated, 'k-o'); % plot interpolated on top of red curve

As a bonus, here's what happens when the signal is not enough bandlimited: I changed the frequency of the first sine from $sin(2\pi 10 t)$ to $sin(2\pi 30 t)$:

Note: The Shannon-Interpolation requires bandlimited signals, so in particular it does not work exactly with time-limeted signals (since these are never band-limited). So, you will see some deviations especially at the edge of the signal.


Doubt: Signal detection in noise Part 1 I have a confusion regarding the matched filter detection technique. I have a binary information source signal $s(t)$ that is corrupted by additive White Gaussian noise $w(t)$ at a particular SNR. The received signal is:
$$x(t) = s(t) + w(t)$$
Then I have created a matched filter $h(t)$ as the time reversal of the source binary signal $s(t)$. How does a matched filter work?

You use the matched filtering technique when you search measurement data for a signal of a given form. You know the signal waveform in advance and either you are not sure if the signal is present in your data, or you see the signal and your task is to compute the signal magnitude and location. In these cases you use the matched filtering technique.
If your measurement data contains an unknown signal buried in noise, you can try and recover the waveform of this signal, with additional assumptions made about the noise nature. Often, this noise is additive white Gaussian noise, AWGN. The signal mixed with AWGN is amenable to processing with filters; filtering gives an approximate waveform of the "pure" signal. The precision of filtering -- the proximity of a recovered signal to the "pure" signal -- varies depending on the signal and noise parameters and the filter type and implementation.
Summing up: you use matched filtering when you know the signal waveform in advance or test you data to establish the presence of hypothetical signal waveforms. The waveform you test against is called template. If you want to extract the signal waveform from noisy data, you do not use matched filtering; you apply filters.


Why is particular solution zero for an impulse excitation signal? We were being taught the impulse response for a series RC Circuit- consisting of simply one resistance, one capacitor, and an impulse excitation all in series.
I get that the homogeneous part of the solution will be of the form:
$$x(t) = Ke^{\frac{-t}{RC}}$$
And we can find the value of $K$ using the initial conditions.
However, when it came to the particular solution part, we were straightaway told that since the excitation is an impulse, the particular solution will be equal to $0$. I'm not able to get why is this the case.
Also, I am not getting exactly how to calculate the particular solution for a given network, say RC network.

The correct form of the statement should be

the particular solution is zero for $t>0$

This is simply the case because the input $\delta(t)$ is zero for $t>0$.
So for $t>0$, the impulse response equals the homogeneous solution. For $t<0$, the impulse response is zero because the system is causal. What remains is the moment $t=0$. The only thing that can happen there is an impulse (or its derivatives but let's not get into that here). So the complete form of the impulse response is
$$h(t)=A_0\delta(t)+y_h(t)\tag{1}$$
where $y_h(t)$ is the homogeneous solution of the differential equation. The constant $A_0$ can be determined by looking at the highest derivative of the input on both sides of the differential equation. For the given system we have
$$h'(t)+\frac{h(t)}{RC}=\frac{\delta(t)}{RC}\tag{2}$$
With $h(t)$ given by $(1)$ we see that on the left-hand side of $(2)$ we have a term $A_0\delta'(t)$ with nothing to match on the right-hand side. Consequently, for this specific differential equation we obtain $A_0=0$.
Note that it's easy to find a system for which $A_0\neq 0$. E.g., for a first-order RC highpass filter we have the following differential equation for the impulse response:
$$h'(t)+\frac{h(t)}{RC}=\delta'(t)\tag{3}$$
Using $(1)$ and matching the terms with $\delta'(t)$ on both sides of $(3)$ we obtain $A_0=1$, resulting in
$$h(t)=\delta(t)-\frac{1}{RC}e^{-t/RC},\qquad t\ge 0\tag{4}$$


Logaritmic amplitude squarewave This morning i was experimenting with GNU Octave.
I created several signals(Sine, Saw tooth, square wave and noise). After this i used the periodogram. in the periodogram the square wave weren't showed.
I tried a test file from where the square wave is broken.

A = 1;
f = 2;
T = 2;
fs = 512;

if isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;
    pkg load signal
    clear isOctave;
end

t = 0:1/fs:T-1/fs ;
sine = A * sin(2*pi*f*t);
sq =  A * square(2*pi*f*t);
F = 0:(fs/(T*fs)):fs/2;
SINE = fft(sine);SINE = abs(SINE(1,1:(T*fs)/2+1));
SQ = fft(sq); SQ = abs(SQ(1,1:(T*fs)/2+1));

figure;
subplot(2,1,1); plot(t,sine, t,sq);
legend('Sine', 'Square','Location','northoutside','Orientation','horizontal')
subplot(2,1,2); plot( F, 20*log10(1/(length(SINE))*SINE), F,20*log10(1/(length(SQ))* SQ));

figure
periodogram(sine,[],length(sine),fs,'power'); hold all
periodogram(sq,[],length(sq),fs,'power')
legend('sine','square')

What i found out is that till the fft everything works as suspected.
But when i take the log of the square wave. the result disappears.

My question is why does this happening?

what cause this problem.
Happens this also in Matlab or python.
What solution fix this problem?

The library function log of zero results in a NaN (not a real number, -inf).  So the log of a square with a base of zero can't be plotted.


Why do we get ripples in the frequency spectrum of a chirp, and how do I remove it? The ripples in the amplitude spectrum of a linear spectrum is undesirable. I desire to obtain a flat spectrum response. Hence, how can I remove these ripples? Also, which factors (like time bandwidth product, chirp rate) are affecting the ripples, and how?

I designed a chirp with a frequency from 10 to 90 Hz in MATLAB. However, I obtained Fresnel ripples (which normally occurs) in the passband, along with overshoots at the band edges. I need to study the application of chirps on a material and use it for defect analysis. For this, I need to optimize the whole energy applied to the chirp. Hence, I want to get rid of these ripples. So, what parameters like chirp rate and time of sweep will affect the system?

Another thing the OP can do is use the "complex" chirp. What you really have to do is (at two different times and both times sweeping all frequencies from -Nyquist through DC to +Nyquist) pass the real part, $x_r(t)=\cos(\pi \beta t^2)$ and then pass the imaginary part $x_i(t)=\sin(\pi \beta t^2)$ of the complex chirp through the material under test and synchronously measure the two real results ($y_r(t)$ and $y_i(t)$), then assemble them as the real and imaginary parts of the response.

Let
$$ \begin{align}
x(t) & =  e^{j \pi \beta t^2} \\
& = \cos(\pi \beta t^2) + j \sin(\pi \beta t^2) \\
& = x_r(t) \qquad + j \ x_i(t) \\
\end{align} $$

Then
$$ X(f) \ \triangleq \ \mathscr{F}\left\{ x(t) \right\} \ = \ \sqrt{\frac{j}{\beta}} e^{-j\frac{\pi}{\beta} f^2} $$

You will see that $ |X(f)| = \sqrt{\frac{1}{\beta}} $ for all $f$.

The respective outputs that are synchronously measured with the two driving signals are $y_r(t)$ and $y_i(t)$:

$$ \begin{align}
y_r(t) & \triangleq h(t) \circledast x_r(t) \\
 & = h(t) \circledast \cos(\pi \beta t^2) \\
\\
y_i(t) & \triangleq h(t) \circledast x_i(t) \\
 & = h(t) \circledast \sin(\pi \beta t^2) \\
\end{align} $$

Then this aggregate complex output, $y(t)$, can be assembled (in the mind of the computer) from the responses of the cosine sweep (the real part) and the sine sweep (the imaginary part):

$$ \begin{align}
y(t) & =  y_r(t) + j \ y_i(t) \\
& = h(t) \circledast x_r(t) + j \ h(t) \circledast x_i(t) \\
& = h(t) \circledast \left(x_r(t) + j \ x_i(t)\right) \\
& = h(t) \circledast x(t) \\
& = \int_{-\infty}^{+\infty} h(u) \, x(t-u) \ du \\
& = \int_{-\infty}^{+\infty} h(u) \, e^{j \pi \beta (t-u)^2} \ du \\
& = \int_{-\infty}^{+\infty} h(u) \, e^{j \pi \beta (t^2-2tu+u^2)} \ du \\
& = e^{j \pi \beta t^2} \int_{-\infty}^{+\infty} \left(h(u) \, e^{j \pi \beta u^2}\right) \, e^{-j 2 \pi \beta t u} \ du \\
& = e^{j \pi \beta t^2} \int_{-\infty}^{+\infty} \tilde{h}(u) \, e^{-j 2 \pi (\beta t) u} \ du \\
& = e^{j \pi \beta t^2} \tilde{H}(\beta t) \\
\end{align} $$

where
$$ \tilde{H}(f) \triangleq \int_{-\infty}^{+\infty} \tilde{h}(t) e^{-j 2 \pi f t} \ dt = \mathscr{F}\left\{\tilde{h}(t) \right\} $$

and
$$ \tilde{h}(t) \triangleq h(t) e^{j \pi \beta t^2} $$

So this is what you get: first consider making this adjustment (dechirping) in the measured output:
$$ y(t) e^{-j \pi \beta t^2} = \tilde{H}(\beta t) $$

Now if the sweep rate $\beta$ is low enough, then $e^{j \pi \beta t^2} \approx 1$ for as long as $h(t)$ is not close to 0. Then $\tilde{h}(t) \approx h(t)$ and $\tilde{H}(f) \approx H(f)$ and $\tilde{H}(\beta t) \approx H(\beta t)$ and what you have measured in the output is another chirp, but with relative magnitude and phase changed by $\tilde{H}(\beta t) \approx H(\beta t)$.

Even if the sweep rate is not so low, you can compute either the true frequency response $H(f)$ by first computing

$$ \begin{align}
\tilde{h}(t) & =  \mathscr{F}^{-1}\left\{ \tilde{H}(f) \right\} \\
& = \int_{-\infty}^{+\infty} \tilde{H}(f) e^{j 2 \pi f t} df \\
& = h(t) e^{j \pi \beta t^2} \\
\end{align} $$

Or from that result you can get the true impulse response $h(t)$ and then $H(f)$ if that's what you want. It should be the same $H(f)$ you would get if you divided $Y(f) = \mathscr{F}\left\{ y(t) \right\} $ by $X(f) = \mathscr{F}\left\{ x(t) \right\} = \sqrt{\frac{j}{\beta}} e^{-j\frac{\pi}{\beta} f^2}$. And there is no division by zero, because there are no ripples in $X(f)$ because the magnitude is a constant $\sqrt{\frac{1}{\beta}}$.

That is how to remove ripples (and any danger of division by zero) in the spectrum of a linearly-swept chirp.


The role of GPS in INS/GPS navigation systems Ideally, a gyroscope and an accelerometer would be enough for a complete navigation solution (attitude + position), using dead reckoning. This comprise the Inertial Navigation System, INS. In non-ideal world, we would couple the INS with GPS. The GPS helps 'fix' the navigation solution at certain times, to help mitigate the drift in the INS output.

My question.
Suppose I have no access to a GPS, but have access to a Spirit level instrument with associated camera and image processing algorithms. This system serves to provide an absolute values of pitch and roll angles, but with slower update rate and lower resolution. Still, it is assumed to be drift-free. Can such a system be used in place of normal GPS to achieve dead reckoning solution?

P.S. The reason for not using GPS is that I am trying to find a self-supporting navigation solution...

By itself, no, it cannot.

Your method of fixing the INS system is pretty clever, but it only accounts for drift in the gyroscope. Accelerometers are also (very) prone to drift, generally showing up as an increasing velocity bias. The typical method to account for accelerometer drift is to use zero-velocity updates. By detecting points when you are not moving, you can periodically remove any non-zero velocities from your model.


How to get Harris corners that are present in watershed segmented Bitwise image(road part)? This is Harris corner image of test image:

This is Bitwise watershed segmented image:

I want to get that Harris corners which will be in segmented part of image(road part).

You have a mask of the road were pixels of the road have value 1 and non road pixels as value 0.
Also you have a list of Harris corners were each corner as a coordinate Xi,Yi

For each Xi,Yi i corners
If mask[Xi,Yi] == 1
keep corner


Is it appropriate to use PSD on the analysis of a plucked guitar string? I'm using python 2.7 to compare tonal differences in guitar stings.

I've found that in matplotlib, there are some great built-in features such as
matplotlib.mlab.magnitude_spectrum
and
matplotlib.mlab.psd

The magnitude_spectrum function seems to be the no brainer, as it is the FFT of the signal, showing each specific harmonic. However, when comparing various guitar strings, it's tricky to pinpoint each peak to compare. The PSD function gives very nice smooth graphs, which is perfect for obvious visual comparisons.

My understanding is that PSD's are useful for random signals but not necessarily known signals. Is this true? And is there any reason I shouldn't be using the PSD for guitar string comparisons?

Here's what the PSD and magnitude_spectrum graphs look like:

Note that these are plucks of the same string and not different strings. I'm trying to make sure that my plucks are consistent prior to comparing various strings.

I'd really recommend going for a short-time fourier transform such as that produced by: matplotlib.mlab.specgram. It gives you the change in the spectrum over time and can be calculated using PSD or the magnitude spectrum.

The reason I'd recommend this is that your guitar signal is not spectrally stationary; your graphs are only a snapshot at one point in time. What gives your guitar string their tonal differences is the way all the harmonics shown above change over the length of the note.

There's a very good summary of PSD vs FFT magnitude here and here.


How to compute bit level soft decisions from an M-ary symbol? When working with a decoder such as LDPC that can take in soft-decision inputs (or channel posterior probabilities as I've seen them called as well), how are these soft-decision values calculated for general M-ary modulations like M-FSK, PSK, or QAM?  I have only been able to find the BPSK case so far in Todd Moon's Error Correction Coding Book. I'm assuming AWGN channels and binary FEC codewords over GF(2).  Sources to information are fine too if the answers to this would be too long.

In each case, the soft-decision output is just the matched filter's output.

For example, say you transmit the bit sequence $1,1$ with codeword $1,0,1,1$. You transmit these bits using, say, QPSK symbols $0.5+0.5j,0.5-0.5j$. After matched filtering in the receiver, you will get two noisy complex numbers, for example $0.55+0.45j,1.2-0.49j$. A hard-decision decoder will turn those numbers into QPSK symbols and then bits before FEC. A soft-decision decoder will use those numbers without any further processing (except maybe some quantization).


MATLAB and sampling rate I have a question about sammpling rate in MATLAB.
I read a .wav file using audioread(). It give me sampled data called my_data, and $88.2 \textrm{ kHz}$ as a samping rate, called Fs.

When I played this data with audioplayer(my_data, Fs), it works well and totally right sounds with original file.

What I tried to do is getting odd-indexed data like this :

my_data = [1 2 3 4 5 6 7 8 9 10]
my_odd_data = [1 3 5 7 9]

I thought that this procedure is like sampling my_data with half of Fs, but it is not true. I want to know why this procedure is not sampling with Fs/2.

As people have said in the comments, to have it sound right you will need to low pass filter before subsampling.

The matlab function decimate does this automatically and gives some control over the type of low pass filter used.


Pitch estimation using the autocorrelation method I have the following code:

// input: speech segment
[y,Fs,bits]= wavread('/media/A03036C33036A068/scilab/30msec_voiced.wav');
max_value=max(abs(y));
y=y/max_value;
t=(1/Fs:1/Fs:(length(y)/Fs))*1000;
subplot(2,1,1);
plot(t,y);
xtitle('A 30 millisecond segment of speech','time in milliseconds');
sum1=0;autocorrelation=0;
for l=0:(length(y)-1)
    sum1=0;
    for u=1:(length(y)-l)
        s=y(u)*y(u+l);
        sum1=sum1+s;
    end
    autocor(l+1)=sum1;
end
kk=(1/Fs:1/Fs:(length(autocor)/Fs))*1000;
subplot(2,1,2);
plot(kk,autocor);
xtitle('Autocorrelation of the 30 millisecond segment of speech','time in milliseconds');
auto=autocor(21:160);
max1=0;
for uu=1:140
    if(auto(uu)>max1)
        max1=auto(uu);
        sample_no=uu;
    end
end
pitch_period_To=(20+sample_no)*(1/Fs)
pitch_freq_Fo=1/pitch_period_To

In this, I am not getting the meaning of pitch_period_To=(20+sample_no)*(1/Fs). Why did they add 20?

An offset of 20 samples is added as 8 lines above (auto=autocor(21:160)) this offset is cut away from the autocorrelation sequence. In the line you mention, this offset has therefore to be considered.

This is done in order to cut away the first maximum (at time lags around zero) and avoid pitch errors caused by picking a delay value that does not make sense for the task of fundamental frequency estimation.

(By the way, your code does not perform pitch estimation but fundamental frequency estimation. There's a difference between these two terms as pitch relates to a perceptual property of a tone while fundamental frequency is a technical property - Wikipedia on "Pitch".)


Is it a good idea to whiten a colored noise in order to apply formulas for AWGN? I have posted a question before regarding the detection of a signal under colored Gaussian noise:
Energy detection in presence of colored Gaussian noise

Theoretical formulas for setting the detection threshold are widely used for AWGN.

Is it a good idea to just perform a pre-whitening to the colored noise and then apply the well-know formulas for energy detection under white Gaussian noise?

Do you have any happy experience using this approach in practice?

The short answer is yes, whitening noise is usually recomended (if you know the spectrum of the noise) but whitening also whitens your signal.  Detecting random signals in colored noise is a well known problem is sonar signal processing.  There are a number of derivation that I like particularly:

Nielsen, Richard O. Sonar signal processing. Artech House, Inc., 1991.

and Burdic

@book{RN17,
   author = {Burdic, William S.},
   title = {Underwater acoustic system analysis},
   publisher = {Prentice Hall},
   ISBN = {0139476075},
   year = {1991},
   type = {Book}
}

Also the updated

@book{van2004detection,
  title={Detection, estimation, and modulation theory, part I: detection, estimation, and linear modulation theory},
  author={Van Trees, Harry L},
  year={2004},
  publisher={John Wiley \& Sons}
}

in section 4.3

and if you can find it, in the old volume 3

If you look, a copy of Eckart's  classic paper can be found at Scripts.


High audio frequencies in television broadcast I would like to know if in television broadcast high end audio frequencies are transmitted e.g. frequencies higher than 14 kHz. What is the cut off audio frequency of transmission in television broadcast?

I am preparing a report (for research) characterizing audio frequency response of various devices e.g. iPhones, Samsung Galaxies, various televisions. While doing this, I got curious about the spectrum output of television. Given a piece of Audio for broadcast transmission, how it's spectrum get altered? At what frequency there is cutoff and specially can the cutoff be below ~20K, the human audible range?

Broadly my interest is in knowing what happens to an audio given to a broadcasting channel and by the time it arrives as an input to a television.

Thanks.

In digital television, I think they use mostly a 48 kHz sampling frequency, so it is up to whether the audio encoder allocates any bits to the ultrasonic frequencies:

Figure 1. Before encoding: Left) almost full scale white noise, Center) almost full scale pink noise, Right) full scale 20-24 kHz band pass white noise.

Figure 2. After encoding in MPEG 1.0 layer II, 224 kbit/s, 48000 Hz stereo (known to be used in DVB-T) and decoding, using mp2enc for encoding and mpg123 for decoding.

With white noise, in the 16-20 kHz range only some bursts of frequency survive, but above that nothing. Starting with "pink" noise, everything above 15 kHz gets stripped away, because more bits are allocated to the higher-intensity lower frequencies. Starting with content only >20 kHz, some ultrasonic frequencies survive, as no bits need to be allocated to lower frequencies.

Also other encoding methods than MPEG 1.0 layer II are in use in broadcast television.


Producing images with subpixel shift in region of interest Let's say for a zebra image, and at certain pixel location, I would like to makes it "vibrate/warp" to and fro in horizontal direction at a certain frequencies in a sequence of images of time t, with subpixel shift. The circled area is the most significant area that I intended to simulate the vibration. Those images will be mainly used for synthetic motion estimation.

In MATLAB, the the zebra image can be produced using the following code:

m=75;
n=75;
stripe_px=1;
single_pattern=[ones(m,stripe_px),zeros(n,stripe_px)];
zebra=repmat(single_pattern,[1 ceil(n/stripe_px*2)]);
zebra=zebra(:,1:n);

The magnitude of the subpixel shift indicates by this map below:

Also can be produce by MATLAB code:

X=repmat(1:m,[n 1]);
Y=repmat((1:n)',[1 m]);
Shape=0.04*( ((((X/m)-0.75).^2).*X) - ((((Y/n)-0.75).^2).*Y) );

Based on the information that I currently have, what's my next step? I was looking at imwarp (link), but I'm not sure if it can shift local pixels, where usually it was meant for whole image transform. And its usage is very vague for me. A simple pseudocode to depict what I'm trying to achieve would be:

omega=3;
T=30;
Fs=20;
dt=1/Fs;
for t=0:dt:T
    vibrate_im=(zebra + Shape) + sin(omega*t);
end

Can a wavelet transform be used for this subtle pixel shifting? If so how do I manipulate its coefficients?

Appreciate any input and help! Thanks in advance.

I tried to distill your question down to a 1-d problem. Does this "moving sine" do ~what you want to do?

t = 0:50;
f = (1/4)/5;
x = sin(2*pi*f*t);
nsteps = 50;
figure(1),
clf
ylim([-1 1])
xlim([0 nsteps])
stepvect = linspace(0,3,nsteps);
stepvect = [stepvect stepvect((end-1):-1:2)];
for step = 1:length(stepvect)
    y = interp1(t, x, t+stepvect(step)*sin(pi*t/50), 'linear');
    plot(t, y)
    grid on
    drawnow;
end


Difference between the use of 'bandwidth' in DSP and computer networking As far as I know, bandwidth is measured in herz or 'radian/second' in digital signal processing or in analog communications. But in computer networking it is measured in 'bits/second'. I don't get  how these two are related. Can anybody please help?

They both define the amount of data you can get through a channel -- the channel bandwidth, together with the SNR, defines the channel capacity $C$, which has the unit of $\frac{\mathrm{bit}}{\mathrm{s}}$ (which, by the way, is simply $1$ binary decision per second), giving you the maximum amount of information you can get through that:

$ C_\mathit{Shannon} = f_\mathrm{bandwidth} \log_2 ( 1+\mathrm{SNR} )$

You see, bandwidth in the signal sense just says "we can see so and so many changes per second", and you can directly translate that to bits you get through that bandwidth, if you know how well your signal gets through.


Discrete Convolution of two piecewise sequences having this specific form Assume I have the following two sequences :
$$
x[n]=\begin{cases}
\alpha&\text{if $a\leq n\leq b$}\\ \\ \tag{1}
0&\text{if otherwise}
\end{cases}
\qquad
\text{and}
\qquad
h[n]=\begin{cases}
\beta&\text{if $c\leq n\leq d$}\\ \\
0&\text{if otherwise}
\end{cases}
$$
I was wondering in this case if the convolution :
$$
(x*h)[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k] \tag{2}
$$
has a specific form applicable for these two piecewise sequences.
I have seen on the internet something like :
$$
(x * h)[n]=\sum_{k=\max (a, n-d)}^{\min (b, n-c)} x[k] h[n-k] . \tag{3}
$$
However, I can not confirm if this is true. I hope someone can provide a proper proof and thank you very much.

Given your Eq-3
$$
(x * h)[n]=\sum_{k=\max (a, n-d)}^{\min (b, n-c)} x[k] h[n-k] . \tag{1}
$$
where $x[n]$ and $h[n]$ are as in Eq-1, then you will have
$$
\sum_{k=\max (a, n-d)}^{\min (b, n-c)} \alpha \beta = \alpha \beta (\min (b, n-c)- \max (a, n-d) + 1). \tag{2}
$$
Where you should evaluate the closed form expression Eq.2 for the non-zero range of output given by the range of $n$ :
$a+c \leq n \leq b+d $


Confusion about decibels I know that $\text{dB}$ actually represents dimensionless value, meaning, it is the ratio of two values that have the same dimension, for example:

$$p_1=10\log\frac{P_1}{P_0} $$

This can be expressed in decibels if $P_0=1\text{ W}$, however, if $P_0=1\text{ mW}$ then unit we use is $\text{dBm}$.

I am wondering, are we allowed to subtract and add values if one of them is in $\text{dB}$ and another is $\text{dBm}$ ?

The reason i am asking this is because i've found an example in one workbook where two values are summed up where the first one was in $\text{dB}$ and another $\text{dBm}$, however, i don't understand how's that possible.

For example, if we say $p_1=10\text{ dB}$ and $p_2=20\text{ dBm}$ if we sum these two up, we end up with $30$ but how is this possible? It's $30$ what? $\text{dB}$ or $\text{dBm}$? I hope someone could clarify me this. Any help appreciated!

Short answer: $p_\textrm{dB} =\alpha\textrm{ dB}$ means multiplication factor (no unit) while $P_\textrm{dBm} =A\textrm{ dBm}$ is truly a power, nothing more!

Then $p_\textrm{dB} + P_\textrm{dBm} = 10\log_{10}(10^{\alpha/10} \times 10^{A/10})$ is a power, hence has unit $\textrm{dBm}$ in logarithm scale.

$p_\textrm{dB} + q_\textrm{dB} = 10\log_{10}(10^{\alpha/10} \times 10^{\beta/10})$ is factor, hence no unit (or $\textrm{dB}$ in logarithm scale).

$P_\textrm{dBm} + Q_\textrm{dBm} = 10\log_{10}(10^{A/10} \times 10^{B/10})$ is a power square ($\textrm{mW}^2$), unit undefined in logarithm scale.


DFT Signal DFT Length N , FFT If we sample a signal let say sine(2 * pi * f) with f=1Hz  and a sampling Frequency of Fs = 8Hz, is it right that the length of the data should be N = Fs/f or multiple of Fs/f like N= d*(Fs/f) with d=1,2,3,.., is that right ?
I read often, that it is better to increase N,
my question : we can only increase N with Nnew= d * Nold, right ?
and not like Nnew= 1,2 * Nold?
for my example , with f=1HZ, Fs= 8Hz, T=1s, Ts=0,125s, N=8
we can increase to  f=1HZ, Fs= 8Hz, T=2s, Ts=0,125s, Nnew=16,
but it would not be right to increase ist to f=1HZ, Fs= 8Hz, T=1,5s, Ts=0,125s, Nnew=12, wich is not right  because  x[n]=x[n+N] doesnt fit, right ?
......
if it is right above, then where do I know how to get the right length N if i down know my signal ?
Do I have to Look the x[n] Series of the Signal one by one to see where N (full Cycle) of the signal ist ?
Thanks

For a time domain signal the length or the number of samples in the sampled signal can be computed as shown in Equation $(1)$
$$
T = NT_s \equiv T = \frac NF_s \implies N = TF_s\tag{1}
$$
Where $T$ is the signal total duration in [seconds] and $F_s$ is the sampling frequency in [Hertz], and $N$ the total number of samples.
The length of the DFT depends on the number of point taken in computing the DFT. For a $N$-point DFT on a signal of $N$ total samples you have $N$ as the length of your DFT.


What does it mean when you zero pad data and why does one do it? I am struggling to understand what effect zero padding data would have on EMG data. How does it improve the signals? Any help in simple terms that's not too mathematical would be greatly appreciated!

There are a few reasons to zero pad, improvement is mostly in terms of advantages in manipulating the signals, not explicit SNR enhancement.

Delay, Advance (by n), and/or Add a signal of length M with respect to a longer signal of length N.  These are operations used in languages with linear algebra calls like Matlab.  These are essentially programming considerations.
LazyMan FFT, Most DSP libraries have FFT libraries which have lengths that are a power of 2.  A common way to perform a FFT on sequences is to zero pad them so that have a length corresponding to a power of 2.  A library like FFTW is much more flexible in the lengths that can be selected but FFTW isn't available for every processor and there is a possibility that the zero pad approach may be faster.
Convolution by FFT (circular convolution  IFFT(FFT(x).$*$ FFT(Y)), where $.*$ denotes element by element by element multiplication ),  When one convolves a sequence of length M with a sequence of length N, the result is a sequence has length M+N-1,  The Element by Element multiplication requires that both sequences be the same length (prior to DFT, hence same after DFT), so the shorter is zero padded so that both sequences have the same length.  There are 2 variation of FFT convolution, one requires that both sequences be zero padded to at least M+N in length.
Zero padding to interpolate a DFT, usually by a whole factor, double, triple, ... The orthogonal bases of a DFT are of length N for a N point DFT, so we can use zero padding to calculate intermediate points of the DFT of a short sequence.  The resolution isn't increased, the magnitude plots are smoother.  One can zero pad in front or in back, the phase is different in the resulting DFT but the magnitudes are identical
Arbitrary Time shift of signal with a window.  In beamforming for instance, a number of sequences are produced by sensors and these sequences need to be time shifted prior to being added.  This can be achieved by zero padding each time sequence, performing a DFT, using the DFT time shift property, and performing the beamforming superposition in the frequency domain, and inverse transforming the delayed and summed frequency domain sequences.
Zero padding is often relatively cheap. A FFT algorithm can be modified so that zero padding avoids math operations relative to a full DFT acting on zeros.
Linearity and Time Invariance,  Zero padding is often used to take advantage of those properties, and are implicit in all the above bullets
While inserting zeros between samples isn't called zero padding, interleaving zeros (with filtering) can be used to interpolate in the time domain.

.


Help with IQ balance equations? I have developed a spectrum scope that plots signal vs. frequency starting from IQ signals.  I currently do not have a way to balance the IQ signals, so am looking at various ways to do so.  I have come across this article, which is reasonably easy for me to follow since I am fairly new to all DSP. The one question I have is concerning the notation of some of the equations:
For example, $a^2 = 2\langle t(i) * t(i)\rangle$. What does $\langle\cdot \rangle$ represent?  I can find no reference for these in mathematical symbol tables.

I wrote code in Java based on that article but it isn't working for me. Alpha (I-data scalar) calculates to be around $100$ and $cos(\omega)$ often calculates to an imaginary number. Did anyone have any luck with this or could point me to a better IQ Correction algorithm?

(To be clear: calculateMatrix, a, c, and d are all instance variables)

public void correctIQ(int[] data) {

    int length = data.length / 2;

    // ( I'(t) & Q'(t) )
    int[] iData = new int[length];
    int[] qData = new int[length];

    // Deinterleave IQ data
    for (int i = 0; i < length; i++) {
        iData[i] = data[i * 2];
        qData[i] = data[i * 2 + 1];
    }

    // ( I''(t) & Q''(t) )
    double[] qDataOffset = new double[length];
    double[] iDataOffset = new double[length];

    // Calculate DC offset
    double qOffset = mean(qData);
    double iOffset = mean(iData);

    // Apply DC offset
    for (int i = 0; i < length; i++) {
        iDataOffset[i] = iData[i] - iOffset;
        qDataOffset[i] = qData[i] - qOffset;
    }

    // Provided equations calculated on frequency change
    if(calculateMatrix) {
        double alpha = Math.sqrt(2.0 * multiplyMean(iDataOffset, iDataOffset));
        double sinw = (2.0 / alpha) * multiplyMean(iDataOffset, qDataOffset);
        double cosw = Math.sqrt(1.0 - (sinw * sinw));

        a = 1.0f / alpha;
        c = -sinw / (alpha * cosw);
        d = 1.0f / cosw;
    }

    // Replace into interleaved array
    for(int i=0;i<length;i++) {
        data[i*2]=(int)((iData[i]-iOffset)*a);
        data[i*2+1]=(int)((iData[i]-iOffset)*c+(qData[i]-qOffset)*d);
    }
}

double mean(int[] values) {
    long sum = 0L;
    for(int i=0;i<values.length;i++) {
        sum+=values[i];
    }
    return ((double)sum)/values.length;
}

double multiplyMean(double[] values1, double[] values2) {
    long sum = 0L;
    for(int i=0;i<values1.length;i++) {
        sum+=values1[i]*values2[i];
    }
    return ((double)sum)/values1.length;
}


What filter is optimal for noise cancellation when working with very low sample count I am working on sound processing android application project in which I have to input sound and remove noise as much as possible at run-time, somewhat like streaming. I am recording video at one place and play at a different place with very low latency.
I tried using Chebyshev filter but found that it didn't work if I had less frame available to process. There was clipping at the boundaries of the sample so when the video was played I could hear hissing like sound.

To demonstrate what actually happened I am attaching a Matlab code that simulates the situation output that I am getting.

Matlab code when run on a long sample give reduced noise output.

clc;clear;close all;
[sample_data, sample_rate] = audioread('save.binary.wav');
Fs = sample_rate;                                       % Sampling Frequency (Hz)
Fn = Fs/2;                                              % Nyquist Frequency (Hz)
Wp = 1000/Fn;                                           % Passband Frequency (Normalised)
Ws = 1010/Fn;                                           % Stopband Frequency (Normalised)
Rp =   1;                                               % Passband Ripple (dB)
Rs = 150;                                               % Stopband Ripple (dB)
[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs);                         % Filter Order
[z,p,k] = cheby2(n,Rs,Ws,'low');                        % Filter Design
[soslp,glp] = zp2sos(z,p,k);                            % Convert To Second-Order-Section For Stability
figure(3)
freqz(soslp, 2^16, Fs)                                  % Filter Bode Plot

filtered_sound = zeros(size(sample_data));
filtered_sound = filtfilt(soslp, glp, sample_data);
sound(filtered_sound, sample_rate)
audiowrite('save.binary.test.wav',filtered_sound,Fs)

While when running on small samples of the big file gives, hissing sound (hear simulated by using small samples in a loop).

clc;clear;close all;
[sample_data, sample_rate] = audioread('save.binary.wav');
Fs = sample_rate;                                       % Sampling Frequency (Hz)
Fn = Fs/2;                                              % Nyquist Frequency (Hz)
Wp = 1000/Fn;                                           % Passband Frequency (Normalised)
Ws = 1010/Fn;                                           % Stopband Frequency (Normalised)
Rp =   1;                                               % Passband Ripple (dB)
Rs = 150;                                               % Stopband Ripple (dB)
[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs);                         % Filter Order
[z,p,k] = cheby2(n,Rs,Ws,'low');                        % Filter Design
[soslp,glp] = zp2sos(z,p,k);                            % Convert To Second-Order-Section For Stability
figure(3)
freqz(soslp, 2^16, Fs)                                  % Filter Bode Plot

BufferSize = 960;
bufferCount = 830400/BufferSize;
filtered_sound = zeros(size(sample_data));

for i = 0:(bufferCount-1)
    filtered_sound(BufferSize*i+(1:BufferSize),:) = filtfilt(soslp, glp, sample_data(BufferSize*i+(1:BufferSize),:));
    sound(filtered_sound(BufferSize*i+(1:BufferSize),:), sample_rate)
end
sound(filtered_sound, sample_rate)
audiowrite('save.binary.test.wav',filtered_sound,Fs)

So what I want to know is which sort of filter is suitable for such an application in which frame length is small due to the requirement of low latency? I tried filters like Butterworth , elliptic, etc . but the results weren't any different.

filtfilt filters backwards in time, so the latency can be very large (infinite in theory), depending on the length of the filter's impulse response above your desired noise floor (or quantization floor).  Thus filtfilt is best used only for offline filtering.

If you use an ordinary forward-in-time IIR filter, you need to save the internal state variables of your IIR filter(s) between each short frame.

An alternative, for a FIR filter (or if an IIR impulse response tapers below some noise floor in a short enough number of samples), is to use some form of zero-padding plus overlap add/save to carry the impulse response forward.


Some questions about the intuition of the DWT Assuming a DWT of a signal of length 8 with Haar filter taps. At the lowest level, I end up with a3 and d3 both of length 1, d2 of length 2 and d1 of length 4 which is the same number of coefficients of the original signal and which I can plot on a dyadic grid.
In contrast to the WPT in the DWT there is only basis representation for a given decomposition level ([a3,d3,d2,d1] for the 3 level decomposition in the example above). From this sequence one can choose to recover the original signal or zero out some coefficients and then reconstruct for denoising purposes and the likes.
The dyadic grid representation is analogous to the CWT and so displays signal energy over a non-uniformly discretized 2-dimensional grid of time and scale/frequency.
I often see in books a representation in which the coefficients are just concatenated and plotted on a 1-dimensional frequency axis just next to the original signal axis.

I find this representation quite confusing since it suggests to directly compare both graphs although they are displayed over time (figure a) and frequency (figure b) respectively.
Now my questions:

What is the intuition behind the stacked frequency plot of the coefficients, i.e. what is the best way to read it? Going from right to left there is a decreasing time-resolution but increasing freq resolution. In reference to the signal above, I would have the 4 d1 coefficients on the most right followed by the 2 d2 coefficients and so on. Is the correct way of interpreting this to say, that in the spectrum of [pi/2T, pi/T] I have frequency information of the original signal at just 4 different points in time (due to downsampling with a factor of 2)  for [pi/4T, pi/2T] one has freq information at only 2 points in time but being narrowed down to half the frequency interval than before and so on?

Does this mean, that the signal representation I described above has now been transformed from a time-series (the original signal) to a pseudo tf-representation and only becomes a pure time-signal again by sending it again through a synthesis filter bank? Can both representations still be considered equivalent? If one would not have the axis description, how would it be possible to recognize, if one is represented the original signal or a decomposition?

As for a WPT there are now a large number of possible signal representations and one has to keep track of the chosen one, if a later reconstruction is desired. Limiting ourselves to a level basis of level 3 in the above example, what is now the meaning of the approx and detail coefficients in every freq subband? For the DWT one can argue, that d3 are the signal details at the coarsest scale that are necessary to go from a3 to a2 and so forth. So the details capture the missed-out fine-grained signal details that were neglected along the decomposition.
Is there a similar meaning to a WPT decomposition? I have a hard time transferring this intuition to the more general signal deconstruction of a WPT

In all these transforms, usually boundary problems arise whenever using a filter that has a larger number of taps than 2. So decomposing a length 8 signal with dB2 filters already yields more than 4 values both for a1 and d1 at the first decomposition level. If I were to plot the signal at this level, would I have to cut off the additionally introduced values somehow or would the signal at this level just naturally comprise then more values than the original time signal?

Thanks a lot for any help to deepen my understanding of this complex matter

First, one should be cautious about processing short signals like these. I am unsure about the length of 3 for a3 and d3. Now, the many questions:

I would not call them "frequency plot", but stacked subband plots. They are just illustrations of the behavior of the wavelet coefficients. For each subband, different scalings may be used. Another rendering is stacking the intensity of coefficients in a image. For the interpretation of the points: the coefficients in each subband at located in time and scale, so the interpretation in $[\pi/2T, \pi/T]$, $[\pi/4T, \pi/2T]$ is only approximative (because of aliasing, subsampling, etc.).
No, one is concatenating chunks of coefficients in different subbands and at different rates
Yes, subbands of a given WPT can generally be interpreted in dyadic portions of the time-scale axis
Haar wavelets (2-tap filters) have the rare property (among standard wavelet filters) they don't overlap at a given scale. For longer filters with 2-band wavelets, those filters overlap, and yield troubles at the extremities of a finite-length signal. There is a lot of works on non-expensive expansions, relaying on the symmetry/antisymmetry of the wavelet.


Binarization of image for varying background and foreground colors for OCR I am trying to develop an automated OCR system for images obtained from the internet. The images mostly consist of computer generated text imposed on a background, so the OCR should produce near perfect results. However, the texts vary in color, and the background its imposed on can be of different colors too or even photos. The OCR tool often generates garbage when the colors are different enough, so I have to binarize the image properly, cleaning up the text content.

I found a promising article here, and implemented the algorithm using Python OpenCV. However, it seems to work poorly for most computer generated text. I also tried some local binarization algorithms, but they also fail to separate most text colors, like red on an orange background, the details being lost once converted to grayscale.

Is there an alternative image processing algorithm for this task ?

This problem is well studied. It is termed scene text recognition. Of course, the state-of-the-art methods involve deep learning. In particular, I like the works, which use a bi-directional LSTM together with features extracted from Convolutional Neural Networks. This allows for end-to-end training and frees the segmentation or binarization requirement. In particular, I like the work:

  Reading Scene Text in Deep Convolutional Sequences -
  Pan He, Weilin Huang, Yu Qiao1, Chen Change Loy2, Xiaoou Tang2
  AAAI 2016


Difference between frame rate and sampling rate? Below are the two different methods of generating time to create an audio sine wave. While framerate is being used as steps in one method, the inverse of sampling rate is being used as steps in the other method. Where is the difference?

Method 1

# Create a sound
framerate = 44100
time = np.linspace(0,5,framerate*5)         #Multiplying framerate

signal = np.sin(2*np.pi*220*time) + np.sin(2*np.pi*224*time)

----------------------------------

Method 2

#Simulation parameters
srate = 500 #Sampling Rate in Hz
time = np.arange(0.0, 2.0, 1.0/srate)       #Dividing Sampling rate

#Sine Wave Parameters
freq = 3         #Frequency in Hz
ampl = 2         #Amplitude in a.u.
phas = np.pi/3   #Phase in Radians

#Generate the Sine Wave
sine_wave = ampl*np.sin(2*np.pi*freq*time + phas )

One method uses np.linspace() and the other one np.arrange() to create the time domain vector.

I suggest reading the documentation

https://docs.scipy.org/doc/numpy/reference/generated/numpy.arange.html
https://docs.scipy.org/doc/numpy/reference/generated/numpy.linspace.html

The first method has a small error, so the frequency will not be 100% correct.

Update

Common error with using linspace: Let's say you want three samples that are 1 apart. linspace(0,3,3) will give you [0 1.5 3], so the spacing is bigger than you want. The correct way is to adjust the upper limit by subtracting 1 interval (or 1/samplerate), so it should be linspace(0,2,3) or

linspace(0,5-1/framerate,framerate*5)


Phase response for conjugate zeros If a second order system has 2 poles/zeros that are conjugate symmetric, how does this affect the phase response? I know that if there are 4 zeros/poles that are conjugate reciprocals, then it is a zero phase system.

If you have conjugate symmetric poles and zeros then the coefficients of the system (i.e., the filter coefficients, or, equivalently, the coefficients of the corresponding difference equation) are real-valued. E.g., if the zeros are $z_0$ and $z_0^*$ then the numerator polynomial of the transfer function $H(z)$ is (up to a scaling factor)

$$N(z)=(z-z_0)(z-z_0^*)=z^2-(z_0+z_0^*)z+|z_0|^2=z^2-2\text{Re}(z_0)z+|z_0|^2\tag{1}$$

Note that the coefficients in $(1)$ are all real-valued. The same is true for the denominator polynomial. So the frequency response (i.e., the transfer function evaluated at $z=e^{j\omega}$) can be written as

$$H(e^{j\omega})=\frac{b_0e^{2j\omega}+b_1e^{j\omega}+b_2}{e^{2j\omega}+a_1e^{j\omega}+a_2}\tag{2}$$

with all coefficients $a_i$ and $b_i$ being real-valued. Consequently, the following equation must hold:

$$H(e^{j\omega})=H^*(e^{-j\omega})\tag{3}$$

If you write $H(e^{j\omega})$ in polar form $H(e^{j\omega})$=$M(\omega)e^{j\phi(\omega)}$ with magnitude $M(\omega)$ and phase $\phi(\omega)$ you get

$$M(\omega)e^{j\phi(\omega)}=M(-\omega)e^{-j\phi(-\omega)}\tag{4}$$

From $(4)$ it follows that for conjugate symmetric poles and zeros the magnitude of the frequency response must be an even function of frequency, and the phase must be an odd function of frequency, i.e., $M(\omega)=M(-\omega)$ and $\phi(\omega)=-\phi(-\omega)$.


Audio fingerprint questions I am developing an application that requires audio fingerprints. I have been reading a lot of articles and PDFs, now i think i have gotten myself confused. Based on my present understanding, I have some questions

After decoding the audio to its raw format, resampling and extracting the mono channel. Do i convert the bytes to integers or floats?
What is convolution and why is it necessary to convolve the audio samples
What is a window and its length
FFT transforms a signal from time domain to frequency domain. Am i correct? If yes, does the frequency components of an audio determine the contents/sound/noise/volume?
What operation can be carried out  on two audios of same content but different bitrate to normalize the bitrate and get the data that they share
What are low pass and high pass filters and how are the derived. What are their uses in relation to audio fingerprint operations
In the text Computer vision for music identification. After convolving the signals with a low pass filter and extracting every 8th sample. Then, a short term fourier transform with a window size of 2048 samples with successive windows offset by 64 samples. Also divide the power between 300Hz and 2000Hz in 33 logarithmically spaced bands.

Will these operations be applied to each of the samples retrieved after the convolution process? Can i get a simpler explanation?

32 learned filters and thresholds are applied to get a 32-bit descriptor for every time step (11.6 ms) of signal. This series of signals are known as the signature.

what are these learned filters and thresholds?

Wow :), that's a lot of questions. I am sure a lot of beginners will find the answers helpful.

Thanks

There are rare situations when processing is done with integers (embedded systems with no floating point units), but usually, audio signal processing is done with floats.
(Handwaving) Convolution is the operation through which a filter is applied to a signal. If the authors of the paper you're looking at found that applying a specific filter to the data (for example to remove noise, any information not processed by the human auditory system, etc.) yielded good results, then you have to trust them on this!
Since audio signals are not stationary, they are processed in short blocks, referred to as "windows" or "frames". The analysis/processing is thus done on short overlapping segments. For example you will process segment 0.00s:0.02s ; then segment 0.01s:0.03s ; then 0.02s:0.04s etc. In this example, the window size is 20ms and the overlap 50%. There is no "universal" window size for audio applications, look at what is used in the paper describing the technique you're implementing.
Not really, but going from the time domain to frequency domain is "one step closer" to what is ultimately happening in the human auditory system.
No operation could do that. If there was one, we would send all our files with the lowest quality and apply this magical operation to recover the highest quality. What this means is that the features used for fingerprinting must be robust enough to capture only the overall "shape" of the sound rather than focus on details - which would be lost by an increased amount of MP3 encoding.
Look at a signal processing textbook. These are very basic signal processing operations and can be used for many things inside or outside the context of a fingerprinting system - reducing sample rate, smoothing temporal sequences, emphasizing details, removing noise...
Let's say you have 800000 samples. You apply the convolution, you get 800000 samples out of it. Now you extract every 8th sample, you get 100000 samples. Then you apply your FFT on samples 0..2047 to get your 33 coefficients ; then on samples 64..2111 etc... so in the end you have a list of roughly 100000 / 64 = 1562 vectors of 33 coefficients.
(Handwaving). I am assuming you are referring to the Ke et al system (later refined in the works of Baluja et al), which is itself inspired by the Viola-Jones face detector. Once you have converted the signal into a sequence of vectors of 33 coefficients, you compute weighted sums of the coefficients in adjacent bins or along the time axis (a filter), and you compare the result to a threshold. This gives you a binary "feature". There are 32 such binary features, each of them computed with a different filter and a different threshold. They are "learned", because they were obtained the following ways: the authors built a database with millions of pairs of songs, some of them identical (but maybe at difference bitrate, volume, noise level etc.), some of them different. Using machine learning techniques (boosting on decision stumps), they evaluated which filters were the best at binning together identical content and separating unrelated content, and which threshold were being used for that. You can find the list of those filters/thresholds in the source code of an open-source implementation of lastfm's Ke-style audio fingerprinter.


Implementation of Fourier Domain Denoising with Hard Threshold I just tried the Fourier denoising method with a hard threshold and my code is as follow:

F = imnoise(phantom(128),'gaussian');
B = fft(F);
s = 0.3;
for i = 1:128
for j = 1:128
    if abs(B(i,j))<s
        B(i,j)=0;
    end
end
end

C=ifft(B);
subplot(121)
imshow(F,[])
subplot(122)
imshow(C,[])

Left is the phantom after adding noise and the right is the denoised image from my code.

The Fourier denoising hard threshold method just uses threshold value to keep high frequency coefficients and the coefficients below the threshold to be 0.

In my code, I just added gaussian noise with default mean 0 and variance 0.01 to my phantom image and then I set the threshold to be 3 times of the standard deviation of the noise, which is 0.3 (this is some good practical threshold I looked at from literature). But why the result is not good? It seems it is still very noisy. Even I change my threshold the result is not good. Why? And I want to my denoised image as close to the original phantom.

Thanks in advance!

Thresholding in the Fourier domain is an archaic method sometimes called spectral subtraction, used for background noise removal in speech. Bad results in image processing can be due to several factors:

Misinterpretation: you say "The Fourier denoising hard threshold method just uses threshold value to keep high frequency coefficients". Not quite. You keep "high amplitude coefficients", no matter low or high frequency.
Test data: the Shepp-Logan model is known to be quite specific. Bad results on it do not imply bad results on real data, yet:
Transform choice: Fourier is very good for close-to-stationnary data. Standard images with smooth parts inside contours are not well represented in the Fourier domain, at least not concentrated enough for acurate data compaction and reconstruction with few (high) coeffcients. Derivatives, or multiscale directional derivatives like wavelets can help better.
Implementation: Fourier implementation via FFT suffers boundary artifacts
coefficients, not talking about the 1D implementation you seem to use (prefer fft2 here)
Coefficient selection: hard and scalar thresholding typically yields more important visual artifacts than "smoother" thresholding and vector or group thresholding.
Formulation: you can generally get better results when you do not use a "parameter-tuned algorithm", but prefer a variational formulation based on your priors, including noise statistics.

So you have plenty of room to improve your results. Two references for multivariate image denoising with automated threshold (for Gaussian noise):

A nonlinear Stein based estimator for multichannel image denoising, 2008
Wavelet transform for the denoising of multivariate images, 2010


Biquad Filter Transfer Function Frequency Response $H(z) = G\cdot\displaystyle\frac{2b_0\left(1-z^{-2}\right)}{1-2a_1z^{-1}-2a_2z^{-2}} + 1$

I have a data sheet with the above transfer function, I need to calculate its frequency response in code. I have the functions to calculate $b_0$, $a_1$, $a_2$ and $G$, but it has been a loooooong time since I have handled $H(z)$ transfer functions and I can't for the life of me remember how to work it out.

I have googled around and it mentions jw this and that and I have tried everything Ive found but I'm not getting a frequency response anything like what is expected.

I have a sampling rate of $48000\textrm{ Hz}$ and want to calculate the frequency response from $0-24000\textrm{ Hz}$ in $1\textrm{ Hz}$ steps.

Any help would be very much appreciated.

To get the frequency response, replace $z$ in $H(z)$ by $e^{j2\pi f/f_s}=\cos(2\pi f/f_s)+j\sin(2\pi f/f_s)$ where $f_s$ is the sampling frequency in Hz and $f$ is the desired frequency in Hz. This will give you the value of the frequency response at frequency $f$. Note that the frequency response is a complex function, so you might want to look at its magnitude.


Is it possible to open up this particular convolution equation? I have three vectors, $x[n]$, $w[n]$, and $z[n]$. I need to compute:

$y[n] = z[n] * \Big(x[n] \cdot w[n] \Big)$.

This is easy to do. However, is there a way I can do this by 'opening' up the bracket?

For example, let us say that I am already given $z[n] * x[n]$. If this was the starting point, what could I do with this and with the three already given vectors to get $y[n]$? The wiki on convolution properties does not seem to have a case for a scenario like this.

Here is a re-write of a previously mistaken answer:

Let us write $y[n] = z[n] * \Big(x[n] \cdot w[n] \Big) = \displaystyle\sum_{m=-M}^{M} {z[n-m] \cdot (x[m] \cdot w[m])}$ and rearrange as $y[n] = \displaystyle\sum_{m=-M}^{M} {(z[n-m] \cdot x[m]) \cdot (w[m])}$. After this point, we cannot take the terms out of the brackets.


importance of using histogram equalization I'm reading opencv tutorials and I'm diving into histogram equalization. i have looked in wikipedia, there is a nice example that sums up exactly the problem:

original:

equalized:

but for getting this result i would take a different approach:

find the minimum and maximum in the original
normalize (remap) everything upon it.

no histogram, no  cumulative distribution function. for sure a more stupid approach but i can't see the difference.. why use histogram equalization? can someone can help me in make the reason out?

The difference is that with your method, if you imagine the histogram, you're simply going to stretch it to span from 0 to 255, but its shape will be preserved. Histogram equalization not only stretches your histogram, but also tries to make it flat, so that you get an approximately even distribution of pixels of every shade of gray.

In terms of why one is better for certain applications than others, that's application specific.

Edit:

Here are some examples from Bruzed:

Contrast stretching (your method):

Histogram equalization:

As you can see, the bulk of the pixels in the original image was gray, represented by a large peak in the middle. When you do contrast stretching, the peak is still there, even though your darkest pixel is now black and you brightest pixel is now white. By contrast (pun intended), using histogram stretching you get a much flatter histogram response. This actually increases the overall contrast of the image.


Filtering Overlapping Frequencies in sound file I have been given a sound file of a plane passing over a rain forest filled with birds. I am supposed to filter out the sound of the plane as it flies over.  I've accomplished this with various types of filters in MATLAB, but I always run into one problem. I can either cut out all of the plane and lose some of the background rain forest noise, or leave (almost) all of the rainforest noise and leave a decent amount of airplane noise. This is the frequency spectrum of the original sound file:

(Sorry for the crappy image.) Note the different colors represent left and right channels. From the plot, the majority of the plane's frequency is in the 0-1000 Hz range. The plane is clearly heard however until I use a high pass filter to remove everything from around 1700 Hz to 0 like so:

I suspect that the plane is causing that spike around 1700 Hz but evidently the rain forest is still producing lower frequency noise than that.

What would be the best method to remove that spike while keeping all that extra rain forest frequency?

EDIT:

Here is the plot of the time domain. Note how the plane gradually gets louder than fades away.

The sampling frequency is 22050 samples/sec.

And the spectrogram.

In my opinion, the best approach for this kind of spectrally-overlapping noise is to use blind source separation techniques, like independent component analysis (ICA) or time-frequency masking, or sparse source decomposition.  Look up the work of Emmanuel Vincent, Kostas Kokkinakis, or Philipos Loizou.  These methods work best when some statistics of the noise are time-independent, and the SNR is small or negative, as appears to be the case in your example.  Vincent has some MATLAB code available on the web.


What does structure matrix do can anybody explain to me what does structure matrix exactly do in image processing? What happens to the pixels?

By structure matrix I mean:

$$S_w[p] = \begin{pmatrix} \sum_{r}w(r)[I_x(p-r)]^2 & \sum_{r}w(r)[I_x(p-r) I_y(p-r)]  \\ \sum_{r}w(r)[I_x(p-r) I_y(p-r)] & \sum_{r}w(r)[I_y(p-r)]^2 \end{pmatrix}$$

At each pixel $p$, this $2\times 2$ matrix summarizes the local behavior  of the image (structure) inside a neighborhood of pixel $p$ (a neighborhood weighted by function $w$.

From this matrix, different quantities can be computed, that provide interesting attributes for image analysis. Structure tensors, or second-order matrices, have been used for image segmentation, especially with geodesic or level-set methods, as they can provide insights to image textures or orientations.


Question about ramp filter used in filtered backprojection Question is this.
First, a ramp filter (in frequency domain) is defined by $H(Q)=|Q|$.
What are the responses of a ramp filter to
(1) a constant function $f(r)=c$ and (2) a sinusoid function $f(r)=\sin(wr)$?
What does the response mean? Following is my work.

My work:

First, take fourier transform of a function $f(r)=c$. It is
$\int_{-\infty}^{\infty}f(r)e^{-2i\pi rQ}dr=c\delta(Q)$. Then multiply ramp filter and take inverse fourier transform. It is $\int_{-\infty}^{\infty}c\delta(Q)|Q|e^{2i\pi irQ}dQ=0$??
Similarly, $\int_{-\infty}^{\infty}\sin(wr)e^{-2i\pi rQ}dr=\frac{\delta(Q-w/2\pi)-\delta(Q+w/2\pi)}{2i}$. So Applying the ramp filter and i.f.t gives
$\frac{w(e^{iwr}-e^{-iwr})}{4i\pi}=\frac{w\sin(wr)}{2\pi}$.

It this right?

To see if your math is correct, it is useful to first understand what is, in general, the effect of a filter on a signal, and then see if you can predict what the theoretical result should look like.

If a filter has frequency response $H(Q)$, this means that its response to an input $e^{j2\pi Q_0 r}$ is the signal $H(Q_0)e^{j2\pi Q_0 r}$. In other words, a sinusoidal input of frequency $Q_0$ produces an output of the same frequency, but with amplitude $|H(Q_0)|$ and phase $\angle H(Q_0)$.

In your first question, the input has frequency $Q_0=0$. The filter's response at that frequency is $|Q_0|=0$. Then, the filter's output will be 0: frequency $Q_0=0$ is completely absorbed by the filter and it does not appear at the output.

In your second question, the input has frequency $Q_0=w/2\pi$. The filter's response at that frequency is $|Q_0|=w/2\pi$. The output, then, should be $\frac{w}{2\pi}\sin(wr)$.

As you can see, following this line of reasoning we obtain the same results as you did by just solving the equations. This corroboration should give you confidence that your results are correct, and most importantly, you can understand what the filter is doing physically.