Deconvolution of an Image Acquired by a Square Uniform Detector So, I acquired some images by scanning a radiation source with a square detector like in the following gif.

Where the dashed grid represents reality, the 3x3 square my detector, and the 4x4 my acquired data.
Obviously, because I'm moving the detector in steps smaller than its size, my image is convoluted. Assuming my detector has uniform sensitivity, how do I go about deconvolving it?

The point spread function should be simple, but I can't seem to find examples anywhere.

(I can work with Python or Matlab, so feel free to give any examples you want)

The data I'm working with consists of gaussian-ish peaks surrounded by noise like so:

Your model is exactly a Convolution with Uniform Kernel where the output is what is called the Valid Part of the Convolution.

In MATLAB lingo it will be using conv2(mA, mK, 'valid').

So the way to solve it will be using a matrix form of the convolution and solving the linear system of equations.

Let's use the Lenna Image as input (Size was reduced for faster calculations):

We have a uniform kernel for the sensor model.
The output of the convolution with uniform kernel is given by:

The output from the sensor is both blurred and smaller (Less 2 rows and 2 columns as it is 3x3 kernel) just as in your model. This is the model of Valid Convolution.

In Matrix form what we have is:

$$ \boldsymbol{b} = K \boldsymbol{a} $$

Where $ \boldsymbol{b} $ is the column stack vector of the output image, $ \boldsymbol{a} $ is the column stack vector of the input image and $ K $ is the convolution operator (Valid Convolution) in matrix form. In the code it is done in the function CreateConvMtx2D().

So, now all we need is to restore the image by solving the Matrix Equation.
Yet the issue is the equation is Underdetermined System and the matrix has high condition number which suggest not to solve this equation directly.

The solution is to use some kind of regularization of the least squares form of the problem:

$$ \arg \min_{\boldsymbol{a}} \frac{1}{2} {\left\| K \boldsymbol{a} - \boldsymbol{b}  \right\|}_{2}^{2} + \lambda r \left( \boldsymbol{a} \right) $$

Where $ r \left( \boldsymbol{a} \right) $ is the regularization term.  In the optimal case the regularization should match the prior knowledge on the problem. For instance, in Image Processing we can assume a Piece Wise Smooth / Constant Model which matches the Total Variation regularization.

Since we have no knowledge here, we will use the classic regularization to handle the Condition Number - Tikhonov Regularization:

$$ \arg \min_{\boldsymbol{a}} \frac{1}{2} {\left\| K \boldsymbol{a} - \boldsymbol{b}  \right\|}_{2}^{2} + \frac{\lambda}{2}  {\left\| \boldsymbol{a} \right\|}_{2}^{2} =
{\left( {K}^{T} K + \lambda I \right)}^{-1} {K}^{T} \boldsymbol{b} $$

The output is given by (For $ \lambda = 0.005 $):

We can see that near the edge we have some artifacts which are due to the fact the system is Underdetermined and we have less equations to describe those pixels.
One can use the $ \lambda $ parameter to balance between how sharp the output is (Yet with artifacts) to how smooth it is, basically governing the level inversion of the system.

I advise playing with the parameter to find the best balance for your case but more than that, find a better regularization. Since the information you're after looks smooth you can use something in that direction.

The full MATLAB code is available on my StackExchange Signal Processing Q63449 GitHub Repository.

Enjoy...


Relation between two k-spaces phase-frequency and spatial frequencies in When I see MRI explained, two types of 2D k-space images seem to be described as if they were the same.

Axes are the two spatial frequencies. This images is directly fourier-transformed into the final image.
Axes are phase and frequency (temporal, I assume). The data read through the RF coil is put directly into this matrix.

What is the relation between 1 and 2? Are they the same?

Example with spatial frequencies:

Example with phase and frequency, which seems to be the "raw data":

"k-phase" and "k-readout", which I think is identical to phase and frequency.

In MRI the term $k$-space usually means that you have spatial frequencies in a matrix. The Fourier Transformation of this matrix is the image.

The terms readout and frequency correspond to each other; they mean the same. The readout direction is a logical direction, along which a gradient is applied during signal acquisition, i.e. during readout. This gradient field causes a frequency spread of the excited magnetization along the direction of this gradient, hence the second name readout for this direction and for the corresponding signal axis in $k$-space.

The term phase is associated with another direction, usually perpendicular to the readout/frequency direction. It is called so, because a gradient field is applied before the signal is read out. This causes a frequency shift, too, but the gradient field is switched off before the read out of the signal, hence in this direction there is only one frequency during read out. What is left from this frequency spread is a phase shift of the magnetization vectors along that direction. The trick here is to repeat the measurement with different field strengths of that phase encoding gradient - that way, this direction is fully encoded.

Of course you are right: The coil that measures the signal does measure it in a time-domain. Keep in mind, that you additionally switch magnetic gradient fields. These fields act upon the magnetization and hence alter the measured signal. Since the temporal behavior of those gradient fields is known and chosen carefully (to fulfill certain criteria such as the Sampling Theorem, etc.), there is a simple relation between the time $t$ of the measured signal and the corresponding $k$-space position: $\vec{k}(t) = \frac{\gamma}{2\pi}\int_0^t \vec{G}(\tau)\,\textrm{d}\tau$.

Note: Even though the phase and readout gradients look completely different in sequence diagrams and appear to be completely different encoding techniques - they are completely the same! Even for physicists who are new to the field of MRI this bears some trouble initially in understanding the image measurement process.


Help in interpreting the auto correlation graph I want to check if a time series is (a) random  (b) independent. For these I am using the autocorrelation (AC). Autocorrelation
refers to the correlation of a time series with its own past and future values.

If random variables are independent then there should be no correlation between them. The reverse is not true.

I have two different vectors as observations obtained from sensor measurement that are a sequence of samples and I treat this as a time series -denoted by variables phi and phit respectively.

Question 1: I cannot interpret the graph, what the spike indicates. Please help
Question 2: For time series generated from linear system in order to check for randomness, what  is the test to perform and what is the value to expect??
Question 3: For time series obtained from linear system in order to check for independence between samples, what test to perform and what is the value to expect?
Question 4: The given time series is obtained from nonlinear dynamical system that is chaotic.  If correlation test is not valid, then what tests should I perform to check for (a) randomness (b) independence?

In my opinion this question is too broad. The first point is that you want to see if a given time series is random or not. I have doubt on how we can define this term here. Also you don't mention any distribution. Normally, we are interested to see if a time-series belongs to a specific distribution or not, for which there are standard techniques.

Regarding dependency, you should note that correlation only captures linear dependencies. You can have random variables that are highly dependent but have very low correlation. So it is important that what specific type of dependency you are looking for.

Edit: That said, referring to the ACF a peak at lag zero is visible. The reason is that at zero lag the time-series shows maximum similarity which is not surprising since it is obviously similar to itself! Note that the amplitude of this peak for a general complex signal $x[n]$ is $$R_{XX}[0]=\sum{x[n]x^*[n]}=\sum{|x[n]|^2}$$ which is the signal's power. You can easily show (theoretically) that $$R_{XX}[k]\le R_{XX}[0], \,\forall k$$

In addition, you should note that the autocorrelation of a white process (i.e. a process with uncorrelated random variables) is a pure delta function. So closer your autocorrelation to a delta means it is more uncorrelated, which is a good thing for a PRNG.

Regarding, your question about mutual information and entropy, yes both of them are used. Entropy is a well known measure for randomness. I cannot give you more details in here but can suggest you looking at this or similar resources.


How do you convert angular frequency to $\textrm{Hz}$ from MATLAB's $\tt freqz$ function? Let's say I have used MATLAB's butter function to generate a 3rd order, low-pass, Buttersworth, digital filter. I have a sample rate of $2000\textrm{ samples/second}$, and I want the cutoff frequency to correspond to an analog frequency of $700\textrm{ Hz}$. Because MATLAB defines its normalized frequencies as percentages of the Nyquist frequency, the $700\textrm{ Hz}$ would turn into $700/1000 = 0.7$. I feed these parameters into the butter function and outcomes my filter coefficients.

I then use the freqz function on these coefficients to obtain the [h, w] vectors. h corresponds to the filter output at a specific frequency, and w is the angular frequency.

I'm having a difficult time understanding the units of the angular frequency. I would typically think of $\omega$ of having units of $\rm rad/s$, but i don't believe this is correct due to us having normalized the frequency earlier.

Would anyone mind explaining to me the units of w outputs by freqz?
And specifically how I would be able to scale that back to an analog frequency in $\textrm{Hz}$?

After a little hard thought, I believe I stumbled across the answer. In order to know how to convert, you really need to pay attention to units.

When you use the 'butter' function, it is necessary that you normalize the frequency. By normalizing the frequency, you are essentially removing the time dependency and are only focused on relative values. Commonly in DSP, this frequency is achieved by dividing the analog frequencies by the sampling rate. In this, the units of the normalized frequency are cycles/sample. This was achieved by dividing the analog frequency in hz (cycles/second) by the sample rate (samples/second). Dividing these yields the units of cycles/sample. Notice now that there is no units of seconds, we have thus removed the time dependency.

However, MATLAB is different in the sense that it likes to normalize the analog frequencies in terms of the Nyquist rate as opposed to the sample rate. The Nyquist rate is 1/2 the sample rate. The reason for this is due sampling theorem which states that the sampling rate must be twice that of the true frequency in order to capture the true frequency. For analog frequencies greater than 1/2 the sampling rate, they will be aliased onto frequencies less than 1/2 the sampling rate. Again, just refer to the sampling theorem.

So now, using MATLAB's convention, the real units of the normalized frequency is (1/2 cycles)/sample. Fine and dandy.

After using the 'butter' function to obtain the filter coefficients, we feed them into the 'freqz' function which will create an h and w vector. h and w are the frequency response and angular frequencies respectively. Assuming you know about frequency responses/bode plots, the h vector should make sense.

However, the real question is regarding the w vector. MATLAB states that the w vector is the angular frequency vector. But what are the units of this angular frequency? Well typically when you think of angular frequency in terms of radians, you think of radians/second. However, since we're working with normalized frequencies with removed time dependencies, the units of normalized angular frequencies is radians/sample. If you look at the range of the w vector, you will notice that is goes from 0 radians/ample to pi radians/sample. This make sense due to the fact that pi radians/sample is equivalent to (1/2 cycles)/sample which covers the entire normalized frequency range up until the normalized Nyquist frequency.

So we're now at the point where we have our filter response magnitude 'h', but it's been evaluated at the weird units of radians/sample. We want to get back to our analog frequency equivalent which has units of Hz (Hz = cycles/second).

It's now simply a game of unit analysis of converting radians/sample -> cycles/second.

That's achieved by doing the following operation to the w vector:

w*Fs/(2*pi) -> units of Hz.

TL;DR
freqz computes the digital filter response in terms of radians/sample. We want to convert this to Hz (which is cycle/second). To do so, multiply by Fs/(2*pi).


Why does FFT output need to conjugate processing? I came across some matlab code as the below,

fftout   = fft(adc_in, 1024);
fftout_p = fftout.*conj(fftout);

I think fftout is formatted as complex form. from here I'm ok.

Next stage that signal is multiplicated by itself which is conjugated.

I think it means that it want to get the real number not complex number.

From here, I'm not sure and want to know that practically what does Imaginary part affect to real signal processing?

Can we ignore like above that?

In real practically engineer field, What do conjugating and Imaginary part mean and affect to the real engineer field?

The real part represents cosine correlation, and the imaginary coefficient represent the sine correlation.  Cosines are even or symmetric functions in the DFT aperture, and sines are odd functions (or anti-symmetric) in the DFT aperture.  So you can only ignore the imaginary coefficients for exactly symmetric input.

Conjugation and multiplication is “just” a mathematical trick (identity) to compute sin^2 + cos^2 which is the hypotenuse squared (Pythagoras).


harmonic algorithm for extending compressed bandwidth of audio file I just found this site and feel at home :)

I've been playin around with python and some compression ideas and was wondering if there's anything like what I am trying to figure out...

I've noticed that mp4's uploaded to instagram go from ~0 to 13100Hz and on Youtube from ~0 to ~14000Hz... which is OK, but what if there was an algorithm to extend or at least soften this highcut, by taking the amplitudes of the frequencies within the respective bandwidths and applying a harmonic function, (possibly with averaging/weights), to add in implied values to the currently non-existant values (i.e. the and additional bandwidth from 14,000 to 20,000Hz would be harmonically synthesized and still have realistic values.) It could be a plug-in or something.

Anyways, I'm going to check out the board some, and hope to hear some smart replies! thanks!

there is also this answer in this comp.dsp thread.

high-pass filter your input signal (keep the unmolested input around for later).  the cutoff frequency should be adjustable.  i would think a 2nd-order HPF with a Q approx= 0.707 (Butterworth) would be good enough.
square that HPF output,
add (maybe with an adjustable gain coefficient) that squared HPF back to the original input and output that audio sample.

see if that brightens up your day.


Is there a way to do a DFT with low latency and fine frequency resolution at low frequencies? I am learning about the DFT in order to design a fundamental pitch detector for the human singing voice. Using the cepstrum technique.

I am concerned about the latency for low frequency detection at <100Hz. Is it always necessary for the DFT algorithm to work with sample length that contains an entire wavelength of the lowest frequency being detected, or can the window be shorter than the wavelength? For example if the singer sings <100 Hz then there is the possibility of the >10ms latency. This latency could get worse with further processing.

I am also concerned about frequency resolution at low frequencies. The standard DFT is inherently linear this causes poorer pitch discrimination at lower frequencies, If frequency index 1 was 50Hz, then the next bins would be 100Hz and 150Hz. The steps are clearly too big for pitch detection. If the index started at a lower frequency to get finer graduations, say 1 Hz steps, then the latency becomes 1 second, if what I said in the previous paragraph is true. Btw 1Hz steps is still pretty poor for pitch detection.

Are there more sophisticated DFT algorithms that address these issues?

I know products like Antares Auto-Tune have pretty fast detection, so this stuff is possible.

The answer is no.  You've explained some of the reasons why a DFT alone is not a good first step when trying to do low latency pitch estimation on low pitched human voices.

Around the first DFT bin, the effects of the window and the complex conjugate image can interfere with just spectral frequency estimation, which a pitch of that frequency may even be missing.

But it may be possible to do pitch estimation with better resolution and latency than a bare DFT by using other methods (time domain, parametric, or composite, look for many other similar Q&As about pitch estimation on this site).


Windowing of signal before FFT I am facing an interesting problem: I have a signal in the time domain. In order to avoid discontinuities at the edges of the signal --- which would lead to artifacts when calculating the FFT --- I apply a window function, like e.g. hann or hamming. Then I perform the FFT of this windowed signal.
To preserve the energy of the signal after the FFT, I calibrate the FFT of the windowed signal with the RMS of the window function. However, the strange thing is that the energy (=area under the curve) of the windowed signal in the frequency domain is substantially different (about a factor of 2) compared to the FFT of the unwindowed signal.
Do you have an idea where this issue could originate from?
Thanks a lot in advance!

So you think that the energy of the windowed signal in the frequency domain equals to the original signal multiplied by a certain value related to RMS of the window function, which is wrong. According to Parseval's theorem, the total energy of a signal can be calculated by summing power-per-sample across time or spectral power across frequency:
$$
\sum_{n=0}^{N-1} |x[n]|^2 = \frac{1}{N}\sum_{k=0}^{N-1}|X[k]|^2
$$
Hence I'm going to prove your statement wrong in the time domain.
The energy of a windowed signal $w[n]x[n]$ is
$$
\sum_{n=0}^{N-1} |w[n]x[n]|^2 = \sum_{n=0}^{N-1} |w[n]|^2|x[n]|^2 \neq C\sum_{n=0}^{N-1} |x[n]|^2
$$
where $C$ is a constant. The above inequality holds true as long as the window function $w[n]$ is not a rectangular window.


Estimating Convolution Input Under the Assumption of Sparsity and Constant Non Zero Values Using Compressive Sensing Approach I was wondering about if there is compressive sensing algorithm to estimate the sparse vector where the number of non-zeros values and amplitude of every non-zeros value are known. For example, assume we have a vector $x$ whose length is $N$x$1$ with only $N/2$ equal non-zeros known values but unknown location of those values. That vector  is convolved with a channel $h$ resulting a vector $y$. it means that:

$y = h⊗x$, where $⊗$ is the convolution operation.

Is it possible to use compressive sensing to estimate the locations of non-zeros values in $x$ based on $y$.

Here is the code of an example where the length of vector $x$ is $32$ and channel $h$ = 16:

clear all; clc;

%%%% Build the sparse vector
X = hadamard(32);
X2 = randi([1 length(X)-1], 1);

x = X(1,:) + X(X2+1,:);         % Here the built sparse vector
x = x / max(x);                 % To make the sparse vector either one or zeros

h = randn(1,16);                % channel
y = conv(x,h);
y = y(1:end-length(h)+1);       % To remove the delay of convolution

Thank you

Basically your problem is called Blind Deconvolution.
It means we want to estimate both the operator and the input given the output.

You model is Linear Time Invariant Operator so we have LTI Blind Deconvolution.
In general blind deconvolution is ill poised problem.
So we need to make assumptions about the model. The more assumptions the better the chance to solve this really hard problem.

What do we have in your case:

The input signal is sparse.
The input signal has 2 values, either zero or other known value.

What's missing is some assumptions on the operator $ h $.

Deconvolution in Image Processing

The field which pushes the deconvolution problem farther and farther is mostly the image processing field.
There are many models of real world images and convolution kernels.

Let's talk about the most common for each:

In most cases the convolution kernel is assumed to be LPF with its sum of coefficients equal 1 and each coefficient is non negative.
In most cases the image is assumed to be "Piece Wise Smooth. Enforcing it using the Total Variation Model which basically says the Gradients are distributed according to Laplace Distribution.

With those 2 models we can model the problem as:

$$\begin{aligned}
\arg \min_{h, x} \quad & \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} + \lambda \operatorname{TV} \left( x \right) \\
\text{subject to} \quad & \sum h = 1 \\
& {h}_{i} \geq 0 \\
\end{aligned}$$

As can be seen this is highly non convex problem. The method used to solve it is by splitting methods.

So we solve it by iterations:

We set $ {h}_{i}^{0} = \frac{1}{N} $, then:

For the estimated signal:

$$\begin{aligned}
{x}^{k + 1} = \arg \min_{x} \quad & \frac{1}{2} {\left\| {h}^{k} \ast x - y \right\|}_{2}^{2} + \lambda \operatorname{TV} \left( x \right) \\
\end{aligned}$$

For the estimated kernel:

$$\begin{aligned}
{h}^{k + 1} = \arg \min_{h} \quad & \frac{1}{2} {\left\| h \ast {x}^{k + 1} - y \right\|}_{2}^{2} \\
\text{subject to} \quad & \sum h = 1 \\
& {h}_{i} \geq 0 \\
\end{aligned}$$

So, in your case we can do the following:

Replace the regularization by the Sparsity Model. Solve the $ x $ iteration by the methods in Thomas' answer (Yaghoobi, Blumensath, Davies, 2007, Quantized Sparse Approximation with Iterative Thresholding for Audio Coding - DOI,  Nagahara, 2015, Discrete Signal Reconstruction by Sum of Absolute Values - DOI). Solve for $ h $ as for Least Squares with Simplex Constraint.
Use model without convolution using Dictionary and use Dictionary Learning Methods like K-SVD. For the signal estimation iteration still you should use the methods above.

Some related questions:

Using Total Variation Denoising to Clean Accelerometer Data.
The Meaning of the Terms Isotropic and Anisotropic in the Total Variation Framework.
Why Sparse Priors Like Total Variation Opts to Concentrate Derivatives at a Small Number of Pixels?
How Can I Use MATLAB to Solve a Total Variation Denoising / Deblurring Problem?
Intuitive Meaning of Regularization in Imaging Inverse Problems.
Estimation / Reconstruction of an Image from Its Missing Data 2D DFT.
Deconvolution of an Image Acquired by a Square Uniform Detector.


separating stationary and non stationary parts of univariate signal Does anyone know if there is a procedure as to separate the stationary and non stationary parts of a univariate signal. I have seen signal source separation and blind separation algorithms (all of them assuming that the signal is a linear superpossition of underlying hidden source signals where the signals are statistically independent - linearly independent - not sure if colored noises are linearly independent) where different signals are mixed in a single combined signal. For this several methods exist as to separate the sources mostly based on either orthogonality, correlation, variance and distances, however all of them based on having already multivariate variables:

Principal components analysis
Singular value decomposition
Independent component analysis
Dependent component analysis
Non-negative matrix factorization
Low-complexity coding and decoding
Stationary subspace analysis
Common spatial pattern

My question is that if you could treat a single signal which is the combination of signal + uncorrelated stationary noise colors (white) + correlated noise colors (pink,brown,etc) as a source problem (implying that there is a mix of signals formed by the stationary noise, non stationary noise and signal) and thus apply a specific blind signal separation algorithm as to separate the sources statistically (for example as in stationary subspace analysis there is identification of the stationary and non stationary parts).

All blind signal separation algos I have seen can only be applied to multiple variables, my question is if there are there any blind separation algorithms that can be applied to univariate series where stationary and non stationary components can be separated as individual signal sources.

EDIT

Following comments from LPT

I know the signal is spread across a wide range of frequencies from the Evoluionary power spectrum. I know uncorrelated stationary white noise affects the signal and that other unidentifed stationary correlated colors do as well but ignore their mix. My idea was to run a wavelet at one level decomposition as to have two range of frequencies and use this as to run one of the analysis above as to determine the stationary and non stationary parts. As the stationary parts will correspond to correlated and uncorrelated noise colors which I will assume to be non gausian and additive. As the nature of the signal is non stationary I would effectively have separated the signal from noise. Then from noise I will run wavelet thresholding as to eliminate white noise based on a MA MAD and remain with correlated noises which via FFT and their power relation understand the type dominating noise color or identify their mix.

All these algorithms are based on the determination of a new basis (i.e. a new set of basis vectors), chosen according to a given constrain (to maximize the variance of each component in PCA for instance).

If you have only one signal, you have only one basis vector. So you cannot use these algorithms because you don't have enough information.

In lower dimensions, you can think of it as the calculation of the arithmetic mean from a set of noisy values corresponding to the same measurement. You need multiple measurement to be sure to get close to the "true" value of the measurement. What you are asking is how to get the "true" value from only one measurement without any repetition. That cannot be done because you don't have enough information even to do anything.

However, depending on you signal and its specific noise, there might be other ways to achieve what you want. Actually, all you have to do is to identify a way to differentiate noise from signal.

Their distributions are obviously different, but these are mere projections, so again, you have too few data to use them.

One working possibility is to filter your data to remove high frequency components (to remove most of the noise) and very low frequency components (to remove pink noise, but usually pink noise is really not a problem). But that would work only if you have an idea of the bandwidth of your signal (the maximum frequency beyond which your signal has no meaning anymore).

What you describe won't work. If I continue my analogy of the average over noisy measurement data, it is like saying that you will split the result of you unique measurement (say, 15 = 12 + 3) and use that result (12 and 3) to calculate the mean (7.5). I hope you understand that this makes no sense (the mean of 15 is at best undefined, or at worst 15). And even if it did, your decomposition would be quite arbitrary (why not use 10+5 or 20-5 or 13264585120-13264585105?). Again, you won't get more information by decomposition of your signal.

Final word:

If you really want us to help you or to give you more details about why this won't work as you expect, please open a new question, provide actual data and explain your setup.


Cascaded one-pole sections and relationship to time constant I have implemented a typical one-pole lowpass filter to smooth out some control signals.
The form of the filter is:
$$
y1 = x0 + b1\centerdot(y1 - x0)
$$
x0 is the current input sample.
y1 is the past output sample.
b1 is a smoothing coefficient.

Now I specify the b1 coefficient in terms of how much time it will take for the filter to converge to a certain percentage of the final value.
For example, the following equation makes sure that it will take 300 ms for the signal level to converge to 1% (0.01) of the final value. SR is the sampling rate of the system.
$$
\large {b1 = e^{\frac{ln(0.01)}{300 * SR * 0.001}}}
$$
…or in a more computationally expensive form:
$$
\large {b1 = 10^{\frac{log_{10}(0.01)}{300 * SR * 0.001}}}
$$

Let's see an example:
$$
\large {b1 = e^{\frac{ln(0.01)}{5}}}
$$

This causes the one pole filter to converge to 1% of the final value in 5 samples.
Let's initialize y1 with 1 and y0 with 0.
Now let's see the first 5 samples of the impulse response of the filter:
$$
[1]\quad y1 = y0 + b1\centerdot(y1 - y0)\rightarrow y1 = 0.3981072\\
[2]\quad y1 = y0 + b1\centerdot(y1 - y0)\rightarrow y1 = 0.1584893\\
[3]\quad y1 = y0 + b1\centerdot(y1 - y0)\rightarrow y1 = 0.0630957\\
[4]\quad y1 = y0 + b1\centerdot(y1 - y0)\rightarrow y1 = 0.0251189\\
[5]\quad y1 = y0 + b1\centerdot(y1 - y0)\rightarrow y1 = 0.0100000\\
$$
This is what we wanted, the signal is now 0.01 compared to the
original 1.
The impulse response of the filter will be exactly x after n samples.

So far, so good…

Now I would like to make the transition even smoother so I started cascading more one-pole sections.
$$
y1_a = x0_a + b1\centerdot(y1_a - x0_a)\\
y1_b = y1_a + b1\centerdot(y1_b - y1_a)\\
y1_c = y1_b + b1\centerdot(y1_c - y1_b)\\
…
$$
This certainly works well and gives the wanted effect, however it also messes up my time constant.
For the same b1 coefficient the more cascaded one-pole sections the longer the filter takes to converge.
So I would like to know what value the b1 coefficient has to take in order to make the filter converge in the same amount of time for more cascaded sections.

More specifically I would like to find an equation that takes:
- time in ms
- sample rate
- number of cascaded one-pole sections

and gives:
- value of the b1 coefficient.

Is anybody so kind to take my hand and lead me through the math?
Thank you so much.

Luigi

The difference equation of your single pole filter is

$$y[n]=x[n]+b(y[n-1]-x[n])\tag{1}$$

where $n$ is the time (sample) index, $x[n]$ is the input, and $y[n]$ is the output. Its transfer function is

$$H(z)=\frac{1-b}{1-bz^{-1}}\tag{2}$$

If you cascade $N$ identical single pole filters, the total transfer function is

$$G(z)=H^N(z)=\frac{(1-b)^N}{(1-bz^{-1})^N}\tag{3}$$

The inverse $\mathcal{Z}$-transform of $(3)$ is the impulse response of the total system:

$$g[n]=(1-b)^N\frac{(n+1)(n+2)\ldots (n+N-1)}{(N-1)!}b^n,\quad n\ge 0,\quad N>1\tag{4}$$

From $(4)$ it's straightforward to compute the desired value of $b$ for $g[n]$ to decay to $0.01$ of its initial value in a given time $T$. The number of sample values corresponding to time $T$ is $K=\text{round}(Tf_s)$, where $f_s$ is the sampling frequency. We now need to solve

$$g[K]=0.01\cdot g[0]$$

for $b$, which results in

$$b=\left(\frac{0.01\cdot (N-1)!}{(K+1)(K+2)\ldots (K+N-1)}\right)^{1/K}\tag{5}$$

For $K=5$, as in your example, and for $N=2$ cascaded filters you get from $(5)$ a value of $b=0.27821$. For $N=5$ you would need to decrease this value to $b=0.15133$.


What is a good way to distinguish aircraft noise from other sounds? The FAA's NextGen routes, which were modeled, not measured, have created a bunch of aircraft noise for tens of thousands of people. I want to make it easy for those people (including me) to buy a cheap USB microphone, plug it into a laptop or PC, detect aircraft noise, correlate that with publicly available flight-path data, and file complaints with the FAA, airport(s), and legislators. With enough actual (as opposed to theoretical) data, we might be able to get them to roll back the route changes. Or not, but it's worth a try.

To do that, I need some help with the "detect aircraft noise" part of the problem. It's not enough to trigger off noise above a certain dB threshold, because lawnmowers, leaf blowers, trucks, etc. also make noise. I can weed out some over-threshold noise that doesn't last long enough or lasts too long (aircraft take 5-10 seconds to pass overhead, vs. minutes for lawnmowers and a couple of seconds for trucks), but that's a pretty crude approach. And I still need to be able to factor multiple noise sources into individual components.

Correlation would probably be easier with a pair or array of microphones, to establish noise direction that could be tied to flight paths, but that adds a lot more hardware and complexity. Ideally we'd just need one microphone.

Aircraft include airplanes (jets, turboprops, and standard propellor-driven), as well as helicopters.

I'm a good coder (Java expert, former C expert), but have zero experience with DSP or training algorithms. Can you help me get started?

You need to train a classifier. The way I would approach it is to allow users to submit audio recordings through your app, and indicate which sections of the recordings contain the plane. This would enable you to quickly crowd source a lot of data with which to build a model. Then you would find the STFT and volume of the signal over time. Finally, train the classifier (such as an SVM or neural network) using the amplitude spectrum and volume features. Once you have trained the classifier, you can use it to detect planes by extracting the same features on live audio (STFT isn't computationally expensive, so it should be feasible to do on the phone).

Here's a backgrounder on Audio signal classification.


How to calculate the total signal strength and phase shift of multiple signals in a reciever I have a simulation to make where I have an array of transmitters that transmit the same signal. At a random point, which I have to consider as a receiver I have to measure the phase shift of the signals and also measure the total signal strength. The requirement is that the transmitter emits the same signal but with different phase offsets.

As far I understood because the transmitters are separated apart, the received signal is already offset by some phase because they different travel distances. Is my assumption right?
I know the formula to find the phase offset given the distance travelled of two signals - how do I calculate it for multiple signals? Or am I understanding the question wrong?
How do I calculate the total signal strength for all signals combined? I know for one signal assuming a free space, line of sight connection I can use friis equation. But how do I use it for multiple signals?

Any help is very much appreciated. I am Software Engineer who have recently taken up Wireless communications course, and I'm not sure if I understand the motivation and meaning behind the simulation right. Thanks!

This is simplified assuming there is not multipath occurring in the transmission (where the receiver would typically receive multiple copies of the same transmit signal at different delays resulting in multipath distortion which are typically addressed using channel estimation and equalization).

The OP's phase shift would all be relative since there is no common reference (clock) with the transmitter, so the comparison would be of the phase of the received signals to each other. A cross-correlation of the received signals in complex notation with the reference signal of what you know was transmitted would provide this.

This is accomplished by doing a multiply and accumulate of the received signal with the complex conjugate of the known ideal transmit signal, repeating the computation below for each possible offset $m$  (thus is the cross-correlation function specifically, correlating at each possible offset in time).

$$XCorr[m] = \sum r[n-m]t^*[n]$$

Where (*) indicates a complex conjugate.

The magnitude of the correlation peak would be proportional to the signal strength, and the shift $m$ is proportional to the delay. With two signals at the same shift $m$ (in the same bin), the phase of the correlation is proportional to the carrier phase, giving a fine delay indication by comparing the phase result from two different receiver signals if they were in the same bin. In this case the signals would need to not be correlated themselves, otherwise the result would be the vector addition of the two (appearing as a single received signal). For this reason the transmit signals from each transmitter would be uncorrelated by separating them in frequency, time or code; and you can identify each one separately in receiver using correlation.

Note the above computation for the cross correlation function can be done directly with FFTs given the FFT correlation property (which results in a circular cross-correlation):

$$XCorr[m] = \text{ifft}\{\text{fft}(r[n])\}\text{fft}^*\{(t[n])\}$$


DTFT fourier transform (modified property) I know there are 3 properties of DTFT that help with my problem

$$ a^{n}u[n]=\frac{1}{1-ae^{-jΩ}} $$
$$ (n+1)a^{n}u[n]=\left(\frac{1}{1-ae^{-jΩ}}\right)^{2} $$
$$ \frac{(n+r-1)!}{n!(r-1)!}a^{n}u[n]=\left(\frac{1}{1-ae^{-jΩ}}\right)^{r} $$

But I cannot find some use between them to calculate the DTFT of the following signal
$$ x[n]=(n+5)\left(\frac{7}{45}\right)^{n}u[n] $$ where $u[n]$ is the unit step function. Can anyone help?

Hint:

Note that you can express your signal as

$$x[n]=(n+1+4)\left(\frac{7}{45}\right)^{n}u[n] = (n+1)\left(\frac{7}{45}\right)^{n}u[n] + 4\left(\frac{7}{45}\right)^{n}u[n] $$

So now you have the sum of two signals:

$$x_1[n]=(n+1)\left(\frac{7}{45}\right)^{n}u[n]$$
$$x_2[n]=4\left(\frac{7}{45}\right)^{n}u[n]$$

Using your properties and the fact that the Fourier transform is linear, you should be able to solve it.


Help implementing audio dynamic range compression I am trying to implement audio dynamic range compression in JavaScript (not using web audio API).

There is a lot of articles for sound technicians, and some high level documentation, but I couldn't find any helpful reference for actually implementing digital dynamic range compression.

From what I understand, there is at least 3 steps in calculating the rectified signal.

computing the input level
computing the gain to apply to the signal
applying the gain

I process the audio in blocks, so for 1) I was thinking of computing the RMS of one block

Any pointers to a good reference? Or anybody willing to explain me a bit the steps required to implement this?

Here are some suggestions:

There are plenty of opensource implementations (Sox, Audacity, etc). Even if you don't understand them, you might be able to translate the code from C to javascript.
I'm not aware of a good explanation of the process online, but there are plenty of books on the subject:

Digital audio signal processing covers this topic and is well written. (As does DAFX, but DAFX is poorly organized and the coverage is less straightforward)
Digital Audio with Java also covers this topic and comes with working Java code that should be easy to translate to other languages, such as javascript. This book has many many flaws, but it is good for someone with no audio programming experience.

The principle is to create an envelope of the signal (controlled by attack and release), shape that envelope using using some transfer function (controlled by ratio, threshold and knee), and then apply that result back to the original signal. A makeup gain stage often follows.

@Deve's answer suggests some possible transfer functions.


How to determine the poles from a graph From my knowledge of stability, I understand that if the function approaches a finite number then the system will be stable. Thus if a system is stable its poles will be on the left of the $j\omega$-axis. This means that $Re(p_{1}, p_{2}) < 0 $. What I don't understand is how you determine whether the Imaginary part is equal to 0 or not.

If you have complex conjugate poles (i.e., $\text{Im}\{p_i\}\neq 0$) you have oscillations in the step response. Only real-valued poles will give you a monotonic step response as shown in the figure.

Take a look at this related answer.


Image Matching on iOS Was recommended to migrate this question from SO.

I am building an iOS app that, as a key feature, incorporates image matching. The problem is the images I need to recognize are small orienteering 10x10 plaques with simple large text on them. They can be quite reflective and will be outside(so the light conditions will be variable). Sample image http://i1274.photobucket.com/albums/y425/Chris_Mitchelmore/2_zpsce84d4f3.png

There will be up to 15 of these types of image in the pool and really all I need to detect is the text, in order to log where the user has been.

The problem I am facing is that with the image matching software I have tried, aurasma and slightly more successfully arlabs, they can't distinguish between them as they are primarily built to work with detailed images.

I need to accurately detect which plaque is being scanned and have considered using gps to refine the selection but the only reliable way I have found is to get the user to manually enter the text. One of the key attractions we have based the product around is being able to detect these images that are already in place and not have to set up any additional material.

Can anyone suggest a piece of software that would work(as is iOS friendly) or a method of detection that would be effective and interactive/pleasing for the user.

Sample environment:
http://www.orienteeringcoach.com/wp-content/uploads/2012/08/startfinishscp.jpeg

The environment can change substantially, basically anywhere a plaque could be positioned they are; fences, walls, and posts in either wooded or open areas, but overwhelmingly outdoors.

I managed to find a solution that is working quite well. Im not fully optimized yet but I think its just tweaking filters, as ill explain later on.

Initially I tried to set up opencv but it was very time consuming and a steep learning curve but it did give me an idea. The key to my problem is really detecting the characters within the image and ignoring the background, which was basically just noise. OCR was designed exactly for this purpose.

I found the free library tesseract (https://github.com/ldiqual/tesseract-ios-lib) easy to use and with plenty of customizability. At first the results were very random but applying sharpening and monochromatic filter and a color invert worked well to clean up the text. Next a marked out a target area on the ui and used that to cut out the rectangle of image to process. The speed of processing is slow on large images and this cut it dramatically. The OCR filter allowed me to restrict allowable characters and as the plaques follow a standard configuration this narrowed down the accuracy.

So far its been successful with the grey background plaques but I havent found the correct filter for the red and white editions. My goal will be to add color detection and remove the need to feed in the data type.


Value of $A_k$ in Fourier series Fourier series in continuous time domain while representing $a_k$ in rectangular form   $$ a_k = B_k + jC_k$$
But when using the value of $a_k$ in the main equation: $$ x(t) = a_0 + 2\sum^{+\infty}_{k\ =\ 1} [B_k\cos k\omega_0t - jC_k\sin k\omega_0t]$$

I want to ask from where $-$(minus) sign comes as there was no minus sign
main Equation
$$ x(t) = a_0 + 2\sum^{+\infty}_{k\ =\ 1} 2\mathrm{Re}\left[ e^j(k \omega_0t +\theta_k)\right]$$

I'm afraid that all the formulas in your question are inaccurate. From what I understand, you're interested in the representation of real-valued signals $x(t)$:

$$x(t)=a_0+2\Re\left\{\sum_{k=1}^{\infty}a_ke^{jk\omega_0t}\right\}=
a_0+2\sum_{k=1}^{\infty}\Re\left\{a_ke^{jk\omega_0t}\right\}\tag{1}$$

With $a_k=B_k+jC_k$ and $e^{jk\omega_0t}=\cos(k\omega_0t)+j\sin(k\omega_0t)$ you get

$$\Re\left\{a_ke^{jk\omega_0t}\right\}=B_k\cos(k\omega_0t)-C_k\sin(k\omega_0t)\tag{2}$$

which you can substitute into Eq. $(1)$ to get the equation in your question (but without the $j$!). The minus sign in $(2)$ appears simply because $j\cdot j=-1$.


MATLAB: Prepare an RGB image for digit and letter recognition I need a Matlab function that prepares an image for digit and letter recognition.
What I need now is to convert the original RGB image to a binary image that every pixel in it is white except the pixels corresponding to letters and digits, as well as all digit and letter must appear chromatic/saturated, i.e. appear full of color.

Here is the code I have been tested. As you can see some pixels of a letter or digit are white.
I = imread('img6.png');   % read the image into the matrix
Ig = rgb2gray(I);

Icon = imadjust(Ig);

subplot(2,2,1)
imshow(Ig)
subplot(2,2,2)
imshow(Icon)
subplot(2,2,3)
imhist(I)
subplot(2,2,4)
imhist(Icon)

1- How we can convert the original image to a high contrast image?
2- How the shadows around letters and digits can be removed?

First, many OCR programs automatically perform image processing, including possible color conversion, image filtration, dilation/erosion, and monochrome conversion, in order to improve the quality of OCR operation. Especially good in this are neural-network enabled cloud-based OCR applications.
If you absolutely need local MATLAB program to perform the image processing supposedly preparing the document for recognition, first check if MATLAB or its addons have some kind of image thresholding function: the function which outputs a monochrome (black/white) image given a scanned document. If there is no ready thresholding function, you can try and code it. An introduction into a computer graphics required to successfully code an image preparation operations "for digit and letter recognition" you will find in https://tesseract-ocr.github.io/tessdoc/ImproveQuality and elsewhere, if you search words and phrases highlighted in this answer.
You may find useful this ref, which, although lacking MATLAB code, contains code snippets written in common computer languages.


Design of discrete time second order all pass filter I want to design a second order  digital all pass filter with the transfer function given by
$$H(z)=\frac{a_2 +a_1 z^{-1}+z^{-2}}{1+a_1 z^{-1}+a_2 z^{-2}}$$
The input to this filter has frequencies ranging from $1$ to $8\ \rm kHz$. The sampling frequency is $20\ \rm kHz$.
How to design the filter coefficients $a_1, a_2$ such that for a particular input frequency, I have to apply 180 degree phase shift? The frequency which is given a 180 degree phase shift will vary. So, how to vary the filter coefficients based on which frequency has to be phase shifted by 180 degree?

HINTS:

What value does $H(z_0)$ have if there's a phase shift of $180$ degrees (for $|z_0|=1$)?
Equate $H(z_0)$ to that value and write down the resulting equation for $a_1$ and $a_2$ and $z_0$.
You want a stable all-pass filter, so choose $a_2$ such that a given desired pole radius $r$, $0<r<1$, is achieved if you assume two complex conjugate poles.
With that value of $a_2$, and with a given frequency $\omega_0$ for which the phase shift is $180$ degrees, solve the equation you got in 2. for $a_1$. The result will have the form $$a_1=f(a_2)\cos(\omega_0)$$ where $f(\cdot)$ is a simple affine function. Note that $\omega_0$ is normalized by the sampling frequency.


Fir filter not working according to frequency response I have designed a FIR filter with following parameters

fs=30,
f_pass1 = 1, stop1 = 0.5, f_pass2 = 3, f_stop2 = 4,
A_pass = 1,A_stop = 40

Here is code :

fs_filter=30;

f_pass1 = 1; %1.2
f_stop1 = 0.5;%0.7

f_pass2 = 3;
f_stop2 = 4; %6

A_pass = 1;
A_stop = 40;

del_f1 = -f_stop1 + f_pass1;
del_f2 = f_stop2-f_pass2;

del_f = min(del_f1,del_f2);

fc1 = f_pass1 -(( del_f)/2);
fc2 = f_pass2 +( (del_f)/2);

w_c1 = (2*pi*fc1)/(fs_filter);
w_c2 = (2*pi*fc2)/(fs_filter);

S_pass = ((10^(A_pass/20))-1)/((10^(A_pass/20))+1);
S_stop = 10^(-A_stop/20);

S = min(S_pass ,S_stop);

A = -20*log10(S)

alpha =0;
D = 0.922;

if A > 21 && A < 50
    alpha = 0.5842*((A-21)^0.4)+0.07886*(A-21);
    D = (A-7.95)/14.35;
elseif A>=50
    alpha = 0.1102*(A-8.7);
    D = (A-7.95)/14.35;
end

N = floor(((D*fs_filter)/del_f)+1);

N = N + mod(N-1,2)

M = (N-1)/2;

coeff= zeros(1,N);
w   = zeros (1,N);
d   = zeros (1,N);

for i=0:N-1
    x_b= sqrt(1-(((i-M)/M)^2));

    w(i+1) = besseli(0,(alpha*x_b))/besseli(0,alpha);

    if i ~= M
    d(i+1) =((w_c2/pi)*(sin(w_c2*(i-M)))/(w_c2*(i-M)))-((w_c1/pi)*(sin(w_c1*(i-M)))/(w_c1*(i-M)));
    coeff(i+1) = w(i+1)*d(i+1);
    else
    d(i+1) = ((w_c2-w_c1)/pi);
    coeff(i+1) = w(i+1)*d(i+1);
    end

end

fft_coeff = abs(fft(coeff));
num_bins = length(fft_coeff);
bins = 0 : num_bins-1;
bin_Hz = bins*fs_filter/num_bins;
N_2 = ceil(num_bins/2);

plot(bin_Hz(1:N_2) , fft_coeff(1:N_2),'x');

but when i am convoluting it with the sin wave of say frequency f ( 0.1 to 5) sampled at fs = 50 the output of the signal is not actually following the response.

Code for generating signal :-

    fs=50;
    f0 = 0.6593;
    cyc = 50;
    t = 0:1/fs:cyc/f0;
    x = sin(2*pi*f0*t);
    y = conv(x,coeff,'valid');

am i missing something ?

Your normalized stop band and pass band edges are:

$$\begin{align}f_{s1}&=0.5/30=0.017\\
f_{p1}&=1/30=0.033\\
f_{p2}&=3/30=0.1\\
f_{s2}&=4/30=0.133\end{align}$$

And the normalized frequency of your sinusoidal input signal is $0.6593/50\approx 0.013$, which is smaller than $f_{s1}$, so it's in the filter's stop band.

Note that you've designed your filter for a certain sampling frequency, and that it is only the normalized band edges that are fixed; the absolute frequencies of the band edges change according to the sampling rate.


Applying 2D Image Convolution in Frequency Domain with Replicate Border Conditions in MATLAB I have created a function that filters an image (250x250) with a gaussian blur kernel (5x5) using FFT and IFFT. I am trying to get my filtered image to equal exactly the filtered image created by the 'imfilter' function using the 'replicate' option. There are small differences in the two filtered images. They look identical visually but upon inspecting the pixel values, they are slightly different.
This is the general logic my code follows:
Zero pad the borders of the input image to size 254x254
Zero pad the kernel to size 254x254
'clamp to edge' the edge pixels of the input image. The outermost input image matrix rows and columns are extended to the borders, replacing the zeros. I am doing this to emulate the 'replicate' option in 'imfilter'.
Use 'fft2' to convert both the input image and kernel to frequency domain
Element wise multiply the two together
Convert the result to the time domain using 'ifft2'.
Crop away the border pixels, converting the filtered image from 254x254 ---> 250x250
What am I doing wrong here? Thanks.

The trick is to properly compensate for the fact that Frequency Domain multiplication applies a convolution with the circular boundary conditions in the spatial domain.
You may use the following code:
clear('all');
close('all');

gaussianKernelStd       = 2;
gaussianKernelRadius    = ceil(5 * gaussianKernelStd);

mI = im2double(imread('cameraman.tif'));
mI = mI(:, :, 1);

numRows = size(mI, 1);
numCols = size(mI, 2);

vX = [-gaussianKernelRadius:gaussianKernelRadius].';
vK = exp(-(vX .* vX) ./ (2 * gaussianKernelStd * gaussianKernelStd));
mK = vK * vK.';
mK = mK ./ sum(mK(:)); %<! The Gaussian Kernel

mIPad = padarray(mI, [gaussianKernelRadius, gaussianKernelRadius], "replicate", 'both'); %<! For the replicate array

mKC = CircularExtension2D(mK, size(mIPad, 1), size(mIPad, 2)); %<! Circular extension for the 2D Kernel

startIdx = gaussianKernelRadius + 1;

mIFiltered      = ifft2(fft2(mIPad) .* fft2(mKC), 'symmetric');
mIFiltered      = mIFiltered(startIdx:(startIdx + numRows - 1), startIdx:(startIdx + numCols - 1)); %<! Removing the padding
mIFilteredRef   = imfilter(mI, mK, 'replicate', 'same', 'conv'); %<! Reference

figure(); imshow(mIFiltered);
figure(); imshow(mIFilteredRef);

max(abs(mIFilteredRef(:) - mIFiltered(:))) %<! Should be very very low

The function CircularExtension2D() is given In my StackExchange Signal Processing Q38542 GitHub Repository.  It was taken from my answer to Applying Image Filtering (Circular Convolution) in Frequency Domain.
The steps the code implements are as following:

Pad the image in order to have Replicate boundary condition convolution.
Convert the spatial domain kernel into a form which matches the image in frequency domain. We assume top left of the image is (0, 0) in spatial domain. So we need the (0, 0) of the kernel to also be in the top left corner.
Apply circular convolution using frequency domain.

As you can see, the result is prefect.
In my answer to How Much Zero Padding Do We Need to Perform Filtering in the Fourier Domain? I implemented a MATLAB Function, ImageFilteringFrequencyDomain(), to apply Frequency Domain Convolution with the border conditions supported in imfilter().
Remark
In your case, the kernel is 5x5 which is very small.
For small kernels and a single image it is better to apply the convolution in the spatial domain. It will be much faster.


How to reconcile "bandwidth" in the Shannon-Hartley Theorem with a spread-spectrum principle I'm trying to understand something about channel "bandwidth" $B$ in the Shannon-Hartley Theorem: $$C = B \log_2 \left( 1 + \dfrac{S}{N} \right) $$

Suppose I have a bitstream encoded as a signal that occupies some bandwidth $B_1$ which is less than channel bandwidth $B$, and I have additive white Gaussian noise $N$. Shannon-Hartley tells me that I'm not using the communication channel efficiently; I should be bumping up the bandwidth of my encoder.

One way of doing that is to use spread-spectrum... but I've read several places (including this paper) that SS doesn't help against wideband noise; it only allows you to increase SNR if you have narrowband noise.

Here are my questions:

If I'm only using bandwidth $B_1$, doesn't that mean I can filter some of the noise out? And if so, would I be able to get $S/N$ down up? Wouldn't that be an alternative to occupying the full channel bandwidth?
Is spread-spectrum a bad idea to try to occupy the full channel bandwidth?

If I'm only using bandwidth B1, doesn't that mean I can filter some of the noise out? And if so, would I be able to get S/N down? Wouldn't that be an alternative to occupying the full channel bandwidth?

Capacity increases linearly with bandwidth, but only logarithmically with SNR. So, increasing SNR by decreasing bandwidth is a bad idea if your noise isn't narrowband.

  Is spread-spectrum a bad idea to try to occupy the full channel bandwidth?

No, but there's more at play here:

Remember, $C$ is just an upper limit assuming you somehow (and it's not generally known how to do that) build the perfect transceiver including perfect channel coding.

Spread-Spectrum techniques can make for very robust transceivers, but they'd typically require $B_1 \ll B$. However, if $B$ is large, you'll likely be met with a non-flat channel, and that means things get a lot harder to receive correctly, and you'll end up building a complicated equalizer to counter ISI.

Especially fast-changing environments, which might be that the receiver, the transmitter, or potential paths change position or phase, getting a channel state estimate sufficient to do the equalization becomes hard enough to require you to add more pilots or otherwise redundant data to your transmission that the increase in bandwidth might not pay.
That's the main reason why for high-bandwidth systems (Wifi, cellular data), the world is leaving spread-spectrum behind in favour of multicarrier systems (typically, OFDM these days, but also FBMC): this has happened to WiFi (802.11b was DSSS, a/g/n/ac/p... use OFDM) and to Cellular standards (3G / UMTS was typically Spread Spectrum, 4G/LTE/4G+/5G is OFDM (and might turn out to do FBMC in the future for specific cases)).

Still, for smaller bandwidth systems, where the channel can be assumed to be tolerably flat, DSSS and alike play a very important role, because they allow you to send at a fixed PSD (say, $x$ dBm / MHz) over a larger bandwidth and thus allow for vastly improved SNR due to processing gain; in effect, you transmit the same bits, but you can put more power into them legally, because you're using a larger bandwidth. Yes, you infer the same increase in noise power; but noise doesn't correlate with your spreading pattern, so you'll have a net gain.

  it only allows you to increase SNR if you have narrowband noise.

Still haven't read that paper, so I don't know if you're misinterpreting the paper or if it's wrong (might simply make different assumptions):

Let's look at simple DSSS.

For every finite-power signal, we interpret its variance as power.

Let's assume we're dealing with a spread sequence $c_i$, $i=0,\ldots, L-1$, with $\left\lvert c_i\right\rvert = 1$. Let the symbols sent be $s$, $|s|=1$. This all works pretty well w.l.o.g., if you consider that the below also works with expected values instead of fixed amplitudes.

Then, what the transmitter sends for a single data symbol is a sequence $(sc_0, sc_1, \ldots, sc_{L-1})$. And at the correct sampling instant, the receiver calculates, from the receive samples $r_i$:

$$R = \sum_{i=0}^{L-1} c_i r_i \text.$$

In the noise-free case:

$$\begin{align}
R_\text{no noise} &= \sum_{i=0}^{L-1} c_i sc_i \\
&= s \sum_{i=0}^{L-1} c_i c_i \\
&= s L c_i^2 \\
&= s L 1^2 \\
&= s L
\text,\end{align}$$

which has the variance of $L^2$, since we scale $s$, which has a variance of $1$, with a constant $L$.

In the noise-only case:

$$\begin{align}
\text{Var}\left(R_\text{only noise}\right) &= \text{Var}\left(\sum_{i=0}^{L-1} c_i n_i \right)&\text{i.i.d}\\
&=L\text{Var}\left(n\right)
\text.\end{align}$$

So, the sum of noise has $L$ times the variance of the individual noise sample. Let's look at the SNR:

$$\text{SNR} = \frac{\text{Var}\left( R_\text{no noise}\right)}{\text{Var}\left( R_\text{only noise}\right)} = \frac{L^2}{L\cdot N} = L \text{SNR}_\text{unspread}\text.$$

So, I'd say, yes, spreading gain works clearly in your favor here.


Simulating BPSK Costas Loop in Matlab I am trying to implement BPSK demodulator in FPGA. First, I am modeling the system in MATLAB. Now, I am trying to design Costas Loop for carrier synchronization. This is the block diagram I'm trying to achieve:

I'm using RRC as a low pass filter. This is the MATLAB code:

load RRC.mat

fc = 0.05; % 0.11  %0.09
N = length(rf_signal);
bb = zeros(1, N);
bb_f = zeros(1, N);
I_f = zeros(1, N);
Q_f = zeros(1, N);
I_r = zeros(1, N);
Q_r = zeros(1, N);
error = zeros(1, N);
PhErr = zeros(1, N);
error_integral = zeros(1, N);
Theta = zeros(1, N);

% Loop Filter Coefficients
BW = 100;  % Hz
loop_theta = 2*pi*BW;
C1 = 4*(loop_theta)^2/(1+sqrt(2)*loop_theta+loop_theta^2);
C2 = 2*sqrt(2)*loop_theta/(1+sqrt(2)*loop_theta+loop_theta^2);

for i = 2:N
    % Downconverting to Baseband
    bb(i) = rf_signal(i).*exp(j*2*pi*fc*i).*exp(j*Theta(i-1));

    % Filtering
    bb_f = filter(RRC, bb(1:i));
    I_f(i) = real(bb_f(i));
    Q_f(i) = imag(bb_f(i));

    % Error
    error(i) = I_f(i).*Q_f(i);

    % Loop Filter
    error_integral(i) = error(i).*C1 + error_integral(i-1);
    PhErr(i) = error(i).*C2 + error_integral(i);

    % Phase Accumulator
    Theta(i) = Theta(i-1) + PhErr(i);

end

figure; subplot 221; plot(real(bb)); title('I Channel')
subplot 222; plot(imag(bb)); title('Q Channel')
subplot 223; plot(I_f); title('I Channel Filter')
subplot 224; plot(Q_f); title('Q Channel Filter')

figure; subplot 311; plot(error); title('Phase Error');
subplot 312; plot(PhErr); title('Loop Filter');
subplot 313; plot(Theta); title('Control Signal \theta');

rf_signal comes from another script that generates BPSK signal.
However, when I close the loop and set frequency exactly as the transmitter's I get these waveforms:

Although the signal demodulates very well when I open the loop. I have tested the loop with a cosine wave and it works just fine by producing a DC signal.

Any advice would be appreciated.

I figured it out. It turned out that I don't really need that 'Phase Accumulator', and I should reverse the control signal in the NCO. So, I use -PhErr(i-1) as a control signal to NCO. Here is the modified code:

% Siraj Muhammad
% 25/3/2015
% BPSK Demodulator
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load RRC.mat

fc = 0.05000001;
phase_offset = pi/7;

N = length(rf_signal);
bb = zeros(1, N);
bb_f = zeros(1, N);
I_f = zeros(1, N);
Q_f = zeros(1, N);
I_r = zeros(1, N);
Q_r = zeros(1, N);
error = zeros(1, N);
PhErr = zeros(1, N);
error_integral = zeros(1, N);

% Loop Filter Coefficients
BW = 100;  % Hz
loop_theta = 2*pi*BW;
C1 = 4*(loop_theta)^2/(1+sqrt(2)*loop_theta+loop_theta^2);
C2 = 2*sqrt(2)*loop_theta/(1+sqrt(2)*loop_theta+loop_theta^2);

for i = 2:N
    % Downconverting to Baseband
    bb(i) = rf_signal(i).*exp(j*2*pi*fc*i+j*phase_offset).*exp(-j*PhErr(i-1));

    % Filtering
    bb_f = filter(RRC, bb(1:i));
    I_f(i) = real(bb_f(i));
    Q_f(i) = imag(bb_f(i));

    % Error
    error(i) = I_f(i).*Q_f(i);

    % Loop Filter
    error_integral(i) = error(i).*C1 + error_integral(i-1);
    PhErr(i) = error(i).*C2 + error_integral(i);
end

figure; subplot 321; plot(real(bb)); title('I Channel')
subplot 322; plot(imag(bb)); title('Q Channel')
subplot 323; plot(I_f); title('I Channel Filter')
subplot 324; plot(Q_f); title('Q Channel Filter')

subplot 325; plot(error); title('Phase Error');
subplot 326; plot(PhErr); title('Loop Filter');

I added a variable phase_offset to test the loop with different offsets, and I can drift the frequency a little bit and the loop will lock very nicely.

I hope this helps somebody as I couldn't find any MATLAB code for Costas Loop online.


Expected Value of a sequence with two random variables If I have a signal of the form $x\left(n\right)=Acos\left(nω+ϕ\right)$
where
$\omega \in \left[\omega -\lambda ,\omega \:+\lambda \:\right]$ is a uniform random variable
and
$\phi $ is also a uniform random variable $\left[-\pi ,\pi \right]$.
I need to find the mean and the autocorrelation.
So we know the mean $E\left[x\left(n\right)\right]\:=\:AE\left[cos\left(n\omega +\phi \right)\right]$.
I think the $E\left[x\left(n\right)\right]\:=0$ regardless of the random variables just because of the periodicity of the cosine sequence.
Now I would also like to find the auto-correlation function
$R\left(n,k\right)=\:E\left[x\left(n\right)x\left(n+k\right)\right]=A^2E\left[cos\left(n\omega +\phi \right)cos\left(\left(n+k\right)\omega \:+\phi \:\right)\right]=\frac{A}{2}^2\left\{E\left[cos\left(2\left(n\omega +\phi \right)+k\omega \right)\right]+E\left[cos\left(k\omega \right)\right]\right\}=\frac{A}{2}^2cos\left(k\omega \right)$
I am not sure if am solving this the right way, can somebody please tell me if am doing something wrong? because I think the random variables must have a role in this

The "periodicity of the cosine" has little to do with the matter at hand. Indeed, unless $\omega/\pi$ is a rational number, it is not even the case that the discrete-time sequence $x[n] = A\cos(n\omega + \phi)$ is a periodic sequence at all. See, for example, this question which asks "Why is $\cos\left(\frac n6\right)$ aperiodic?" Now, coming to the meat of the question and correcting meaning-less notation such as $\omega \in [\omega-\lambda,\omega+\lambda]$ to $\omega \in [\omega_0-\lambda,\omega_0+\lambda]$, we have from LOTUS, the law of the unconscious statistician, that
$$E\left[x[n]\right] = E[A\cos(n\omega + \phi)] = A
\int_{-\pi}^{\pi}\int_{\omega_0-\lambda}^{\omega_0+\lambda}\cos(n\omega + \phi)f_{\omega,\phi}(\omega,\phi)\, \mathrm d\omega \,\mathrm d\phi \tag{1}$$ where $f_{\omega,\phi}(\omega,\phi)$ is the joint density of $\omega $ and $\phi$. The OP has not stated that $\omega$ and $\phi$ are independent random variables, but if they are assumed to be independent (as is most likely the intention), then the joint density in $(1)$ factors into the product of the marginal densities (which are given to be uniform on $[\omega_0-\lambda,\omega_0+\lambda]$ and $[-\pi,\pi]$ respectively,) and so $(1)$ can be written as $$E\left[x[n]\right] =  A
\int_{-\pi}^{\pi}\frac{1}{2\pi}\left[\frac{1}{2\lambda}\int_{\omega_0-\lambda}^{\omega_0+\lambda} \cos(n\omega + \phi) \, \mathrm d\omega\right] \,\mathrm d\phi. \tag{2}$$
Similarly, the autocorrelation function is given by
$$E\left[x[n]x[m]\right] =  A
\int_{-\pi}^{\pi}\frac{1}{2\pi}\left[\frac{1}{2\lambda}\int_{\omega_0-\lambda}^{\omega_0+\lambda} \cos(n\omega + \phi)\cos(m\omega + \phi) \, \mathrm d\omega\right] \,\mathrm d\phi. \tag{2}$$


Feature vector for music recommender system For building a music recommender system(content based) which signal features should i use as feature factor ?? Will it be peridogram,spectrogram or fft or anything else ? This system will recommend music corresponding how much they are similiar. So which features correspond strongly to how much signals are similiar ? Assume correlation measures signal similiarity (but i havenot seen any implementation based on correlation on signals).

For a big data, deep learning approach, you could try constructing an autoencoder from something such as WaveNet to achieve dimensionality reduction to get a feature vector. Alternatively, you could take STFT (short time Fourier transform) spectra of the songs over time and feed them into an RNN, again in autoencoder fashion. My intuition though is that, due to the immense complexity of musical genres, and current limitations of consumer hardware, any approach based on the raw audio alone is going to be very disappointing (unless you have a GPU farm).

Really grossly, you may get some results by finding the BPM/tempo, the spectral density, overall entropy, variance in entropy over time.

Actually, look at some of the feature extractors used in this Kaggle EEG contest. This may give you some rudimentary results, especially combined with XGB (extreme gradient boosting)

I believe most music recommendation services work like search engines, using an algorithm similar to PageRank to create a graph which associates artists.


What is the Autocorrelation of the Impulse Response if the Magnitude of Frequency Response is >1 So by the DTFT pairs, if the magnitude of the frequency response of a signal is 1, then the autocorrelation is the Kroneker Delta Function.

What if I find that the magnitude of the frequency response of the system is greater than one?

Is the autocorrelation just the constant multiplied by the Kroneker Delta Function?

Make sure not to confuse the power spectral density (PSD) with the frequency response.

PSD is the Fourier transform of the autocorrelation function.
Frequency response is the Fourier transform of impulse response.

That aside, your question is about the Fourier transform of a constant $c$.
Since $$\mathcal{F}\{1\}=\delta(\omega),$$
considering linearity of the Fourier transform we have $$\mathcal{F}\{(c)1)\}=(c)\delta(\omega)$$


How to extract features from nonstationary signal using SVD (singular value decomposition)? I am trying to extract features from a nonstationary signal (frequencies and time locations). For this purpose, I computed a wavelet transform (complex morlet) of the signal in MATLAB and obtained the scalogram as follows:

Than I tried to decompose the scalogram matrix by SVD, in order to localize the peaks and the spreads of energy in the image. The problem is that every matrix obtained from the decomposition (one for every singular value) contains informations about more than one single atom (peak in the image). I've read many articles about rotating the SVD basis in order to decompose the original matrix in matrices connected each to a single atom, but I'm not able to apply it in practice.
If we consider the standard SVD:

A=U*S*transpose(V)

a source in particular:

  [Bashor & Kareem, Efficacy of Time-Frequency Domain System
  Identification Scheme Using Transformed Singular Value Decomposition]

talks about a new SVD given by:

A=Y*Z*transpose(X)

where Y, Z, and X are the transformed singular vector and singular value matrices, defined as:

Y=U*C;
X=V*D;
Z=transpose(C)*S*D;

and the trasformation matrices are found by:

E[Y]=transpose(C)*M*C

"where C are the eigenvectors of M that maximize the mean and M is a matrix of the singular vectors". My problem is to understand the mean of the latter sentence and to apply it in MATLAB. Please help me!

What you're describing is a spectrum estimation based on signal representation in a vector space spanned by eigenvectors of a particular matrix generated from your signal observation.

That reminds me extremely much of MUSIC, and ESPRIT. I think the best way to go here is to read up on spectral estimation using MUSIC, which is based on the idea that there's an autocovariance matrix (which can be easily estimated by the dyadic product of a receive signal vector with itself, for example), and that autocovariance matrix can, using an EVD (which is a special case of the SVD for quadratic matrices), be diagonalized like you do it, into a an upper, lower and a diagonal matrix. Using the magnitudes of the eigenvalues (which happen to be the entries of the diagonal matrix), you can select the eigentvectors to span a sub-space that only contains autocovariance Eigenvectors (i.e. periodicities!) that are actually in the signal; the energy in the complementary sub-space is pretty much only noise.

Using that technique, you can directly get a function that you can sample to get a pseudospectrum (MUSIC), or through elegant rotations, the set of strongest frequencies in the signal (ESPRIT).

So, I don't get the full context of the "Eigenvectors that maximize the mean", but maybe you don't have to – you can walk on the shoulders of other SVD-based spectral estimation papers, and maybe that'll help you understand the one you're currently reading (or you just build a scalogram matrix and hit it with MUSIC as if it was an autocovariance matrix estimate, and see what happens – I love just experimenting with algorithms :) ).


Proving that a product of matrices invertible Given $R_x$ a Positive Definite (PD) covariance matrix of size $M\times M$ and $C$ a full rank $M \times N$ matrix, I want to prove that $C^* R_x^{-1} C$ is invertible to derive the Linearly Constrained Minimum Variance Beamforming.

My ideas so far:

Since $R_x$ commutes with its adjoint, it can be written using the eigendecomposition $R_x =  U \Lambda U^*$
$R_x$ is PD, $R_x^{-1}$ is also PD since $\lambda_i > 0 \Longrightarrow \lambda_i^{-1} > 0$
Matrix $U$ defines an orthonormal basis and $C$ is full rank. Can view $U^* C$ as a projection onto this basis.

First, this question is probably better for MATH.SE, but I'll give it a shot.

It's been a long long time since I did this stuff.

If $N > M$:

1) $C^*$ has rank M.

2) $R_x^{-1}$ has rank M. (or it wouldn't exist as an inverse)

3) Therefore  $C^*R_x^{-1}$ has rank M, since $R_x^{-1}$ is a full rank square matrix.

4) $C$ has rank M

5) $C^*R_x^{-1}C$ has rank less than or equal to M and is NxN.

6) $C^*R_x^{-1}C$ is a not a full rank square matrix and is thus not invertible.

If $N \le M$:

1) $C^*$ has rank N.

2) $R_x^{-1}$ has rank M. (or it wouldn't exist as an inverse)

3) Therefore  $C^*R_x^{-1}$ has rank N.

4) $C$ has rank N

5) $C^*R_x^{-1}C$ has rank less than or equal to N and is NxN.

6) $C^*R_x^{-1}C$ could be full rank square matrix and is thus could be invertible.

At this point, I am not sure what is required for $C^*R_x^{-1}C$ to have rank N.  I have narrowed the results, but not fully answered your question.

Gosh, I hope I have this right.

https://math.stackexchange.com/questions/1524444/connection-between-rank-and-positive-definiteness

https://math.stackexchange.com/questions/272049/rank-of-matrix-ab-when-a-and-b-have-full-rank

Starting with another disclaimer:  Ordinarily I don't answer questions here unless I am fairly rock solid in my understanding of the subject matter.  In this case, as the edit history shows, I am muddling about quite a bit.

I had to look up the Frobenius inequality and I understand how you get $ M \leq \text{rank}(C^* R_x^{-1} C) $ from it in the $N \ge M$ case.  I still don't see how it could be greater than $M$.

Oriol B inserted:

  Regarding your question about the Frobenius inequality, I derived it as follows:
  $\text{rank}(C^* R_x^{-1} C) \geq \text{rank}(C^*R_x^{-1}) + \text{rank}(R_x^{-1}C) - \text{rank}(R_x^{-1})=\text{rank}(C^*R_x^{-1}) + \text{rank}(R_x^{-1}C) - M $. And by now using Sylvester's rank inequality, we get $\text{rank}(C^*R_x^{-1}C) \geq \text{rank}(C^*R_x^{-1}) + \text{rank}(R_x^{-1}C) - M \geq rank(C^*) + \text{rank}(R_x^{-1}) - M + \text{rank}(R_x^{-1}) + \text{rank}(C) - M - M = \text{rank}(C^*) + \text{rank}(C) - M = M $
  and the last inequality follows by the assumption.

I do feel confident in my $N>M$ argument, but thanks to this addition, you can change #5 to saying the rank of $C^* R_x^{-1} C$ is $M$.  However, if $ N > M $, then the rank is still smaller than the dimension and it is still not invertible.  It should be with $N=M$.

I'm still mucking with the $N \le M$ case.  Using Frobenius you get $ \text{rank}(C^* R_x^{-1} C) \ge 2N - M $.  For $N=M$ it looks like the result will have rank N and be NxN so it will be invertible.

Since the $N<M$ case is of no interest to you I am stopping for a while on this.

https://artofproblemsolving.com/community/c364309h1480887_rank_inequalities_and_some_consequences

Another update:

I think I found the first flaw in your argument.  I mean, there has to be one in order for mine to be correct.  ;-)

$C$ is $M$x$N$ and $y$ is $N$x$1$ so $Cy$ is $M$x$1$.

$Cy$ can be interpreted as the linear combination of $N$ $M$ length vectors (the columns of C).  If $N>M$, your assertion that Cy can't be zero is false.

I'm not sure where your misconception arises.  Here is a concrete counter example:

Suppose $M=2,N=3$

$$
\begin{bmatrix}
1 & 0 & 1 \\
0 & 1 & 1
\end{bmatrix}
\begin{bmatrix}
1 \\
1 \\
-1
\end{bmatrix}
=
\begin{bmatrix}
0 \\
0
\end{bmatrix}
$$

$$ Cy = 0 $$

Clearly, $y \ne 0$ yet $Cy$ is.  The rank of $C$ is 2.

Response to OriolB's "Of course your assumption ..." comment:

By "your assumption" you mean my second case of $N \le M$.  I didn't get that, sorry.

It seems that you have removed the "less than or" part of step 5, thus proving it invertible.  I had left it hanging since you had said you were only interested in $N>M$ which I proved was not invertible in the first case.

For the $N>M$ case, perhaps a pseudo-inverse will suffice for your greater purposes.


How to detect lines arranged in a visual pattern in an image? I have a black and white image like this:

Visually, there is a clear pattern of black lines evenly separated by somewhat thicker white lines. The problem I need to solve is to automatically detect the black lines and generate the corresponding vector lines for each. What software/library/algorithm can I use to accomplish this?

If it is a true bw-image, i.e. it contains binary values only, you can use the Hough transform for circles (or a special one for ellipses). In the Hough space you will find maxima for the coordinates of the center and the radius or radii, respectively. If necessary, you can use a statistical test to proof the maxima since there are offen multiple local maxima.


C++ based FIR filter design using GNU Octave remez function I am designing a system where I want a FIR filter coefficient generator feature.
I try the remez.cc source file provided in the Octave Signal package.

The source code for remez.cc is here. I have removed the octave interface part (line 757 - end) and the library header (line 34) and created a simple header with only remez function at line 592 then compiled the source code into a static library.

#pragma once
// remez.h
/********************
* remez
*=======
* Calculates the optimal (in the Chebyshev/minimax sense)
* FIR filter impulse response given a set of band edges,
* the desired response on those bands, and the weight given to
* the error in those bands.
*
* INPUT:
* ------
* int     numtaps     - Number of filter coefficients
* int     numband     - Number of bands in filter specification
* double  bands[]     - User-specified band edges [2 * numband]
* double  des[]       - User-specified band responses [numband]
* double  weight[]    - User-specified error weights [numband]
* int     type        - Type of filter
*
* OUTPUT:
* -------
* double h[]      - Impulse response of final filter [numtaps]
* returns         - true on success, false on failure to converge
********************/
int remez(double h[], int numtaps,
          int numband, const double bands[],
          const double des[], const double weight[],
          int type, int griddensity);

Then I followed an Octave example for filter design and implemented a C++ main function for experimenting how to call the remez function. Here is my main function.

#include <remez.h>
#include <cmath>
#include <vector>
#include <iostream>

using namespace std;

int main()
{
    int numtaps = 52;
    vector<double> h(numtaps);
    vector<double> bands = { 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 };
    int numbands = bands.size() / 2;
    vector<double> des = { 1, 1, 0, 0, 0.5, 0.5, 0, 0, 1, 1 };  // designed response for each band

    double ripple1 = 1 - pow(10, -0.3 / 20);
    double att2 = pow(10, -60 / 20);
    double ripple3 = (1 - pow(10, -0.2 / 20)) * 0.5;
    double att4 = pow(10, -70 / 20);
    double ripple5 = 1 - pow(10, -0.1 / 20);
    vector<double> weight = { 1 / ripple1, 1 / att2, 1 / ripple3, 1 / att4, 1 / ripple5 };

    int err = remez(h.data(), numtaps, numbands, bands.data(), des.data(), weight.data(), 1, 16);
    for (auto &c : h)
        cout << c << ", ";
    cout << endl;
    cout << "error code: " << err << endl;
    return 0;
}

This function call produce some weird coefficients. All values inside h vector are 10e144 something. It is clearly an overflow so I debugged the remez.cc and traced what cause this.

The remez function documentation specified that des array should be of size numband. The des array in remez function is passed to CreateDenseGrid function on line 86 for grid generation. The documentation for CreateDenseGrid on line 73 indicate that its parameter des should be size of 2*numband. Clearly since my des array is only of size numband some out of bound memory access happened on line 110.

I changed the des vector to be {1, 1, 0, 0, 0.5, 0.5, 0, 0, 1, 1} and tries again. The coefficients looks nicer, but when I plot the frequency response in Octave it is not anything close to my intended filter design.

I get this weird frequency response:

Anyone could help me with this issue? Is anything wrong with my remez function call? Any suggestion is appreciated here. Thank you very much.

By the way. I am using vs2015, windows7 64-bits, Debug X86 build for compiler settings.

There are some known bugs (such as bug #38134 and bug #48078 on savannah) with the current implementation of remez found in Octave's signal package (see remez.cc on sourceforge), not to mention that the source to link to seems to lag a little from that latest version (you are thus lacking some critical fixes that have been included to address various stability and convergence issues).

Since you appear to be reasonably comfortable with modifying sources, you may want to try out the latest remez.cc source with the submitted patch attached to bug #38134, or alternatively the newly submitted firpm replacement to remez.


Realisation of a passive RC low pass filter I have an unusual low-pass filter problem that I would like to realise with a passive RC network (Note 1).

It's easy enough to find prototype filters for passive RLC networks, or for the various active RC filters.

I am having difficulty finding prototype filters for the classic filter types (e.g. Butterworth, Chebychev, Elliptic) that are implemented using such a network.  The closest I have managed to find is this which recommends equal R and C values for each passive in the network:

Hence, I would like to ask

Is there a good reference that covers this sort of design?
Are there any "tricks" or "gotchas" for this design?
What is the design procedure to implement one of the classic filters?

Note 1 - Context:

The design is for a second-order low-pass anti-aliasing filter. A first order doesn't attenuate an annoying periodic glitch enough - and I may yet need a higher-order filter. Board space is very limited and the ADC has a built-in buffer that I would like to make use of.  The attenuation is desirable - the signal is larger than the ADC would normally accept.  Inductors are completely undesirable for the application.

Note that your desired passive RC network is really just the cascade of multiple first-order RC filters, each of which have the transfer function $\frac{1}{1+RCs}$.  By choosing the R and C values to use in each stage, you can place pole locations where you want them in the $s$-plane.

Here's the problem: recall that the pole locations occur where the denominator of the transfer function is equal to zero. Therefore, for a particular RC stage, its corresponding pole location is:

$$
s = -\frac{1}{RC}
$$

Note that this value is purely real. This makes such a topology infeasible for implementing a prototype such as the Butterworth filter. Its poles are located in a semicircle in the left half of the (complex) $s$-plane. Since you can only place poles on the real axis using your RC stage topology, you cannot implement a Butterworth network. You end up having the same problems with other types like Chebyshev or elliptic structures.

The solution is to use additional elements like inductors or amplifiers; these provide more degrees of freedom in the resulting transfer function to place poles at other places in the $s$-plane. That's why all the references you can find use them!


Voice masking for removing voice and not birdsong First the problem I'm tyring to solve:
I'd like to find a way of removing any potential noise from audio files recorded in noisy environments where we are trying to record birdsong. This is for the sake of privacy, we cannot record anything unless we're sure we aren't recording private conversations.

So, what is the potential for voice masking in this problem? I've done research but it is hard to understand the real potential without learning a lot of background science first. So I'm hoping someone in the community can help.

Basically I'd like to know if voice masking techniques are able to be used on a random sound file to mask all potential voice so it is not understandable... at the same time as ensuring that we aren't masking other sounds, especially birdsong.
Is this doable?

Sound masking in time-frequency is a very complex problem, especially to do it in real-time while collecting data. For ensuring privacy I would suggest using a simpler approach.

Sparse random sampling

Instead of recording continuously for hours (potentially picking up whole conversations), record only short segments spaced from each-other in time. They need to be long enough to capture events in question (birdcall/birdsong), but short enough to not be able to capture much coherent conversation.

For instance record 0.1-1 second segments every 10-100 seconds. This can be done on a fixed schedule, or using reservoir sampling.

This approach has the advantage of not introducing any bias in the selection. It reduces the amount of data, and limits analysis to short time windows. It may also meet the privacy protection requirements.

Drop samples with Voice Activity Detection

While recording, consider segments of 1-60 seconds. Run a standard Voice Activity Detection (VAD) algorithm on each of them. If voice is detected, discard the entire segment. Need to use a low threshold for the VAD to get a low amount of false negatives (sample actually has speech).

Sampling only when the VAD does not trigger may introduce bias in the data. For instance if there is humans talk at certain times of day. Or if some birdsong is confused as human voice by the VAD.


What is unit of time average of product of two signals? Consider the signals $e(t)$ and $h(t)$ with units $\left[\frac{V}{m}\right]$ and $\left[\frac{A}{m}\right]$ and their Fourier transforms $E(\omega)$ and $H(\omega)$ with units $\left[\frac{Vs}{m}\right]$ and $\left[\frac{As}{m}\right]$. The time avaverage P of the product of the two signals (~power spectral density), $e(t)h(t)$, in function of the Fourier transforms $E(\omega)$ and $H(\omega)$ can be written as
\begin{align}
P &= <P'(t)> = \lim\limits_{T\to \infty} \frac{1}{T} \int\limits_{-T}^{T} e(t)h(t) dt \\
&= \lim\limits_{T\to \infty} \frac{1}{T} \int\limits_{-T}^{T} \int\limits_{-\infty}^{\infty} E(\omega)e^{i\omega t} d\omega  \int\limits_{-\infty}^{\infty} H(\omega')e^{i\omega' t} d\omega' dt \\
&=  \int\limits_{-\infty}^{\infty}  \int\limits_{-\infty}^{\infty}E(\omega) H(\omega')\left[ \lim\limits_{T\to \infty} \frac{1}{T} \int\limits_{-T}^{T}e^{i(\omega+\omega')t} dt \right] d\omega d\omega'  \\
&=  2 \int\limits_{-\infty}^{\infty}  \int\limits_{-\infty}^{\infty}E(\omega) H(\omega')\delta(\omega'+\omega) d\omega d\omega' \\
&=  2 \int\limits_{-\infty}^{\infty}  E(\omega) H(-\omega)d\omega \\
\end{align}
The question: the unit of the time average $P$ should be $\left[\frac{VA}{m^2}\right]$ whereas the unit of the final result is $\left[\frac{VA}{m}s\right]$. How can this be and how to interpret it?
Example: if you know the Fourier transforms $E(\omega)$ and $H(\omega)$, and you want to calculate the average $P$, then you will use $P=2 \int\limits_{-\infty}^{\infty}  E(\omega) H(-\omega)d\omega$ but this will give another unit as expected?! So, then you have the wrong result?

What is unit of time average of product of two signals?

It's the product of the units of the signal. In my example is $V \cdot A = W$ which is power and in yours it's  $V/m \cdot A/m = W/m^2$ which is Intensity.

The spectrum of the product and the product of the spectra are completely different things.
Let's look at a simple example. Let's look at a cirucit with a voltage source and single impedance $Z$. For simplicity we will use the complex admittance $Y = 1/Z = Y_m \cdot e^{j \varphi}$ where $Y_m$ is the magnitude and $\varphi$ is the phase of the admittance.
Let's drive this with a voltage $u(t) = u_0 \cdot cos(\omega t)$. The current waveform will then be
$$i(t) = u_0 \cdot Y_m \cdot cos(\omega t + \varphi) $$
The instantaneous power is
$$p(t) = u(t) \cdot i(t) = u_0^2 \cdot Y_m\cdot [cos(\omega t) \cdot cos(\omega t + \varphi)]$$
Using the trig identitiy we get
$$p(t) = u_0^2\cdot Y_m/2 \cdot [cos(\varphi)+cos(2 \cdot \omega \cdot t + \varphi)]$$
The time averaged power would simply be
$$<p(t)> = u_0^2\cdot Y_m/2 \cdot cos(\varphi)$$
The spectrum of the instantaneous power has two frequency compontents: One at $0Hz$ and one at $2 \omega$. The $0Hz$ component is the time averaged power and the actual physical power that would be dissipated in the impedance. It obviously depends on the phase of the impedance. For example for a capacitor or inductor, the dissipated power is 0 and the instantaneous power spectrum ONLY has a reactive componet at $2 \cdot \omega$.
Note that the spectrum of the instantaneous power is a VERY different from the power spectrum of both voltage and the current , i.e. only one component at $\omega$
The assumption that you can get the spectrum of the power my muliplying the spectra of the field quantities is just plain wrong. The math doesn't work out either: multiplication in the time domain is CONVOLUTION in the frequency domain, not multiplication. If you convolve the units work out as well
$$P(\omega ) = \int_{s}^{}U(s) \cdot I(s-\omega)ds$$
which has the correct units of $W/Hz$.


Sample text to collect all possible English biphones for Text-To-Speech Since the list of phonemes in English is fixed, it should be possible to come up with a sample text(s) to collect all possible biphones for text-to-speech synthesis.

Does anyone have a sample text for this purpose?

The "Harvard sentences" corpus is a classic:

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


Need help finding the frequency response So I have this system:
$$
y(n) = \frac{1}{N+1}\sum_{k=-N}^{N}\left(1-\frac{|k|}{N+1}\right)x(n-k)
$$

Now for the impulse response I'm getting:

\begin{cases}\frac{1}{N+1}\left(1-\frac{|k|}{N+1}\right) &\text{for}\quad|n|\leq N\\ 0& \text{otherwise}
\end{cases}

But now I'm supposed to find the frequency response and supposedly a closed form expression can be found for it.

So the frequency response is just the Fourier transform of the impulse response so I have this function:

$$
F(s)=\frac{1}{p(N+1)}\sum_{k=-N}^{N}\left(1-\frac{|k|}{N+1}\right)e^{\frac{-2\pi isn}{p}}$$

From here I used Euler's formula and the $isin$ parts cancel themselves and then I was able to use Legendre's summation formulas to find a super messy closed form but I think that can't be right.  I think instead I have set up the frequency response wrong.  Could some one help me?  Thanks.

Is this for homework?

If so I won't directly give an answer, just a hint about the general approach that can be used to solve such questions without a single derivation:

Plot the impulse response. How does it look? Try building this by combining or modifying (through addition, dilatation, translation, multiplication, convolution, derivation, integration...) very simple functions the DFT of which is known (complex exponentials, constants, rectangles, diracs, Heaviside step). From then, use identities on the DFT to get the DFT of your impulse response.


Combining Bandpass filters on 1D signal to obtain better SNR? I have a signal which I'm trying to extract information from (specifically, features to recognize using DTW). Parts of it have very clear and unique features which are unfortunately only sometimes superimposed with a higher frequency low power interference. Below is an example of a square wave feature.

In this picture,

Black line represents the original signal.
Blue line is a lowpass Butterworth filter (I have added 3 to the value so that the graph looks more legible).
Red line is the signal minus the result of a highpass filter (again offset by +5). And the
green line is a simple box smoothing filter (specifically [1,3,5,7,9,11,9,7,5,3,1]) on the original signal.

The blue line in the second graph is the high pass filter's result.

My questions are: is there any way to disentangle what appears to be a very crisp square wave from a very periodic hum? This hum isn't mains hum, by the way, and it isn't always the same strength or frequency. When it comes on, it'll remain constant all along, but the next time it comes on, it might be a different frequency.

I'm a bit disappointed that I can't get a result that's better than my box filter without severely distorting the resulting signal into something that barely resembles a square signal. Already the low pass filter result (blue graph) is starting to move well past square.

I've gone so far as to thinking of picking the dominant high end frequency, and creating an artificial sine wave to try and "manually remove"  the hum. Is this completely misguided? And what hope do I have of getting a clear signal if I won't have the luxury of fine tuning my filters at design time?

I appreciate any general recommendations.

Edit 2:

Here's a functional equivalent to my problem in python. Any comments as to why I'm not obtaining a clean output signal? I've played around with the cutoff values but I'm not very familiar with signal processing so I'm not sure if I'm approaching it right.

length = 1200
frequency = .00259
interference = 0.0237
interference_amplitude = 0.2

def iirfilter_filter(data):
        nyq = 1/2
        bottom = interference/nyq
        top = interference*10/nyq
        b, a = iirfilter(10, [bottom,top], rs=60, btype="bandstop" , analog=False , ftype="cheby2" )
        return lfilter(b,a, data)

time_series = np.linspace(1, length, length )

sine_wave = np.sin(2 * np.pi * frequency * time_series )
square_wave = (sine_wave / abs(sine_wave))

sine_wave = sine_wave + (interference_amplitude * np.sin( 2 * np.pi * interference * time_series ) ) - 1
square_wave = square_wave + (interference_amplitude * np.sin( 2 * np.pi * interference * time_series ) ) + 1

noisy_sine_wave = sine_wave + ( numpy.random.rand(length) * interference_amplitude)
noisy_square_wave = square_wave + ( numpy.random.rand(length) * interference_amplitude)

a = plt.subplot(2,2,1)
b = plt.subplot(2,2,2)
c = plt.subplot(2,2,3)
d = plt.subplot(2,2,4)

a.plot( sine_wave )
a.plot( square_wave )
b.plot( iirfilter_filter( sine_wave ) )
b.plot( iirfilter_filter( square_wave ) )
c.plot( noisy_sine_wave )
c.plot( noisy_square_wave )
d.plot( iirfilter_filter( noisy_sine_wave ) )
d.plot( iirfilter_filter( noisy_square_wave ) )

plt.show()

Output:

If you know that the interference is only in a given frequency range, you might use a notch filter to filter out this interference. If it is in several frequency ranges, you can either filter with different filters in parallel and use the output that works best for you, or filter the signal serially through all filters (but this might distort the signal shape a lot, depending on the Q-factors and number of filters). Here I show an example with a single filter and a known interferer frequency.

clear all;

Fs = 1000; % 100Hz sampling frequency
t = -5:1/Fs:5;

f_sig = 0.5;

f_intf = 15;
A_intf = 0.05;

% Generate two exemplary signals and pollute them with interference
signal1 = double(sin(2*pi*f_sig*t) > 0);
signal2 = cumsum(signal1 - 0.5) / (Fs/f_sig);
noise1 = A_intf * sin(2*pi*(f_intf)*t);
noise2 = A_intf * sin(2*pi*(f_intf)*t);

received1 = signal1 + noise1;
received2 = signal2 + noise2;

% create a notch filter to filter out this intreference.
w0 = f_intf/(Fs/2);
bw = w0/1; % 1 is the Q-factor here, you can play with it to see the effect
[num, den] = iirnotch(w0, bw);

filtered1 = filter(num, den, received1);
filtered2 = filter(num, den, received2);

subplot(2,2,1);
plot(t, received1)
title('received 1');
grid;

subplot(2,2,2);
plot(t, received2);
title('received2');

subplot(2,2,3);
plot(filtered1);
title('filtered 1');

subplot(2,2,4);
plot(filtered2);
title('filtered2');

Edit: I have ported your code to python, using the correct filter. Here's the filter design function that was taken from scipy v0.19 (not yet officially released), you can find the source here: Scipy Source

def _design_notch_peak_filter(w0, Q, ftype):
    """
    Design notch or peak digital filter.
    Parameters
    ----------
    w0 : float
        Normalized frequency to remove from a signal. It is a
        scalar that must satisfy  ``0 < w0 < 1``, with ``w0 = 1``
        corresponding to half of the sampling frequency.
    Q : float
        Quality factor. Dimensionless parameter that characterizes
        notch filter -3 dB bandwidth ``bw`` relative to its center
        frequency, ``Q = w0/bw``.
    ftype : str
        The type of IIR filter to design:
            - notch filter : ``notch``
            - peak filter  : ``peak``
    Returns
    -------
    b, a : ndarray, ndarray
        Numerator (``b``) and denominator (``a``) polynomials
        of the IIR filter.
    """

    # Guarantee that the inputs are floats
    w0 = float(w0)
    Q = float(Q)

    # Checks if w0 is within the range
    if w0 > 1.0 or w0 < 0.0:
        raise ValueError("w0 should be such that 0 < w0 < 1")

    # Get bandwidth
    bw = w0/Q

    # Normalize inputs
    bw = bw*np.pi
    w0 = w0*np.pi

    # Compute -3dB atenuation
    gb = 1/np.sqrt(2)

    if ftype == "notch":
        # Compute beta: formula 11.3.4 (p.575) from reference [1]
        beta = (np.sqrt(1.0-gb**2.0)/gb)*np.tan(bw/2.0)
    elif ftype == "peak":
        # Compute beta: formula 11.3.19 (p.579) from reference [1]
        beta = (gb/np.sqrt(1.0-gb**2.0))*np.tan(bw/2.0)
    else:
        raise ValueError("Unknown ftype.")

    # Compute gain: formula 11.3.6 (p.575) from reference [1]
    gain = 1.0/(1.0+beta)

    # Compute numerator b and denominator a
    # formulas 11.3.7 (p.575) and 11.3.21 (p.579)
    # from reference [1]
    if ftype == "notch":
        b = gain*np.array([1.0, -2.0*np.cos(w0), 1.0])
    else:
        b = (1.0-gain)*np.array([1.0, 0.0, -1.0])
    a = np.array([1.0, -2.0*gain*np.cos(w0), (2.0*gain-1.0)])

    return b, a

Here's your corrected code using this function:

length = 1200
frequency = .00259
interference = 0.0537
interference_amplitude = 0.2

def iirfilter_filter(data):
        nyq = 1/2.
        bottom = interference/nyq
        top = interference*10/nyq
        #b, a = iirfilter(10, [bottom,top], rs=60, btype="bandstop" , analog=False , ftype="cheby2" )
        b, a = _design_notch_peak_filter(interference*2, 1, 'notch')
        return lfilter(b,a, data)

time_series = np.linspace(1, length, length )

sine_wave = np.sin(2 * np.pi * frequency * time_series )
square_wave = (sine_wave / abs(sine_wave))

sine_wave = sine_wave + (interference_amplitude * np.sin( 2 * np.pi * interference * time_series ) ) - 1
square_wave = square_wave + (interference_amplitude * np.sin( 2 * np.pi * interference * time_series ) ) + 1

noisy_sine_wave = sine_wave + ( numpy.random.randn(length) * interference_amplitude)
noisy_square_wave = square_wave + ( numpy.random.randn(length) * interference_amplitude)

a = plt.subplot(2,2,1)
b = plt.subplot(2,2,2)
c = plt.subplot(2,2,3)
d = plt.subplot(2,2,4)

a.plot( sine_wave )
a.plot( square_wave )
b.plot( iirfilter_filter( sine_wave ) )
b.plot( iirfilter_filter( square_wave ) )
c.plot( noisy_sine_wave )
c.plot( noisy_square_wave )
d.plot( iirfilter_filter( noisy_sine_wave ) )
d.plot( iirfilter_filter( noisy_square_wave ) )

plt.show()

And the output:

Note that you cannot filter the white noise from randn (corrected from rand), as it contains all frequencies and the notch filter only removes a single frequency.