texts: Best Approach to Rank Complex Magnitude Comparision Problem

texts:
This is in reference to the responses for an efficient algorithm for the comparison of bounded fixed point complex numbers at this post.

Efficient Magnitude Comparison for Complex Numbers

See the details of that post for the objectives of the problem. I am now determining my approach to ranking the algorithms to determine who best met the objectives I was seeking, and welcome debate on the ranking approach before I dive in.

Key qualifying factors:

I will baseline a standard CORDIC approach (rotate both vectors to real axis and compare absolute magnitude), as well as what can be done with the use of multiplier operations. The algorithm must be more efficient than either of these approaches (using the equivalent score for the multipliers - see below).

Algorithm must be 100% correct for for magnitude differences within $|z_2- z_1| \gt e$ for any e, where Where $z_n = \sqrt{I_n^2 + Q_n^2}$ with I and Q as bounded signed integers and e is any positive real number>0. (it is understood that it will likely take more operations as e decreases; in fact it would be attractive to be more efficient for large e). The CORDIC is a good example of this, you can select any error bound at the expense of the number of iterations needed.

Acceptable answers may return incorrect results for $|z_2- z_1| \le e$, but a bonus score is included for implementations that provide equivalence resolutions given by the following definitions, with a higher score for tight equivalence:

Loose Equivalence:

$|z_1| \gt |z_2| + e$ returns 1

$|z_1| \lt |z_2| -e$ returns -1

$|z_2- z_1| \le e$ returns 0

Tight Binary Equivalence:

$|z_1| > |z_2|$  returns 1

$|z_1| < |z_2|$  returns -1

Tight Ternary Equivalence:

$|z_1| > |z_2|$  returns 1

$|z_1| < |z_2|$  returns -1

$|z_1| = |z_2|$  returns 0

The function prototype is

result = CompareMagntitudes( I1, Q1, I2, Q2 )

With return values of either $-1,0,1$ corresponding to $<, =, > $ of the first compared to the second (and $0, 1$ for $<, \ge$ for binary solutions).

Test cases will be run with bit ranges from $b = 8$ to $b = 32$ bits for I and Q but nothing in the algorithm should prevent implementation with any size b.

Consider the following closely spaced complex points A, B, C, D as depicted in the diagram below ($A=3+j4$, $B=4+j4$, $C=3+j5$, $D4+j5$).

The true radii are A = 5, B =5.66, C = 5.83 and D = 6.403. In this case, the algorithm should provide a solution to resolve all 4 with 100% confidence (setting e to be $e \le 0.17$ corresponding to the minimum distance between B and C), however it would be acceptable and even beneficial if the algorithm allowed for larger e with an associated efficiency gain.

For example if $e=0.5$ then the following results must be returned using the format $f(z_1,z_2)$ with relation to the function prototype given above:

$f(A,B) \rightarrow -1$

$f(C,A) \rightarrow 1$

$f(A,D) \rightarrow -1$

$f(B,D) \rightarrow -1$

Since all those points have a difference in magnitude from the origin > 0.5.

However the following would be acceptable:

$f(B,C) \rightarrow X$

Where X can be 1, 0 or -1 since B and C have a difference in magnitude from teh origin < 0.5.

The algorithm must be implementable with only the equivalent of standard Boolean operations, binary shifts and compares. Look-up tables if used would add to  buffer size considerations in the scoring.

QUESTION: Please suggest / justify alternative metrics (including alternate scoring the starting numbers I list in my answer, or completely different approaches. It is understood that ultimately there is a trade space and can't be a one size fits all simple score, so consideration is weighted toward the objectives in the original question.


texts: Question about SNR in continuous and discrete domain

texts:
How to make sure the discrete domain's SNR equals to the corresponding SNR of continuous domain?

  If the definition of SNR in continuous domain is:

  $$SNR=\frac{\int_T|s(t)|^2dt}{noiseVariance} $$

  if we sampled $s(t)$ ($T_s$ is the sample period) then what is the
  corresponding SNR of the discrete domain? Is

  $$SNR=\frac{\sum_{n=1}^{N}|s(nT_s)|^2dt}{noiseVariance} $$

  ?


texts:
I have some samples that are read from a square pulse wave and I'm looking for a way to find the characteristics of this wave.

It's a pulse wave, let $T_1$ be the non-pulse duration and $T_2$ to be the pulse duration, usually $T_1$ is bigger than $T_2$.

I would like to know the base frequency, e.g., the frequency that includes $T_1+T_2$.
Can I do this using an FFT?


texts: Image processing and object detection in sales force automation

texts:
We have a specific requirement in our sales force automation product.

Requirement in merchandising module
There is one shelf check activity in merchandising where salesman is capturing which products are available in shelf and how many quantity is available , so he can replenish remaining quantity.

Now we want to automate this process where salesman can capture image of shelf and we want to do image processing where it will check which products are available in shelf and how many quantity available.

We found below 2 ways by which it can be achieved

 - Image Processing though OpenCV in android
 - Image Processing through back end using python which we will send image to
   server after capturing for image processing

Please suggest me which way is better to proceed

Regards,
Bindal


texts: Phase correction algorithm: minimisation phase error

texts:
Thank you for your participation in my discussion in advance.
I am working on the implementation of a phase correction.
Currently, I have finished the algorithm. I think it is a standard algorithm everyone of you uses.
Actually, it is the standard algorithm to correct a phase; I think.

Does anyone have an idea of how can I minimize the difference between angles of adjacent samples/symbols?

Had someone already an experience with such simulation and probably had still the simulation script and can share it with me?


texts: Use polyphase representation to construct a perfect reconstruction filter bank

texts:
I'm posting this question in response to Question about zero-phase filter bank

This is a homework question for which I found an answer, but I didn't do it the way I was supposed to. I'm wondering what the correct approach is. I've been told to use polyphase representation to solve the problem but I couldn't come up with a way to do it.

Let $H_0(z) = 1 + 2z^{-1} + 3z^{-2} + 2z^{-3} + z^{-4}$ and $H_1(z) = H_0(-z)$. Find causal FIR filters $F_0(z)$ and $F_1(z)$ such that $\hat{x}(n)$ agrees with $x(n)$ except for a possible delay and (nonzero) scale.

I have found a solution with no delay and unity scale (see the linked question) but I didn't use polyphase representation. I tried a few times to work it using polyphase but I kept ending up with IIR solutions. What approach should be used that takes advantage of polyphase?


texts:
How to compute Image gradients?

In X Direction - Kernel : -1 0 1

In Y Direction - Kernel :
-1
0
1

Let's assume my input image is: unsigned char* imgImage;

Could any one please share code. I am newbie to image processing.
Is it possible to do in OpenCV with this same kernel?

Thanks a lot in advance!


texts: How to represent ratio between harmonics in audio spectrum

texts: I was provided a report where the amplitudes of spectral harmonics have been extracted from several audio files. The amplitudes have been extracted using Audacity's spectrum visualizer, where the spectrum is provided in dB. I want to compute the ratio between the amplitude of the fundamental and the amplitude of each harmonic. Does it make sense to perform a ratio in dB or should I use the difference between amplitudes? Thank you.


texts: Frequency filtering with a DFT and meaning of removing complex conjugates

texts:
I've seen many implementations to filter out frequency components of some time domain signal by performing a DFT, zeroing the unwanted frequency bins, and performing the IDFT to get the filtered signal, but more than a few times I've seen the results of the DFT treated as frequency components from 0 to the sampling rate instead of 0 to half the sampling rate then the complex conjugates. Take for example this bandpass implementation by a group at MIT: https://github.com/diego898/matlab-utils/blob/master/toolbox/EVM_Matlab/ideal_bandpassing.m

It's defining the frequency bins and corresponding bandpass filter as

Freq = 1:n;
Freq = (Freq-1)/n*samplingRate;
mask = Freq > wl & Freq < wh;

n being the number of samples and wl and wh being the low and high frequencies respectively.

My first intuition is that that's incorrect, but it makes me wonder if I'm missing something since I've seen it done this way multiple times.

But the more interesting question to me is, assuming it is wrong: what physical meaning, if any, does removing one but not both of these complex conjugates have? For example, if I used the above bandpass implementation to make a low-pass filter which filtered out one half the sampling frequency, I would keep the DC component and all the first complex values. But then I would remove all of their conjugates. While I don't actually think that's the low pass filter intended, it does still have the effect of smoothing out the signal, but how would you describe that? I don't imagine it's still considered a temporal filter in that case.


texts: Is there a Python method to calculate moving RMS of a signal? How to calculate the moving RMS using the discrete signal formula?

texts:
I'm trying to calculate a moving RMS of an acceleration signal.
One of the answers in the post, shows the formula for calculating moving RMS for discrete signal.
The formula they suggest:

I have a few questions on this formula. Since $i$ would be, for most part of the summation, less than $N$(no of samples), therefore k would be negative and the negative index does not make sense to me. That's why I'm having trouble understanding this formula.
Also, is there a python method to find moving RMS of a signal like MATLAB has dsp.MovingRMS?


texts: Measure resonant frequencies from system output data

texts:
Suppose I have a real (physical) dynamical system with one capacitive sensor, and I also have the ability to extract data for this system in csv form. The system data output is shown below (graphs) for various sampling frequencies.
As the system is real time and its output depends on the immediate environment (temperature, pressure, etc..) I need to be able to determine if the sampling frequency is hitting any harmonics on-the-fly. All I have is data.

How, in general, can I determine resonant frequencies from data?


texts: Can the overlap-add method be used with varying cut-off frequency filter?

texts:
I'm writing a software analogue synth. I have been using the Overlap-add method to cheaply implement a low-pass filter using FFT.

My low-pass filter has had a constant impulse response; a constant cutoff frequency. Now I'd like to make it possible to vary the cut-off frequency while a note is playing (this is a common feature of analogue synths).

I believe this is impossible while still using the circular convolution theorem trick (that let's us use FFT to implement FIR-filters).

Am I right in thinking so?
If so, is there some other trick I can use to efficiently implement a band-pass filter whose cut-off frequency varies ("sweeps")?


texts: I have a 4096 sequence of FFT coefficient. All of which a real values. My sampling frequency is also 4096 samples per seconds. Now I want to make a plot. I know my Y axis would be my FFT magnitude values, but what would my x- axis be. Would it be my sampling frequency, that is 1,2,3,4.........4096?


texts: Integrator and the multiplier in the Optimal Binary Detector

texts:
The optimal binary detector for asymmetric, equiprobable PAM pulse normally have structure as shown below,

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

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


texts: What's the best way to compare the same object from different "crops" or angles within multiple images?

texts:
I read a similar question, but I'm not sure the techniques discussed will help me.

Take for example the two following images:

I'm trying to determine if the above images are sourced from the same image as well as whether the objects in the following image are the same or similar enough to safely consider them the same.

Thanks in advance
Joe


texts:
R  has a function called filtfilt that performs zero-phase digital filtering by processing the input data in both the forward and reverse directions.  But it suffers from having steep transition bands. Specifically, the output signal exhibits excessive transient response.  The version of filtfilt in Octave apparently handles initial conditions better in a way such that startup and end transients are minimized.

It would be great to have the Octave version of filtfilt coded in R.  Someone has already tried, but it generated an error (https://dankelley.github.io/r/2014/02/19/filtfilt.html). Has anyone attempted to code this? Using the same data, results should be identical in R as in Octave


texts:
In a video sequence, I would like to track an object. Right now, my detection algorithm can correctly locate the object in approximately 10% of the frames, and incorrectly find something else in 5% of the frames.

The incorrect detections are scattered around in the image, while my object movements are constrained by physics.

I've tried implementing a simple tracker based on distance from last detection (frames and pixels), but since my object moves in a 3d world, the velocity can vary a lot.

I think there is a rather simple approach for tracking this, but I don't have a lot of prior tracking knowledge. Where would I start?

(I appologize if this is not for DSP, please redirect me in that case)


texts:
Black rectangles on these charts show same signal, which is actually color sequence from a image, taken from different positions.

It is easy to see that all images contain same pattern. Is here some algorithm which can detect these similarities, ideally implemented in Python ?

What about cut signals in multiple windows and try computing cross corelation on them ?


texts: Is there an algorithm in Matlab that could track moving objects that present change of color?

texts:
I'm studying the photomechanical behavior of photochromic thin films. They react to specific light wavelengths exhibiting bending movement. I use a yellow light, which turns the film from blue to white. As it turns white, the light reflection on the thin film surface intensifies. I want to track this movement through timelapse videos for that I've been doing all the processing in MatLab. To track the thin film movement on the video, I've been using the tracker function found in the Computer Vision Toolbox. So far, I'm not getting exact results due to the change of light intensity, because the algorithm finds it difficult to track the points at each frame (sometimes the tracking points disappear).
Is there a way I can process the video in such a way I can balance the light intensity in each frame?
Is there another MatLab function I can use for my application?
I attached the first and last frame on the video to show how the light changes on the thin film surface:
First frame
Last frame


texts:
How to calculate the variance of noise samples modeled as follows:

$n_a(t)$ is a
Gaussian zero-mean white noise process with (two-sided) power spectral density $\frac{N_0}{2}$. $n_a(t)$ is passed through an anti-alisasing filter to get the output $n_f(t)$. How to calculate the variance of the noise samples $n[j]$ in terms of $N_0$ and $B$, where $n[j]$=$n_f(jT_s)$ and $T_s$ is the sampling period.


texts: What is the error rate in compressed sensing?

texts:
Let $x \in \mathbb{R}^n$ be a $k$-sparse vector. Given $A \in \mathbb{R}^{m \times n}$, we have a measurement vector $y$ given by

$$y=Ax$$

Let $\hat{x}$ be defined as follows

$$\hat{x}:=\arg\min_{z\in\mathbb{R}^n}\|z\|_1,\ \mbox{s.t}\ Az=y$$

Then, which of these two would be called error rate?

$\|\hat{x}-x\|_1$, or,
$\dfrac{\|\hat{x}-x\|_1}{\|x\|_1}$.


texts: Matlab: What should be the BER performance for BPSK using Constant Modulus Algorithm equalizer

texts:
I need help in plotting the Bit error curve or the symbol error curve for BPSK modulation scheme for varying Signal to Noise ratios or Eb/N0. The plot should show the simulated versus the theoretical curve, but I cannot figure out how to mitigate the problems when using the Constant Modulus Algorithm as an Equalizer which are

(1) I want to use the filter function instead of conv in order to model a moving average channel model, chanOut = filter(ht,1,s). But, when I use filter, I am getting an error. How can I use filter function here?

(2) Bit error rate calculation

I am still unable to use filter and unsure if BER curve is proper or not.
Below is the code I wrote:

UPDATE after the answer posted, I have included all the suggestions. The error has been mitigated but now the curve for CMA is coming out to be flat instead of the one shown in the answer. Would like to know what I should do to get similar curve.


texts:
I am detecting a contour of a smooth metal. I do binarization first then use openCV contour function. But sometimes the metal has defects such as dots or dusts at edge. This leads a unsmooth edge. Is there any way I can detect these defects with openCV function or some math calculation?
This is the contour result.

This is zoomed images for defected contour result.


texts: Cross correlation and cross power spectrum for signal + noise

texts:
Quick question about my understanding of cross correlations and cross spectral densities.

Let $C_{ab}$ be the cross correlation between the signals $a$ and $b$, i.e., $C_{ab}(\tau) = \langle a(t)b(t+\tau)\rangle$. Let $S_{ab}$ be the cross spectral density of the signals $a$ and $b$, i.e., $S_{ab}=\mathcal{F}(a)^* \mathcal{F}(b)$ and by the convolution theorem/cross correlation theorem $S_{ab} = \mathcal{F}(C_{ab}(\tau))$.

Say $x(t)$ is some signal, and $\nu(t)$ is some uncorrelated noise signal. Let $s(t)=x(t)+\nu(t)$
Thus, $C_{xs} = C_{xx}+C_{x\nu}$ from the linearity of the expectation value being calculated in the cross correlation. By initial assumption that $x$ and $\nu$ are uncorrelated, we get $C_{xs}=C_{xx}$. (This calculation so far also appears here on page 5, for example so I think it is correct).

Now, if we take the fourier transform of both sides of this, we get $S_{xs}=S_{xx}$. Is this correct so far? But does this not mean that $\mathcal{F}(x)^* \mathcal{F}(s) = \mathcal{F}(x)^*\mathcal{F}(x)$? and hence $\mathcal{F}(s)=\mathcal{F}(x)$, and hence even $s(t)=x(t)$? Clearly there is something wrong here. What mistake am I making?

Thanks!

EDIT: On a similar note, if $C_{x\nu}=0$, does that not mean that $\mathcal{F}(x)^*\mathcal{F}(\nu)=$ and hence (by taking the absolute value of both sides and squaring) that $S_{xx} S_{\nu\nu}=0$? This is patently false. At each frequency this would require either the power in my signal to be zero, or the power in the noise to be zero. Even if my noise is white noise then this is clearly false. I'm doing something fundamentally wrong here. What? Any help?


texts:
I'm writing a software that read engine rpm during acceleration, using a recorded audio.
To do that, I divide the record in small time-steps, and for each one I do an FFT, to know the rpm in that timesteps.
Because, in the time-steps, the rpm is not constant, but vary from minimum to maximum (because of the acceleration), I'm doing some test using a pure sine wave with frequency varying from 300 to 545 rad/s (48 to 87 Hz) in 4s time.
In the real audio, I take a time-step of about 0.5s, but now I'm studying what happen with longer time-step to understand how FFT works in variable-frequency audio.

This is the FFT of that signal:

As you can see, the two peaks are at 52 and 125 Hz. Now, I'm wondering why there is 125Hz... I tought I will found more or less two peak at 50-90 Hz.
So, why? And how can I estimate the frequency range during an acceleration?


texts:
I've got images of blood drops in grayscale and I want to generate 3D data where the first two dimensions represent the area of interest and the 3rd dimension represents the pixel intensity of that area. Here's an example image:

My requirement is as shown in D,E,F while the raw images are A,B,C. Also, what algorithms can be used for detection of the plateau areas as detected in G,H,I?


texts:
I am trying to understand as much as possible the aircraft replies of ADS-B signals. It uses a PPM (pulse position modulation) scheme with a carrier frequency of 1090 plus minus 1Mhz. You can find all the important specifications about ADS-B reply signals on this Link, by clicking on annex 10 volume IV page 35-36 if you are interested.
I think I understand almost everything except things related to the "spectrum":

What is the reply spectrum? how does it affect the signals frequency in practice?

How do I interpret the following plot?

I plotted the amplitude alongside the phase and noticed that for a lot of signals of my dataset: there is a clear pattern where if amplitude increase, then phase decrease and vice-versa. Even though for ADS-B signals, the behaviour of the phase is unspecified. Could this be related to the spectrum and potential small frequency changes when the transmitter changes amplitude?


texts: In LTE system SNR is -1 dB and it is below the noise figure. How can you detect the signal?


texts: Kalman Filter: Why $ Q $ Discrete Is Defined as $\int_0^Te^{\mathbf{A}\tau} Q e^{\mathbf{A}^T \tau} d\tau$?

texts: I would like to ask, why in the transformation to the discretization, $\mathbf{Q}$ is obtained from the expression containing the integral (image attached), what is the theory behind it?


texts:
I posted this in stackoverflow but a poster suggested I ask the question here too...

Basically, I'm getting a USB audio stream in the form of a unsigned 8-bit integer pointer. I want to convert whatever bit rate it's using (whether 16 or 24 bit audio) and scale it up to 32 bits. Once it's scaled up to 32 bits, I'm placing it in a microcontroller's peripheral that will output the data in a 32 bit register. I would also like to do this in the reverse order (getting input from a codec in 32-bit words and then sending it to my USB Host in either 16 or 24 bit audio streams).

Just wanted to open a discussion on how people bit scale these streams and what techniques/algorithms they use.

tl;dr - I want to be able to scale up 16/24 bit audio streams to 32 bit words and vice versa for scaling down. I know the process for scaling down as such:

// For 32 to 16 bits
uint16_t adjusted_bit = uint32_bit_data >> 16;
// For 32 to 24 bits
uint32_t adjusted bit = uint32_bit_data >> 8;

Thanks!


texts: Confusion about ensembles and averages in autocorrelation matrices

texts:
I just noticed that until now I often don't cared about the scaling of the autocorrelation Matrix in Matlab. I then checked with the book Statistical Digital Signal Processing from M. H. Hayes. There it is defined like this:

First of all, why does the author note it like $x(0)x^*(0)$? Written like this wouldm't this mean first element of the vector x times first element of the vector x? And the expected value is then what $\frac{x(0)x^*(0)}{2}$?

$r_x(\tau)$ can be calculated as $E[X(t)X(t+\tau)]$ This is more what I would expect. A pointwhise vector multiplication of a $\tau$ shifted vector.

Here comes my actual question. For e.g. $x = [1 ,2, 3]^T$, $r_x(0)$ should become $\frac{1^2+2^2+3^2}{3}=\frac{14}{3}$ in my oppinion. Matlab shows for xcorr(x)

3.0000    8.0000   14.0000    8.0000    3.0000

$r_x(0)$ should be the value in the middle. But why don't they take the expected value? On this page they suggest that the autocorrelation matrix can be calculated with

z=autocorr(x)

and

Rxx=toeplitz(z)

Instead of autocorr, xcorr(x) should work too stated other posts. As expected this results in

Rx =

    3.0000    8.0000   14.0000    8.0000    3.0000
    8.0000    3.0000    8.0000   14.0000    8.0000
   14.0000    8.0000    3.0000    8.0000   14.0000
    8.0000   14.0000    8.0000    3.0000    8.0000
    3.0000    8.0000   14.0000    8.0000    3.0000

Again, no expected values. And this solution is also wrong because it is shifted. $r_x(0)$ is not in the upper left corner. I could shift it by hand, but is that how autocorrelation matrixes are calculated in Matlab? I checked with some other sources in the internet. Often the algorithm to calculate an entry of an autocorrelation contains just the sum.


texts: discrete Haar wavelet transform, fast and efficient method?

texts:
I'm working on my own implementation of the discrete Haar wavelet transform, I understand the wavelet theory and how to construct the Haar matrix of size N to perform the transform, but obviously there is a problem using the Haar matrix in application - it's simply too big.

I am working on an application that applies the Haar transform to an audio signal. If I sample the signal at 44.1 kHz, even a one second recording would require a 2^16 by 2^16 Haar matrix to do the transform in one step, which is obviously impractical and a hardware constrained machine such as a phone wouldn't have the capability to hold such a matrix in memory.

The other method I've seen used is that a 2X2 Haar matrix is applied to the entire signal iteratively, and the results are stored in two arrays - one array holding the "average" Haar coefficients (first element of the output vector) and the other holding the "difference" coefficients (second element of the output vector). The process is then repeated over the "average" coefficients - as these coefficients are essentially the result of a lowpass filter.  Each time the process is repeated the number of elements needed to process in the following step is halved, until only one lowpass coefficient is left. This seems fine, and definitely works but when I started to think about it, it just seems that it would be really slow.

My main question is, what's a good way to implement a fast and efficient haar transform? Or a practical way to apply one of these two methods?

PS: I've been learning about this completely on my own, so if I made some mistakes in my explanation or way of thinking about this stuff let me know.

PSS: I've never done any kind of audio processing before, so if you know about some other filters that I should apply to the raw signal before doing a wavelet transform let me know!

PSSS: I know that the Haar wavelet may not be the best for this type of signal processing, but when I tried to learn about using other DWTs such as the Daubechies wavelets, the literature seemed very confusing, or at least was directed at more advanced readers. If anyone could point me in the direction of how to implement other DWTs, that would be great.


texts: Is it correct to call a Savitsky-Golay filter of degree 0 a simple moving average?

texts:
I had this question after seeing that polynomial regressions fit polynomial functions of different degrees to a time-series, since the mean of a time series is a constant and that a constant is also a polynomial of degree 0 i wondered if the SG filter who also fit polynomial functions to a signal over a period window can be called a moving average.

So is it correct to say that a Savitsky-Golay Filter that try to fit a polynomial of order 0 to a signal is considered a moving average filter ?


texts: How to choose the "best" measurment (from a given set) as input for a kalman filter?

texts:
Problem:
For an object tracking scenario with multiple objects and multiple tracks, I want to choose the "best" assignment object<->track. Therefore, I need a parameter indicating how suitable an object is for a track.
(Each object refers to a measurement and each track is a Kalman filter)

Idea:
Given the classical Kalman filter with the matrices state $x$ and measurement transition $H$
We can calculate the residual $y_t$ between the predicted measurement $x_{pred}$ and the measurement $z$ by:

$y_t = z-H x_{pred}$

However, this will be a vector. Is there a general approach to reduce this to one parameter? (e.g distance when having just position information). I assume that some kind of normalisation is needed to get all the residuals for each state in the same range? Building the l2-norm would be my first approach.

Question:
Which parameter is suitable to indicate how "good" the measurement fits to the Kalman filter?

Is the residual suitable? If so, how to handle different value ranges, and how to "compress" the vector to just one value?


texts: Practical cross-spectrum estimation using Blackman-Tukey approach

texts:
I would like to estimate the cross-spectrum of two signals using the (lag-windowed) Blackman-Tukey approach but I'm having difficulties with proper practical implementation. As defined in equation 2.8.32 of the the book Spectral Analysis of Signals (2005, Prentice Hall) by Stoica&Moses, the Blackman-Tukey periodogram estimate of the cross-spectrum between two stationary signals $ y(t)$ and $ u(t)$ (for the purpose of my question we may also assume that both signals are real-valued) is given by

\begin{align*}
\hat{\phi}_{yu}(\omega) = \sum_{k = -M}^{M}  w(k) \ \hat{r}_{yu}(k) e^{-i\omega k}
\end{align*}

where $ w$ denotes some appropriate lag-window (symmetric with $ 2M'+1$ non-zero observations, $ M' < M$ ) and $ \hat{r}_{yu}$ estimated cross-covariance function with M lags (also in total $ 2M+1$ observations). I try to implement the estimation via Discrete Fourier Transform (DFT) using a typical Fast Fourier Transform (FFT) algorithm available in most statistical packages, say MATLAB's fft function (documentation here).

What I am not sure about is how to perform the DFT on windowed cross-covariance values since index $k$ is taking both negative and positive values. When using MATLAB's fft function we feed in a vector representing the signal (with some length $N$) we want to transform and the algorithm performs the transform with index $k$ running from $1$ to $N$. As I understand it, if we were dealing with a DFT of an autocovariance function, say $ \hat{r}_{yy}$, then the negative indices would be no problem as $\hat{r}_{yy} $ is symmetric so the "two sides" around $ k = 0$ are essentially the same and one could just perform an "one-sided" DFT by feeding values $ w(m) \ \hat{r}_{yy}(m) \ , \ m = 0,1,...,M$ to the fft algorithm. As symmetry need not be the case with cross-covariances, my understanding is that this approach does not work. My initial solution was just to cram the windowed cross-covariance values

$$  w(k) \ \hat{r}(k) \ , \ k = -M,...,-1,0,1,...M $$

into a vector of length $ 2M + 1$ and feed it to the fft function as it is. However, I fear there are problems with this approach due to exponents of the DFT taking up only negative values. My question is therefore how should one perform the DFT with a vector of (windowed) cross-covariance values in order to obtain the desired estimate of the cross-spectrum?

Just to be clear I explicitly want to use the Blackman-Tukey / cross-covariance approach and not e.g. Welch's approach in estimation. Help is greatly appreciated!


texts: Why are the Fourier transforms of autocorrelation and cross-correlation different?

texts:
I'm almost definitely misunderstanding this, but as far as I can tell we can equate the Fourier transform of the auto-correlation function $R_{xx}(\tau)$ with power spectral density (PSD) like so:

$\hat{R}(\omega) = S_{xx}(\omega) = \lim_{T \to \infty} \textbf{E}[|\hat{x}(\omega)|^{2}]$

Whereas the Fourier transform of the cross correlation is just the complex product of the two functions in frequency space:

$\mathcal{F}\{x(t)∗y(t)\}=\mathcal{F}\{x(t)\}\mathcal{F}\{y(t)\}^{*} $

Why does the expected value/mean come into play with the auto-correlation/PSD?


texts: Analysis and algorithmic complexity of the Constant Q transform

texts:
Here I am again (Implementation of the constant Q transform + property questions) :)

I was reading this paper on a constant-q transform algorithm, where they managed to turn the following equation for the constant-q transform:

into this:

(where X[k] is the fourier transform of x[n], and K*[k, kcq] is the complex conjugate of the fourier transform of the w[n, kcq])
- claiming it to be an efficient algorithm.

While I'm still new to algorithmic complexity, I tried to derive the big-O complexity. I'm assuming the transforms of the windows (ie. K*[], the kernels) are pre-calculated. Let N be the size of the transform, and K the number of components (kernels / filters / results). Would this be correct?

So thats the FFT (N * log2(N)) plus K loops of N multiplies (the convolution) + additional stuff done at each loop (done K times).
Also, since the kernels are precalculated, that incurs both a certain amount of memory storage needed + a lot of memory access during the calculation. With their implementation, they store the kernels in full size (that is, N), so the required storage memory for the calculation would be K kernels of N size.

I'm defining M() to be the amount of memory used during an calculation:

That is, we need to load the whole signal (N), also K number of N-sized complex (2) kernels. That is quite different from algorithms like the standard fft and direct evaluation of the constant-q transform which requires no working storage, and as of such are defined to be:

-Since they only need to load the signal. If we define the complexity of the FFT and the direct evaluation of the constant-q transform to the following:

(This assumes direct evaluation of the CQ-transform which incurs at least NK operations, together with NK operations calculating the discrete window for each component.)

we can create a small comparison table for an example. If the transform is of size 2048, and the wanted number of components is 1048 (N = 2048, K = 1048), we get the following results:

It would seem the cqfft algorithm (theirs) incurs a rather extreme memory cost (that in practical appliances could dominate computational speed) for a mere ~2x increase in computation speed. Is this correctly derived?

Here's a reference-pseudocode of their implementation i wrote, as i understand it. If i have misunderstood their algorithm, it will probably show here:

std::vector<complex<float>> fft(float[], size_t);

class kTransform
{
    vector<complex<float>> cqfft(vector<float> signal)
    {
        vector<complex<float>> fftdata = fft(signal, size);
        for (int kernel = 0; kernel < kernels.size(); ++kernel)
        {
            complex<float> accumulator;
            for (int k = 0; k < size; ++k)
            {
                accumulator += fftdata[k] * conj(kernels[kernel][k]);
            }
            fftdata[kernel] = accumulator / size;
        }
        return fftdata;
    }

    void compileKernels(vector<float> freqs, size_t dftSize)
    {
        auto numKernels = freqs.size();
        kernels.resize(numKernels);
        size = dftSize;
        for (int kernel = 0; kernel < numKernels; ++kernel)
        {
            kernels[kernel] = fft(Window(freqs[kernel]));
        }
    }
    size_t size;
    vector<vector<complex<float>>> kernels;
};

Thank you for any insight. As a bonus question, does anyone know what the complexity of a fast implementation of the Gabor transform is?


texts: Python: Dynamic Time Warping, what actually is a 'similarity score' and how to compute it?


texts:
What type of sensors have complex Gaussian noise?

Do strain gauges also exhibit such noise?


texts: Why coherence function calculation needs averaging of different segments of the signal?

texts:
I am trying to understand coherence function calculation. I am working with a non-linear system without any noise. Coherence function is defined as:

I understand that coherence is 1 for N=1, but I don't understand why coherence is different with segmentation? As far as understand with segmentation without any noise in system, each segment should give same values and them coherence will go to 1 again. I think I am missing something. I am quite new to this subject.
Could anyone help me with this? Is there any good reference related to this where it has proved that why coherence is different with segmentation?
I have gone through the following questions:
Optimal segment length for coherence estimation
and Magnitude-squared Coherence calculation inconsistence
but I still it is not clear to me.
Thanks,
Swati


texts: What are the advantages, if any, of derivative sampling?

texts:
In Five short stories about the cardinal series $[1]$, the author makes the following comment:

  Interestingly enough, Shannon goes on to mention that other sets of
  data can also be used to determine the band-limited signal--for
  example, the values of ƒ and its first derivative at every other
  sample point, the values of ƒ and its first and second derivatives at
  every third sample point, and so on.

The paper mentions some historical developments, but I'm curious what the "killer apps" are for derivative sampling.  Does it go by any other names?  Are there further generalizations of this approach?

A simple overview, or pointers to some references would be great.

--

J. R. Higgins, Five short stories about the cardinal series, Bull. Amer. Math. Soc. (N. S.) 12 (1985), no. 1, 45-89. http://bit.ly/plioNg


texts: autoregressive moving average code implementation

texts:
I am new to DSP and i am trying to take a wav (human speech) file and apply ARMA filter and plot its PSD graph in python.

I see that there are a lot of AR implementations but almost none ARMA. I found only statsmodel that implemented ARMA but it requires the data to be stationary (but sound with lags of 30 is almost stationary)

Does anyone have a working piece of code that i can use to plot?
If not, i will be happy to get some reference on how to implement it.

Thanks in advance!


texts: Wave Shaper Function From A Waveshaped Singal

texts: Is there a way to work out an unknown waveshaper function, that was used to distort and map the values of an input signal to new output signal, form the output if we only have input signal and the output signal but missing the waveshaper function ?


texts: Help with finishing this integral, to obtain the power spectral density of a pure cosine wave

texts:
I am trying to evaluate the power spectral density $S_{xx}(f)$ of a cosine signal $x(t) = A\cos(2\pi f_0t)$, by starting from its definition for deterministic power signals
$$S_{xx}(f) = \lim_{T\rightarrow\infty}\frac{1}{T}\left|X_{T}(f) \right|^2\;, \tag{1}$$
where $X_T(f)$ is the Fourier transform of a signal $x(t)$ which has been truncated to a time window $T$, given by $$X_T(f) = \int_{-T/2}^{T/2}x(t)\;e^{-2\pi if t}\;dt \;.\tag{2}$$

I know that the answer should be $$S_{xx}(f) = \frac{A^2}{4}\bigg[\delta(f-f_o) + \delta(f+f_o)\bigg]\;,\tag{3}$$ but I am stuck with my working. Can anyone check it so far, and help me fill in the gaps? Perhaps my mathematics is not as good as it ought to be!

\begin{align}
  S_{xx}(f) &= \lim_{T\rightarrow\infty}\frac{1}{T}\left|X_{T}(f) \right|^2\ \tag{4}  \\
            &= \lim_{T\rightarrow\infty}\frac{1}{T}\left[ \int_{-T/2}^{T/2}x(t)\;e^{-2\pi if t}\;dt  \int_{-T/2}^{T/2}x^\ast(t')\;e^{2\pi if t'}\;dt' \right] \tag{5} \\
            &= \lim_{T\rightarrow\infty}\frac{1}{T}\left[ \int_{-T/2}^{T/2}\int_{-T/2}^{T/2} A^2\cos(2\pi f_0t)\cos(2\pi f_0t')\;e^{2\pi if(t'-t)}\;dt \;dt'\right] \tag{6}  \\
            &= \lim_{T\rightarrow\infty}\frac{A^2}{2T}\left[ \int_{-T/2}^{T/2}\int_{-T/2}^{T/2} \Big[\cos\big(2\pi f_0(t+t')\big) + \cos\big(2\pi f_0(t'-t)\big)\Big]\;e^{2\pi if(t'-t)}\;dt \;dt'\right] \tag{7}
\end{align}
where the last line is obtained by using the product-to-sum identity for cosines. Now, I'm a little stuck. I can make the substitution $\tau = t'-t$, which simplifes the exponential

$$  S_{xx}(f) = \lim_{T\rightarrow\infty}\frac{A^2}{2T}\left[ \int_{-T/2-t}^{T/2-t}\int_{-T/2}^{T/2} \Big[\cos\big(2\pi f_0(2t+\tau)\big) + \cos\big(2\pi f_0\tau\big)\Big]\;e^{2\pi if\tau}\;dt \;d\tau\right] \tag{8}  $$

but now I am not sure how to proceed. I can also replace my cosines with complex exponentials to obtain

$$\lim_{T\rightarrow\infty}\frac{A^2}{4T}\left[ \int_{-T/2}^{T/2}\int_{-T/2}^{T/2} \Big[ e^{2\pi if_0(t+t')} + e^{2\pi if_0(t-t')} +e^{2\pi if_0(t'-t)}+e^{-2\pi if_0(t+t')}    \Big]\;e^{2\pi if(t'-t)}\;dt \;dt'\right]$$

Can someone teach me how to finish this? Ideally I would like it if someone could show me how to evaluate it in terms of $T$, and then show how the $\delta$ functions arise as we finally take the limit $T\rightarrow\infty$.


texts: Loudness detection for natural conversation

texts:
I'm working on a project where I want to be able to analyze natural conversation as an audio file and determine the points in the conversation where the speech is louder than normal i.e. be able to classify the points in the audio file that contain louder than normal voice.

I would like to do the same for quieter than normal voice.

Is there any good way to do this?

Thank you!


texts: Gardner Symbol Timing error detector gain

texts:
I am trying to design the Symbol timing recovery block for QPSK receiver.
I have set the loop bandwidth to be 0.005, and damping factor to be 1. I am following Digital Communications: A Discrete-time Approach by Michael Rice, for loop filter and Gardner TED design. According to S-curve for gardner TED, the gain Kp comes out to be around 0.1751. I calculated other parameters such as K0 and K1 for the loop filter. However, the loop is not converging, infact its "mu" value is showing a behavior of converging and then diverging.

I think I have calculated the Kp (gain of Gardner TED) wrong. can somebody guide about the value of Kp and gardner TED recovery approach?


texts: Inverse Z-transform mystic simplification

texts:
I have the following expression:

$$X(z) = \frac{16}{15}\frac{1}{1-\frac14z^{-1}} - \frac{16}{15}\frac{1}{1-4z^{-1}}$$

According to my understanding this should become:

$$x(n) = \frac{16}{15}\left(\frac14\right)^n  u(n) - \frac{16}{15} 4^n  u(n)$$

But according to my source it is:

$$x(n) = \frac{16}{15}\left(\frac14\right)^n  u(n) + \frac{16}{15} 4^n  u(-n-1)$$

Are these expressions equal?

If they are, how can $-u(n) = u(-n-1)$ ?

If not, why?

I used $a^nu(n)\Longleftrightarrow  \frac{1}{1-az^{-1}}$.


texts: Criteria to distinguish FIR and IIR filters from difference equation?

texts:
What criteria should be used to safely decide if it is an IIR filter or FIR filter from a difference equation? FIR filter are always stable, meaning all poles are inside the unit circle AND have pole-zero cancellations i.e have equal number of poles and zeros at the same point in the z-plane. IIR filters are conditionally stable, i.e always have poles on/outside the unit circle. Is this good enough to distinguish from the difference equation? what other criteria should be considered?

PS: Please don't say IIR filters have recursive difference equations, and FIR doesn't. It is not true.


texts: Dealing with Random Errors in Sampling Rates?

texts:
I wonder what generalization I can make about what happens when I introduce some error in my sampler spacing.

For example, If I sample once every second I would want my sampler to take a signal at 1s,2s,3s,4s.

Unfortunately, my real world sampler looks like more like 1.01s, 1.99s, 2.99s, 4.01s.

I want to make a statement about the Fourier domain. I expect my high frequencies to be ambiguous while my lower ones should not to change, too much. Is this correct?

I understand how to do this in terms of a uniform change, for example 1.01, 2.02, 3.03 (which is just a change in sampling rate), but I'm wondering how I can describe this with some knowledge of my error distribution, or for my application even understanding how to describe this knowing the mean of the error would be wonderful.

Can somebody point me in the right direction?


texts:
In the image, the signal is a pseudorandom sequence with values oscillating between -1 and +1. It is periodic (in this case with period 0.001). I don't see how the sequence yields a fourier transform whose real and imaginary parts are given on the right side.
Any help on how to compute this would be really appreciated.


texts: Why does changing the frequency of a single component change total sum of amplitudes in my FFT?

texts:
I have the following code which generates a single simple sine wave with amplitude $1$ and then calculates the Fourier transform, normalises it, and finally takes the sum of all the components.

However I've observed something strange, if I generate a $10\textrm{ Hz}$ wave and a $100\textrm{ Hz}$ wave. The amplitudes are pretty damn similar as expected and they are  0.9999483519908307 and 0.9997854546246567 respectively.

However the total sum of all the components after normalisation are 1.01723486798 and 1.17263231047 respectively, which is approximately a $15\%$ difference! Why is this occurring?

Below are the plots of the spectrums and the python code below.

mport matplotlib.pyplot as plt
import numpy as np
from decimal import Decimal

#number of samples
iNumSamples=10000

#period of samples (i.e. sampling at 10Khz)
Ts=1/iNumSamples

#time axis used to sample 10,000 times for a second
t = np.linspace(0, 1, iNumSamples)

#10Hz Signal
signal1=np.sin(t*20*np.pi)

#100Hz Signal
signal2=np.sin(t*200*np.pi)

#creates a new figure object
freqFig=plt.figure(2)
#change title
freqFig.canvas.set_window_title('Frequency domain')

#fast fourier transform of 10Hz signal
absoluteValues1=abs(np.fft.rfft(signal1))

#fast fourier transform of 100Hz
absoluteValues2=abs(np.fft.rfft(signal2))

#find the frequencies of the rtansform
frequencies=np.fft.rfftfreq(iNumSamples, d=Ts)

##Normalises (rescales) the fft in a sensible way so that the corresponding freq component of a sine wave is 1 if the sine wave amplitude is 1
##non zero freq amplitudes (even when using rfft) are halfed as there are negative and positive components, hence first multiple by 2

#first create empty array of correct size
normalisedValues1=np.zeros(absoluteValues1.size, dtype=Decimal)
normalisedValues2=np.zeros(absoluteValues2.size, dtype=Decimal)

#normalise 0 component
normalisedValues1[0]=absoluteValues1[0]/iNumSamples
normalisedValues2[0]=absoluteValues2[0]/iNumSamples
#normalise non zero components
normalisedValues1[1:]=absoluteValues1[1:]*(2/iNumSamples)
normalisedValues2[1:]=absoluteValues2[1:]*(2/iNumSamples)

print("peak value signal 1:")
print(np.amax(normalisedValues1))
print("peak value signal 2:")
print(np.amax(normalisedValues2))

print("total sum of all componenets in spectrum of 10 Hz")
print(np.sum(normalisedValues1))

print("total sum of all componenets in spectrum of 100 Hz")
print(np.sum(normalisedValues2))

#creates a new figure object
freqFig=plt.figure(1)
#change title
freqFig.canvas.set_window_title('Frequency domain')
plt.plot(frequencies,normalisedValues1)

#creates a new figure object
freqFig=plt.figure(2)
#change title
freqFig.canvas.set_window_title('Frequency domain')
plt.plot(frequencies,normalisedValues2)

plt.show()


texts: Improving spectrogram resolution in Python?

texts:
I'm using the specgram() function in matplotlib to generate spectrograms of speech wave files in Python, but the output is always of vastly inferior quality to what my normal transcription software, Praat, can generate. For example, the following call:

specgram(
    fromstring(spf.readframes(-1), 'Int16'),
    Fs=framerate,
    cmap=cm.gray_r,
)

Generates this:

While Praat, working on the same audio sample with the following settings:

View range: 0-8000Hz
Window length: 0.005s
Dynamic range: 70dB
Time steps: 1000
Frequency steps: 250
Window shape: Gaussian

Generates this:

What am I doing wrong? I've tried fiddling with all the specgram() parameters, but nothing seems to improve the resolution. I have virtually no experience with FFTs.


texts: Evaluation of Jacobian for Extended Kalman Filter

texts:
For the non-additive noise case,
\begin{equation} x_k = f(x_{k-1}, u_{k-1}, \xi_{k-1}) \\
y_k = h(x_k, \nu_k)
\end{equation}

the EKF takes into account the jacobian wrt to the noise terms $ L_{k-1} = \frac{\partial f}{\partial \xi} |_{\hat{x}_{k-1|k-1}, u_{k-1}} $ and $ M_{k} = \frac{\partial h}{\partial \nu} |_{\hat{x}_{k|k-1}} $

I know that jacobian wrt to state: $ A = \frac{\partial f}{\partial x}$ and $ H = \frac{\partial h}{\partial x}$ are evaluated at means of noises i.e. at $\xi = 0, \nu = 0$ resp.

But I'm not clear where to evaluate $L, M$ ? The above expression for $L_{k-1} , M_k$ tells only the $x, u$ but not the noise terms. Should  $\xi = 0, \nu = 0$ or some random sample?


texts: Decimation with IIR filter altering amplitude of signal

texts:
I'm downsampling a signal through decimate (scipy implementation), but the result has a different amplitude than the original signal, which if I understood correctly, shouldn't happen. The original signal is sampled at 128Hz, and I want to downsample it to 1Hz. Decimation code is:

    sig_d = signal.decimate(sig, 8)
    sig_d = signal.decimate(sig_d, 8)
    sig_d = signal.decimate(sig_d, 2)

I then checked in the first 1 minute window of the original and decimated signals:

    df = pd.DataFrame(sig[0:7680])
    print(df.describe())
    df = pd.DataFrame(sig_d[0:60])
    print(df.describe())

and this is the output statistics:

          Original
count  7680.000000
mean     94.191307
std       0.551882
min      92.580971
25%      94.028488
50%      94.086906
75%      94.794744
max      94.969666

       Decimated
count  60.000000
mean   91.011716
std     0.516294
min    89.658212
25%    90.847055
50%    90.886592
75%    91.589404
max    91.689528

As you can see, there is a huge decrease in amplitude for the same time window in the decimated signal. Per the statistics, it doesn't seem to be related to peak removal, since this window in the original signal has very small variability, and not really any pronounced peaks. I did plot the two signals to check if I wasn't messing up the alignment, and they were aligned perfectly, their shape in this time window is the same (besides sampling rate), but the decimated signal seemed to be "pushed down" in amplitude. Unfortunately, I cannot post the plots (for proxy reasons).

This seems to happen only with IIR filter(though FIR has other issues). Is this normal behavior of decimation? Or am I doing something wrong?


texts:
I have been using the Canny edge detection function in OpenCV to detect the edge of an elliptic annulus (light object, dark background). In theory there should be two edges. However, it sometimes does not all edges on the outer ring.

So if I radially plot my detected points it looks like this (x: angle, y: radius from center of mass)

There is a wavy line which is the outer border of my object and a few outliers which are the inner border (you can tell because they look like they are systematically offset). I don't know the exact breadth of the annulus.

Q: Is there any way to programmatically get rid of these "obvious" outliers so that the exact outer contour is interpolated?


texts: Is splitting single white gaussian noise samples-stream in multiple sample-streams still gaussian noises?

texts:
Suppose I have a random generator that generates gaussian distributed samples (mean = 0, variance = sigma) at a rate of 4*N samples per second.

I want to generate 4 different streams from this sequence (each at a rate of N samples per second) by taking 1 out 4 of the generated samples.

Will each stream still be gaussian distributed samples with mean = 0 and variance = sigma ?

NB: I want to see if I can save resources on a FPGA board by creating only a single noise generator to create multiple streams (at the expense of having a higher clock rate) for this generator.


texts: What is the unit for power spectrum of acceleration signals in log scale

texts:
I have a signal recorded by an accelerometer in m/s^2. It is basically accelerations over a period of time. I have calculated the power spectrum of the signal over in a certain frequency band and then transferred into a 20log(x) scale for plotting. I am wondering what should be unit for it?

Since it has been transferred into a log scale, I could use dB, but the standard formula for dB calculation is something like 20log(x/x0), where x0 is the reference and cancels out the unit leaving us with just dB. I am not sure if I can the same thing in the above case since I don't have any reference power? Can I just assume something like x0=1 and report the numbers in dB?

If not, then I am not sure what should accompany dB. In one paper, I have seen for accelerations reported in gravities 1g = 9.8m/s^2,dB g^2/2 was used as the unit for the power spectrum? This makes me guess maybe I should use dB m/2s^2 or weirdly dB m^2/2s^4? Any ideas are greatly appreciated.


texts: Creation of Brownian noise as random walk

texts:
I am trying to create brown noise by integration of white noise (random walk).

I created Gaussian white noise. The histogram shows Gaussian distribution and Welch (and also FFT) shows flat spectrum.

Now I try to integrate it

$x(k+1) = x(k) + q(k)$

where $q$ stands for white noise and $x$ for resulting Brownian noise.

But the resulting spectrum is

According Wikipedia I would expect that there should be straight line in the plot.

What I am doing wrong?


texts: Confusion in proof of Inverse Laplace Transform

texts: For the proof of inverse Laplace transform, we change the integral from $\omega$ to $s$. I want to know the reason why we need to change the integral?


texts: FFT of a pure sine wave is not a single peak

texts:
For fun, I computed the FFT of a pure sine wave. I chose chose the sample length to be an even multiple of the signal period, so that I don't see windowing-effects. Here is the Matlab code that I used:

% Signal properties
f_signal = 1000;

% Sampling properties
f_sample = 4000;
T_total = 100;

% Calculate signal vector
t = 0:1/f_sample:T_total;
N = length(t);
V = sin(2*pi*f_signal*t);

% Calculate FFT
f_spectrum = (0:N-1)/N*f_sample;
V_spectrum = abs(fft(V))/(N/2);

% Plot
figure(2)
semilogy(f_spectrum, V_spectrum)
set(gca, 'XLim', [f_signal-5, f_signal+5])
xlabel('Frequency in Hz')
ylabel('Normalized amplitude')

But instead of a single peak, I get the following result:

This strikes me as odd for two reasons:

The effect is really huge. The center peak has an amplitude of 0.9 instead of 1, and the next bin to the right has an amplitude of 0.3 instead of 0.
While I would expect that limited numerical precision broadens my peak a bit, I am quite surprised that the effect is so silky smooth. I always thought of rounding errors as "noise" of the least significant digit. Noise usually leads to jarry spectra.

I observed that if I choose a slightly lower frequency (999.9975) then the effect is minimal. What is going on here?


texts: Why Zero Padding in the Center of the DFT Interpolates / Upsamples the Signal (Sinc Interpolation / DFT Interpolation / Periodic Interpolation)

texts:
I'm experimenting with the Inverse Discrete Fourier Transform. Starting from the two-cycles continuous $x(t)$ signal below:

I have the discrete signal $x(n) = \{ 1, 0, -1, 0, 1, 0, -1, 0 \}$ leading to the 8 points DFT $X_0(n) = \{ 0, 0, 4, 0, 0, 0, 4, 0 \}$

Now, if I use the IDFT on $X_0$, I obtain $x_0$ looking like that (the blue curve is the real part of the IDFT $Re[x_0(t)]$):

Here $x_0(n) = x(n)$ for integer values ($n = 0, 1, 2, \dots 7$). I understand I do not get back the continuous $x(t)$ function because of aliasing. I would explain that by saying on $x_0(t)$, there is a second signal, above the Nyquist frequency, and "riding" the "carrier"1

I read about zero padding, so I tried to add $k$ extra zeros in the middle of $X_0(t)$ and now, if I perform the IDFT, I obtain the following results:

1. With 8 extra values

I have $X_{k=8}(t) = \{ 0, 0, 4, 0, \underbrace{\mathbf{0, 0, 0, 0, 0, 0, 0, 0}}_\text{8 extra bins},  0, 0, 4, 0 \}$, leading to $x_{k=8}(t)$:

2. With 24 extra values

I have $X_{k=24}(t) = \{ 0, 0, 4, 0, \underbrace{\mathbf{0, 0, 0, 0, \dots, 0, 0, 0}}_\text{24 extra bins},  0, 0, 4, 0 \}$, leading to $x_{k=24}(t)$:

3. The problem

Adding more bins proportionally decreased the amplitude of the rebuild signal and increase the frequency of the signal riding the carrier. I think I understand both phenomenons.

However, I can't find an intuitive way of explaining why, at integer positions (the orange dots), $x_k(n)$ reproduces more and more accurately the original shape of the $x(t)$ signal.

1Do not hesitate to edit the question if I don't use the correct vocabulary here.


texts: Calculation of the clutter component of the space-time snapshot

texts:
I've encountered a silly problem during the implementation of SMI beamformer. I have KASSPER project datacube from IEEE dataport as my data and I want to calculate covariance matrix from that datacube. According to J.Guerci book "Space Time Adaptive Processing for Radar", The total space-time clutter return from a given iso-range can be calculated from formula:
$$ X_{c}= \sum_{i=1}^{N_c} \gamma_iv_i $$
where $ \gamma $ i is a complex scalar RV that accounts for the amplitude and phase
of the i-th clutter patch, and Nc is the total number of clutter patches in
the iso-range ring. For example for antenna with angles from -90 to 90 $ N_c = 180$, so there are 180 clutter paches and 180 complex amplitudes $ \gamma $ if I'm understanding right.
So, my question is how or where I can calculate/find those amplitudes $ \gamma $? In Ward book "Space Time Adaptive Processing for Airborne Radar" are formulas for estimation those amplitudes, but I want (if it's possible) extract those information from data, which is radar datacube. Is it possible with one datacube or I will need 180 datacubes for each one of 180 azimuths ? What I'm missing here?
Thank you for help.


texts: Why is the median filter called 'edge-preserving'?

texts: I want to know about median filtering. Recently, I read that it is best for preserving edges. But how does it preserve edges? I need complete knowledge of median filters. Why it is non-linear filter and why it is best, compared to linear filter?


texts: Optimum Filter Signal Detection for Non AWGN Channels

texts:
I have been reading this question and it confirms that the matched filter is the maximum-likelihood receiver in the presence of additive white Gaussian noise. So  in the AWGN channel it maximizes the signal-to-noise ratio at the decision instant at the end of each symbol.

However,according to Wikipedia

  AWGN channel is not a good model for most terrestrial links because of
  multipath, terrain blocking, interference, etc.

I would like to know:

If the transmission channel could not be considered as affected by AWGN. Which filter would yield the optimum SNR? I'm particularly  interested in the noise present in sound channel (sound channel through air not underwater).
Is it possible to "adapt" matched filter to non AWGN channels?


texts:
I have an image processing algorithm which uses Gaussian and Laplacian Pyramid Decomposition.
At each level of the pyramid the algorithm decimates the previous level by a factor of 2.
The LPF the algorithm uses is Gaussian.

I was wondering whether I could get the same results without the decimation step (By altering the LPF).

As far as I know, the decimation step is equivalent of cutting by half the bandwidth of the LPF.

What I tried, with no success, is to double the STD of the LPF at each step and skip the decimation.

Is there another / better way to do it?
Thanks.


texts: Error probability for orthogonal signaling with unequal symbol probabilities

texts:
I'm trying to derive the equation of error probability in orthogonal signaling. Here what I've tried:
Let $\vec{s_1} = (\sqrt{E_b} , 0)$ and $\vec{s_2} = (0, \sqrt{E_b} )$ be vector representation of $s_1$ and $s_2$. We know that $$f(\vec{y}|\vec{s_m}) = (\frac{1}{\sqrt{\pi N_0}})^2e^{-\frac{|\vec{y} - \vec{s_m}|^2}{N_0}} , \ m = 1,2$$ According to the Bayes' law, we can write $$P(\vec{s_m}|\vec{y}) = \frac{f(\vec{y}|\vec{s_m})P(\vec{s_m})}{f(\vec{y})}, m = 1,2$$Let $P(\vec{s_1}) = p_1$ and $P(\vec{s_1}) = p_2$. Using MAP rule, we want to find $m$ such that $P(\vec{s_m}|\vec{y})$ be maximized. So $\vec{s_1}$ is chosen if $$P(\vec{s_1}|\vec{y}) \gt P(\vec{s_2}|\vec{y})$$ and otherwise $\vec{s_2}$. After simplification we get $$\frac{N_0}{\sqrt{E_b}} \ln(\frac{p_1}{p_2}) \gt -2(y_1 - y_2)$$ Let $a = \frac{-N_0 \ln(\frac{p_1}{p_2})}{2\sqrt{E_b}}$ and note that $y_1 - y_2$ is a Gaussian random variable with mean $\sqrt{E_b}$ and variance $\sigma^2 = N_0$. So the probability of error is $$P_e = p_1\int_{-\infty}^a \frac{1}{\sqrt{2\pi N_0}}e^{-\frac{(z-\sqrt{E_b})^2}{2N_0}}dz + p_2\int_{a}^{+\infty} \frac{1}{\sqrt{2\pi N_0}}e^{-\frac{(z-\sqrt{E_b})^2}{2N_0}}dz$$ which simplifies to $$P_e = p_1\int_{-\infty}^{\frac{a - \sqrt{E_b}}{\sqrt{N_0}}} \frac{1}{\sqrt{2\pi}}e^{-\frac{u^2}{2}}du + p_2\int_{\frac{a - \sqrt{E_b}}{\sqrt{N_0}}}^{+\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{u^2}{2}}du = p_1Q(\frac{\sqrt{E_b} - a}{\sqrt{N_0}}) + p_2Q(\frac{a - \sqrt{E_b} }{\sqrt{N_0}})$$ Is this result correct? I couldn't find a source which includes error probability for the orthogonal signaling when prior probabilities are unequal.
Edit: Thanks to Dilip Sarwate's comment, the second exponent should be
$$-\frac{(z+\sqrt{E_b})^2}{2N_0}$$ So the final answer is $$P_e = p_1Q(\frac{\sqrt{E_b} - a}{\sqrt{N_0}}) + p_2Q(\frac{\sqrt{E_b} + a }{\sqrt{N_0}})$$ Where $a = \frac{-N_0 \ln(\frac{p_1}{p_2})}{2\sqrt{E_b}}$. We can simplify this formula and get $$P_e = p_1Q(\sqrt{\frac{E_b}{N_0}} - \frac{\sqrt{N_0}}{2\sqrt{E_b}}\ln(\frac{p_2}{p_1})) + p_2Q(\sqrt{\frac{E_b}{N_0}} +\frac{\sqrt{N_0}}{2\sqrt{E_b}}\ln(\frac{p_2}{p_1}))$$ which agrees with Dilip Sarwate's answer.


texts: How would I represent an image using a basis $A$, given a sparse or compressible image?

texts:
Looking for a Basis representation:
Given a sparse or compressible image M, how would one represent this image using a basis A?

If I let m = vec(M)

m = Ax

where A is my "representation" basis and x is the vector of coefficients in the A domain, what is the method for actually getting A?


texts:
I'm new into DSP and Filters in general. Its the first step ahead I do from learning and using a 1-pole IIR basic Low Pass Filter.

I'm about to implement for my VST audio application (using C++ and IPlug Framework) a RBJ's Peak filter.

By steps (from that document), the formula that I'll use within my Process() function will be y[n], thus (as code):

inline double Process(double input) {
    mBuffer += (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] - (a1/a0)*y[n-1] - (a2/a0)*y[n-2];
    return mBuffer;
}

Since I'll use peakingEQ, I need to calculate these coefficients:

b0 =   1 + alpha*A
b1 =  -2*cos(w0)
b2 =   1 - alpha*A
a0 =   1 + alpha/A
a1 =  -2*cos(w0)
a2 =   1 - alpha/A

which all depends by user defined parameters and intermediate variables.

My questions are:

What's the difference between using Q (thus sin(w0)/(2*Q)) or BW  (thus sin(w0)*sinh( ln(2)/2 * BW * w0/sin(w0))) on calculating alpha? Isn't Bandwidth = Q when reasoing as Peak filter?
What's the meaning of H(s) function (which is (s^2 + s*(A/Q) + 1) / (s^2 + s/(A*Q) + 1) for the peaking filter)? I don't see where to "insert" it in the formulas.

I hope you can help me to undestand these two dubts.


texts: Cut-off frequencies for fractional sample rate adjustment

texts:
We have a signal sampled att 22 kHz that we want to interpolate to 40 kHz.

So we can do this by upsampling by a factor 20, then downsampling by a factor 11. My question regards the choice of cut-off frequencies after the upsampling and the downsampling.

The old Nyquist frequency was at $11$ kHz so it would make sense to me to low-pass frequencies above $11$ kHz after upsampling - we're basically replicating our signal 20 times across the new spectra, right? The correct cut-off is apparently 20 kHz, but I can't quite wrap my head around that fact.

We then downsample by a factor 11 to 40 kHz. Now the correct LP cut-off is 11 kHz. I am seriously confused. Can anyone shed some wisdom here?


texts:
As I understand aliasing is because you sample a high frequency with a low sample rate, and the sample points don't have enough information about the high frequency component. Here is an example of aliasing:

w = 5.0;(*signal frequency*)
dt = 0.66125;(*sampling rate*)
f[x_] := Sin[w x]
ls = Table[f[x], {x, dt, 200 dt, dt}];(*sample points*)

ListPlot[ls, PlotRange -> All, Joined -> True,
 DataRange -> {dt, 200 dt}]

we can see in the plot that there is beating.
And the Fourier transform gives only one frequency(except the symmetric negative frequency one).

DFT[A_, ht_] :=
 RotateRight[
  ht/Sqrt[2 \[Pi]]*
   Fourier[RotateLeft[A, Length[A]/2 - 1],
    FourierParameters -> {1, 1}], Length[A]/2 - 1]

ListPlot[Abs[DFT[ls, dt]]^2, PlotRange -> All, Joined -> True,
 DataRange -> {-1/dt/2, 1/dt/2}, Axes -> False, Frame -> True]

question:

Why there is the beating but only one peak in the Fourier transform? Since there is only one frequency in the sample, what is it beating against?
How to calculate the beating frequency and the aliasing frequency?

Update 1

As chirlu points out there is not beating, if it is the beating because of adjacent frequencies, we would expect to see two separated peaks if we increase the Fourier transform resolution. Here I increase the frequency resolution by including more sample points, and we can see it is a single peak, not adjacent peaks.

ls = Table[f[x], {x, dt, 800 dt, dt}];(*sample points*)ListPlot[ls,
 PlotRange -> All, Joined -> True, DataRange -> {dt, 200 dt}]

Update 2

Here is an example using sound:

Play[Sin[3000 t], {t, 0, 2}, SampleRate -> 6000]

Play[Sin[3000 t], {t, 0, 2}, SampleRate -> 950]

So use a lower sample rate, we can actually hear the beating!


texts: Calculate a DFT of a simple finite signal {+1, −1, +1, . . .})

texts: Hey i was wondering whether approach presented in this video (6:10) is a good approach to solve this task. Cause it seems kinda tricky and i was wondering whether it can be done faster.


texts: Checking if a lecturer can be heard clearly from the back of a room

texts:
I would like to indicate to a human speaker standing at the front of a classroom if he/she is speaking loudly enough to be heard from the back. From some research, i have gathered that I am looking to calculate the Signal to Noise Ratio or something similar. It looks like this sort of thing can be calculated in MATLAB or preferably Python.

What types of inputs are needed? Is it possible to figure this out all from a mic in the back of the room? Or do I need a reference signal. For example, if i have a mic in the back of the room (connected via Bluetooth or WiFi) getting a noisy signal, will I also need a mic on (or near) the speaker so that i have something to compare to?

In the end, this would be running on a raspberry pi. Can anyone point me towards helpful libraries, modules, or tools for this type of thing?

This seemed like the best forum to ask this question. If there is a better place for it, let me know.


texts: Link between LTI systems and square-integrable functions eigenfunctions

texts:
I'm trying to study on my own a little bit of functional anaylsis.
I try to summarize the main concepts (without being rigorous) that I learnt. Correct me where and if I'm wrong.

In Hilbert spaces it is always possible to define a basis made of eigenfunctions and elements of the space can be rewritten using that basis. But how to find this basis?
A self-adjoint operator $L$ on a space $H$, with the eigenvalue problem $Lx_{n}=\lambda x_{n}$, defines a basis of orthonormal eigenfunctions for the space objects and for the same operator expansion too. Moreover, there are particular linear self-adjoint operators called Sturm-Liouville operators that, defining an eigenvalue problem, generate too a basis for the corresponding space I'm working with.

In particular, there are Sturm-Liouville "singular" problems that are defined over a space $H$ that is not "finite" but extends from $-\infty$ to $+\infty$: one of them, and for an engineer the most important, is $L^{2}(-\infty,+\infty)$. In this space, the operator $\frac{d^{2}}{dx^{2}}$ defines a Sturm Liouville singular problem and the solution of the eigenvalue problem is an infinite set of complex exponentials. From this theory we get the a continous spectrum for each signal of $L^{2}$ in terms of complex exponentials, i.e. a Fourier transform. I think that we can say that complex exponentials compose the spectrum of the operator $\frac{d^{2}}{dx^{2}}$ and of $L^{2}$ signals.

In engineering class, we met a particular class of operator called LTI (linear time-invariant) systems that have a peculiar characteristics: complex exponentials are eigenfunctions of these kind of operators.

My questions are:

1) Is there some kind of link between the operator $\frac{d^{2}}{dx^{2}}$ spectrum and LTI operators? (PS: I'm not looking for the demonstration on wikipedia, but a justification or a reference coming from spectral theory)

2) I know that LTI systems can be characterized with the convolution integral: maybe the convolution integral defines the same eigenfunction basis with the corresponding eigenvalue problem?

3) I also know that some LTI systems can be written as $$\sum^{N}_{k=0}  a_{k}\frac{d^k}{dt^k}y(t)=\sum^{M}_{k=0} b_{k} \frac{d^k}{dt^k}x(t)$$, where basically we have TWO operators acting on two functions $x(t)$ and $y(t)$: these two functions can be written in the same eigenfunction basis?

Sorry for the long post and thank you in advance if somebody can help me with these concepts or giving references! My usual references are Hanson-Yakovlev Operator Theory, Snider's PDE book and Oppenheim's Signals and Systems.


texts:
In a project, I have real-time voice audio being delivered in buffers of around 128 samples at a time. I want to calculate the autocovariance of the signal and use it to control an AR filter.

To do this, I use a Fourier transform and an Inverse Fourier transform.

Edit: Between the forwards and inverse FFT, I am basically just squaring (complex norm) and scaling the signal.

However, with such small buffer sizes I think the frequency resolution will be insufficient. An obvious choice is to fill up a larger, internal buffer and calculate a sliding DFT on the incoming data.

But I was wondering if there are specific algorithms for treating a whole block of incoming data, instead of doing it per individual sample?


texts: Recursive Least Square Adaptive Linear Equalizer

texts:
For the adaptive filter to work properly, a desired signal d(n) needs to be provided. The output from the equalizer y(n) is subtracted from d(n) to produce an error signal, which is used to adjust the filter weights.

The adaptive filter is located on the receiver side, so how to obtain the desired signal, and use it at the receiver?
When using a training sequence, the desired sequence can easily be located at the receiver. What will be the situation with other unknown sequences?
What is the difference between least mean square and recursive least squares adaptive linear filters?


texts: Unscented Kalman Filter - estimation of two parameters

texts:
I am working with Unscented Kalman Filter (UKF) in thermal modelling of a box.

I have 3 state variables ($T$ temperature, heating $h_h$ and cooling rates $h_c$) in my model and I am observing just one of them ($T$).
I am working with the filterpy python library following this example.
My main equation is the following:

$\hat{T} = T + (- h_c(T-T_e) + h_h ) \Delta t$

where $T_e$ is the external temperature - not measured - just assumed.

If I initialise the state matrix correctly (i.e. I provide an approximation of the simulated $h_h$ and $h_c$ parameters) the filter converges and provides a correct estimation of $h_h$ and $h_c$ coefficients.

If I initialise the state matrix providing just correct $h_c$ and I put a coarse guess of $h_h$, the filter estimates correctly the temperature but the cooling/heating coefficients converge to some nonsense numbers (from the physical point of view).

Is there any way of putting constraints on the state variables in the UKF?

Thx


texts: Symbol timing recovery: Polyphase vs piecewise linear interpolation

texts:
A symbol timing recovery scheme shown below has been successfully implemented in C++. Different TEDs (Mueller & Mueller, Early-Late, Maximum Likelihood, Gardner, Zero-Crossing, etc) are included in the implementation. I use a cubic interpolator I found in Michael Rice's book. The implementation performs sufficiently well, at least in simulations. I have seen that other implementations tend to use the PFB interpolator. Out of curiosity, does it bring any significant performance improvement over the cubic interpolator? The cubic interpolator is quite straightforward to implement and its performance seems to be quite modest.

The cubic interpolator from Michael Rice's book. $m(k)$ and $\mu (k)$ represent basepoint index (ideal sample) and fractional delay (between 0 and 1) respectively.


texts: Are there methods to measure how the change in one signal will affect another signal?

texts:
I have two distinct signals and I would like to compute their correlation. However, I think that just computing the crosscorrelation does not work here, since the signals are from different sources — also different amplitude, and length.

I would like to see that, for example, if the variance of the first goes up, the same happens for the second one. So, I plotted the variance, but I didn't see much. Are there any other characteristics I could look at? Energy?

By the way, I do it in Python.


texts: How do I show that two sets of spectra are different?

texts:
I essentially have 6 spectral curves, 3 repeats for two samples. I want to be able to show that two samples are distinct (if they are). I have tried spearman's correlation already, which does produce p values however I feel like I've very much hacked them and that the p value doesn't tell me if they curves are different like an anova test would say. One idea I had was doing PCA and anova on the first PC. Would this be valid?

Any suggestions welcome.


texts:
I am trying to prove that for all $x,y \in l^{2}(\mathbb{R})$, it holds that $$\left\langle\left(\downarrow M\right)\left[x\right],y \right\rangle = \left\langle x,\left(\uparrow M\right)\left[y\right]  \right\rangle \text{.}$$
How should I proceed to prove this?


texts:
I am trying to investigate the use of NIR cameras in eye-tracking research e.g. the use of eye-tracking with infra-red glint (Purkinje Image). What I would like to know is what Make and Model of camera you have used.

I understand the choice of camera will be dictated by the nature of your research. This is not important for now. I am particularly interested in off-the-shelf cameras rather than custom built hardware.

I would like to know the make and model so that I can investigate how best to support different hardware for an app that I am creating. I am trying as much as possible to make the app open-hardware.

The app is focused on the research community, hence I would to get an idea of the cameras you have used in your eye-tracking research.

I will then investigate the camera specs online to determine things like what sort of programming interface the camera supports, it's features such as frame rate, hardware region of interest etc.


texts:
I'm a little confused about the sample rate:

If I have an audio file that was recorded with a sample rate of 48kHz, this means that there are 48,000 values per second in this file, right?

So when I read this file with a smaller sample rate of lets say 44.1 kHz I am summarizing every 44,100 values to one second so effectively I make the audio slower and thus longer than it actually is. Am I getting this correctly or do I have a total misunderstanding?


texts: Does it make sense to use a wavelet that is equal to a sine of one period?

texts: Would I get same results with STFT method (short-time Fourier transform) if I try this?


texts: Iterative Kalman filters and system parameters estimation

texts:
i am working recently on a project in which i want to implement a Kalman filter as being an observer, and i couple this observer with a state feedback controller that produces control actions necessary to control a mobile robot.

The problem is that the governing equations for the mobile robot's motion contain errors, and this will affect the control actions produced, and at the same time the miss modelling error covariance matrix is hard to get as well as the measurement error covariance matrix. I though about using some arbitrary values for these matrices and go on with the filtration process, but i thought that it would be much better if there is anything like a kalman filter formulation that has the ability to estimate closer values for system parameters so better control actions can be obtained.

Thank you in advance.


texts: What's the difference between "Magnitude Spectrum" and "FFT"?

texts: I am particularly talking about the matplotlib's axes.magnitude_spectrum() function that plots the magnitude spectrum. Is this a representation of the discrete Fourier transform?


texts: Advice on filtering a CW signal for processing

texts:
This question is related to my other recent questions.  In the last question, I learned that I should filter about 4 * WPM (words per minute) around the carrier frequency if I want to recover Morse code.  One problem is, learning the words per minute.  I filter because I want to differentiate times when the system is keyed, to times that are just noise.  Filtering requires WPM, WPM seems to require filtering.

Now I am simply going through the whole audio, and if a window has an element above a certain threshold, the whole window is marked as being 'on' or 'keyed'.  The window is big enough for a whole alternation of the carrier to be examined (sampling frequency over carrier freq, or fs/fc).

Are there better ways to approach this?  I have been working on this for awhile now.  I am doing this in Python.


texts: How do you properly organize data to compute multiple (independent) recursive filters at the same time taking advantage of SIMD instructions?

texts:
I'm processing multiple (independent) Exponential Moving Average 1-Pole filters on different parameters I have within my Audio application, with the intent of smooth each param value at audio rate:

for (int i = 0; i < mParams.GetSize(); i++) {
    mParams.Get(i)->SmoothBlock(blockSize);
}

...

inline void SmoothBlock(int blockSize) {
    double inputA0 = mValue * a0;

    for (int sampleIndex = 0; sampleIndex < blockSize; sampleIndex++) {
        mSmoothedValues[sampleIndex] = z1 = inputA0 + z1 * b1;
    }
}

I'd like to take advantage of CPU SIMD instructions, processing them in parallel, but I'm not really sure how I can achieve this.

Any tips?


texts: Why does the frequency window affect the inverse fourier transform oscillation frequency?

texts:
I have an oscillating pulse in the frequency domain that I would like to inverse Fourier transform. My signal looks like:

Which was coded in MATLAB using the following code:

sampleRate = 1000;
freq = 500 : (1/sampleRate) : 1500;
intensityFreq = exp(-((freq-1000).^2));
signalFreq = sqrt(intensityFreq).*exp(-1i*10*(freq-1000));
plot(freq,signalFreq)

When I inverse Fourier transform, it shows a similar oscillating pulse (in time) whose oscillations depends on the frequency range I used in the above code. For example, when:

freq = 500 : (1/sampleRate) : 1500;

is used (as in the code shown above), the inverse fourier transform looks like:

But when:

freq = 980 : (1/sampleRate) : 1020;

Is used instead, the inverse Fourier transform looks like:

This makes no sense to me, since the frequency range shouldn't matter since signalFreq goes to 0 around these points anyway. Why does the altered frequency range affect the oscillations in the inverse Fourier transform so much? How do I obtain the "actual" time domain signal that doesn't depend on frequency windows?

Any help is much appreciated. If this post is not clear I would be happy to provide more detail.

EDIT: The full code used to find the inverse Fourier transform is:

sampleRate = 1000;
freq = 500 : (1/sampleRate) : 1500;
intensityFreq = exp(-((freq-1000).^2));
signalFreq = sqrt(intensityFreq).*exp(-1i*10*(freq-1000));
Y = ifft(signalFreq);
plot(real(Y))


texts: Bandwidth expansion by computing missing harmonics in music?

texts:
In old recordings or highly compressed audio files (lossy), the higher frequencies are lost, and we hear a distinctive muffling effect. (Just like if the song was played through a wall)

Is it possible to compute the missing harmonics using any method (fourier, bilinear, etc.) and add those back to the sound/song? What is currently feasible?

*It is true that those frequencies are lost, and I am not proposing to boost the higher frequencies (as it would raise the noise floor), but I thought there would be some way to approximate the missing frequencies, since we know what all the instruments should sound like (we know their corresponding harmonics)...*

I'm well aware of the fact "garbage in, garbage out", but even so, our brain itself have very advanced processing techniques. (3D visualization using only 2D sources in vision and audition) Wound't it be possible to fool the brain in "hearing" a better quality sound, and give old recordings back their original sharpness?

Edit: To be more clear, I'm proposing to make audio "sound" better, not to recover lost information.


texts:
I'm using MATLAB, and have experience with it but not with this kind of problem. I'm trying to identify white lines in noisy 2D images such as the one below. I can separate the brighter ones just by thresholding, but the dimmer ones (like in the bottom left) are tricky. Any suggestions? I guess I need some kind of ridge detection. It's worth noting that the background isn't entirely uniform... the brightness of the dim filaments might even be lower than the background noise in other parts of the image. So maybe some kind of local filtering instead of uniform background subtraction. As you can see I don't know where to begin.

Ideally I would will want to produce a skeleton describing all the lines for further analysis (measuring lengths). Of course this involves dealing with cross-overs and stuff too, but for now if I can just detect all the lines that would be a good start.


texts: Why idft(dft(a) * dft(b)) not equal to convolve(a, b)?

texts:
I'm a little confused... I always thought the DFT of a convolution was equal to a product of DFTs, but when I tried this in Python:

from scipy import *

a = [1+0j, 2+0j]
b = [4+0j, 5+0j]

print list(ifft(fft(a) * fft(b)))
print list(convolve(a, b))

I got back:

[(14+0j), (13+0j)]
[(4+0j), (13+0j), (10+0j)]

Why are they not the same thing?


texts: Matlab signal processing removing white noise from signal

texts:
What is the way to address the removal of white noise from a signal?

Extra Details
Hi I also agree the details are as follow:                                                                                The signal is a acoustic or a sound wave. The the power associated with the frequency of the noise is distributed evenly across the whole frequency range (white noise) The actual signal characteristic is not known as its a inverse problem. But the actual frequency of the signal is a part of the total frequency. Power associated with the signal is higher than the noise can be seen in a spectogram.
Actuallty I also cant figure out how to extract the coefficients of the power/frequency  (db/Hz) from the spectogram. Can you suggest a way?
 see the graphs in the link

http://www.mathworks.com/help/signal/ref/spectrogram.html?refresh=true

Actual problem is to find the time of arrival of the actual signal


texts: Employing MIMO under the two-ray propagation model

texts:
I'm a little confused about how 'rich' the scattering environment should be for MIMO to work (i.e. for the channel matrix to be full rank). For instance, let's say we have a Tx and Rx with two antennas each, and under the two-ray model (or the flat/plane earth model), one path is direct from the Tx to the Rx and another path follows a reflection from the ground.

Can MIMO be used here? In other words, can two separate spatial streams be obtained?


texts: Maximum likelihood estimation complexity computation

texts:
I have a basic question about maximum likelihood (ML) estimator and its implementation.

I am trying to simulate a communication system, while using ML at the receiver side to find the transmit sequence. The received signal is given by:

$y=Hx+n$

$y, H, x$ and $n$ indicates the received signal, MIMO channel of $N_t$x$N_r$ transmit and receive antennas, transmitted signal and AWGN, respectively.

I'm trying to detecto $x$ based on ML detector, so the estimated signal $x$ can be written :

where $h_{ln}$ represent the $l^{th}$ column of the channel $H$.  I have read in some books and papers, such as in This paper  that the complexity computation of the above equation is $6N_rN_tM$, where M represent the constellation order of signal, where it was mentioned that $h_{ln}x$ requires $4N_r$ for $Hx$ and $2N_r$ for $|y_n-h_{ln}x|^2$.

My question, in $4N_r$ and $2N_r$ what do the $4$ and $2$ represent?? In other words, the comlexity computation doesn't depend on $x$ too? for example, if we up-sampled the signal $x$ in transmitter or we spread it using Walsh code, the complexity computation won't change as long as we are using the same number of $N_t$ and $N_r$ ??

thank you


texts: "removing" parts of an image with matlab

texts:
for an assignment, I'm supposed to work with matlab, with the image below.
The task is to convert this image to a grayscale version (which wasn't a problem) and then to alter the image with imadjust so that only the writing and the arrows (not the shells) are visible. Since both the shells and arrows are yellow, I'm wondering how that would work. I tried different low high respectively in- and outputs for imadjust and also different gamma values, but none of that works.

Thanks!

edit: here's the grayscale version I got. I used imadjust (img,[0.9;1],[],0.15)


texts: Detect human speech in realtime audio on mobile phones

texts:
I am looking to develop an Android app. As part of the functionality the app would require to randomly sample 3-5 seconds of audio and classify it as containing human speech or not. I understand that this concept is called Voice Activity Detection?

What would be the best way to implement this on a Mobile phone. I developed a basic system using energy based features and thresholds. I am hoping to find something less susceptible to noise, probably using features such as MFCC or formants? I did go through a number of papers, but most of them would require me to collect data and train models. Is there any library or framework I could use which would work in realtime?


texts:
So I just started reading on FIR and IIR filters and am amazed by how "simple" the theory seems, so far.

But what confuses me is, why does filtering work by creating a weighted sum of previous samples?
What intuition makes one think this can produce desired filtering effects?
It seems a bit unintuitive to me eventhough anyone can verify that summing delayed signals together does produce comb filtering. But desirable filtering? Why?


texts:
I'm trying to create a basic GNURadio flow graph that will transmit a short message in order to receive an appropriate reply back.

I currently have the following flowgraph:

I have a program that sends UDP messages to the UDP Source and receives the UDP messages from the UDP Sink

I am able to receive a message from the UDP Sink only if I loop the message I send from the program a couple of times (about 250 at least). Because of this I receive multiple repeated messages, and I wish to receive just one. One way of solving this is through programming, but it is proving to be quite challenging (but should be possible).
Is there a proper way through altering the flow graph parameters with which we can transmit a message without looping or repeating?


texts:
I have this matrix

The following points

R1 = (1,7) = 0
R2 = (3,2) = 6
R3 = (5,5) = 3
R4 = (7,3) = 1

And using this equation |g(seedi) – g(pixel)| < T
for T=1

And this is the solution for threshold = 1

I don't understand why all the ones for example, belong to R1 = α

Based on the equation we have |0 - 1| < 1, which means it shouldn't become α.


texts: Matlab - How to create an autocorrelogram using a spike-train

texts:
I have already posted this question on MATLAB's question forum.  However, I haven't received any responses, and am hoping to find some folks who may be a bit more signal processing savvy on here.  Here goes... I'm wondering if someone could help clear up some confusion I am having surrounding the xcorr() function since I am struggling to parse how it works exactly.

Here is the matlab description:

  r = xcorr(x) returns the autocorrelation sequence of x.

I am trying to make an autocorrelogram out of spike train data. Thus, I am inputting timestamps into xcorr() like so... xcorr(spiketimes). I am operating under the assumption that autocorrelating time-stamps for events (spikes) should lead to an autocorrelogram representative of the the autocorrelation between interspike intervals. This autocorrelogram, I believe, should provide the same info as an autocorrelogram of 0's and 1's (0 = no spike, 1 = spike) inputted into a vector who's length is the size of the recording (in milliseconds).  When I run that line of code, however, I find that there's a peak towards the length of the timestamps vector. This doesn't make sense. The peak should be at lag zero (where the vector is 100% correlated with itself). What am I doing wrong here? Or is my understanding of an autocorrelogram incorrect? Does anyone here have any experience using xcorr(), and have any advice?

Also, this is less important and a bit technical, but for fellow matlab users, I'm not totally sure how to adjust bin width. If anyone has some pointers, I'd appreciate it. Finally, the fact that I can't center the plot that's created from xcorr() at zero is driving me crazy. How can I fix this?


texts: Don't understand coefficients in MATLAB "butter" low pass result

texts:
I don't unstand the nature of the output in Matlab's "butter" command. Calling

[b,a] = butter(4, 0.5, 'low');

gives me two vectors, which the documentation says are b for numerator, and a for denominator, of the transfer function.

However plotting the frequency response of this transfer function with freqz() does not give me what I expected.

If I plot

freqz(b./a)

I see:

This is not a low pass filter, this is some sort of band stop filter suppressing the middle frequencies. However if I plot

freqz(b)

I see what looks like the low pass filter I was hoping for:

So if I want to use a Butterworth filter to low pass my data, do I only use b? And if so what is the point of a?


texts: Understanding the relationship between the convolutional feature map and the region proposals in Fast R-CNN

texts:
In the paper Fast R-CNN available here, I am trying to understand the relationship between the region proposals and the convolutional feature map.

So from what I understand, Fast R-CNN takes in an image and a set of region proposals for object detection, with the goal of taking each proposal and giving it a classification from a set of classes, and a bounding box (or none in the case of background regions). More processing can be done after with the region proposals to merge them etc. but in general this is how it works.

What I don't understand is how the region proposals relate to the convolutional feature map. In section 2.2 the authors state that they test their architectures by pre-training networks with some image net architectures, which we will modify by replacing the max pooling layer by a RoI pooling layer, followed by some fully connected sibling layers for regression and classification. My question is what is the relationship between a region proposal, say a region $(x,y, h ,w)$ given by the top left corner and a height and width, and the convolutional feature map, which may have spatial dimensions (and definitely more channels) that differ from the original image, (see the image below)? It would make sense certainly if the feature map had the same spatial dimensions as the original image, but otherwise I'm not sure.


texts: Dynamic range compression in the frequency domain?

texts:
I'm trying to do frequency-band specific dynamic range compression (like multi-band compression) and was thinking of whether dynamic range compression is doable in the frequency domain or whether it's done in the "typical" filtering, gain reduction and then mixing (recombining filtered and GR'ed signals) way using time-domain signals?

Also, what are the benefits/drawbacks over time-domain?


texts: Having Glitches trying to use a Sinc Function Interpolation on a Signal, by chunks of 1024 samples

texts:
I am trying to interpolate a signal using a third party lib that performs Sinc Function Interpolation on a signal. It works well if I input all the samples at once but introduce phase distortion if I try to perform by block of 1024 samples.

Some context:
I am trying to increase or reduce the playback speed of an audio loop.
The audio loop is a stereo interleaved samples with a total of 352,800 samples (176,400 stero frames).

It works flawlessly when I input all the 352,800 samples at once, but this can't be performed real time.
Instead, I try to feed it with chunks of 1024 samples as requested by the audio interface.

To do that, I first calculate how much I have to feed the interpolation function to get an interpolated result that fits this requirement:

// how much samples we need to prune in the original samples
let to_take = asked_size * playback_rate;

Then I pass theses samples to the interp function that spits out excaclty 'asked_size' samples, theses samples have been interpolated.

The interpolation function takes the following parametters:

samples: the samples buffer to be interpolated (1024*playback_rate)
sinc function (this is the sinc function as in the lib).
playback_rate
delay (which I set to 0.0)
wing (how many samples to looks from each sides for interpolation / convolution).

When the playback rate is 1, everything is smooth, but when I move this value, the global signal get distorded, because there is gaps between the first sample of the next buffer and the last sample of the previous buffer that I send to audio.

I suspect there is something to do with the phase, or something to do with the fact that the interp function miss some contexts on previous chunks of samples (blocks of 1024) for the outpout (concatenation of all the blocks of 1024) to be seemsless, and this results in nasty distorted audio.

How should I correct this behaviour ?

EDIT: Added implementation in Rust:

// block loop
pub fn next_block(&mut self, size: usize) -> Vec<Stereo<f32>> {

  // how much we need to prune
  let to_take = ((size as f64) * self.playback_rate) as usize;

  // to interleaved samples (original frame format is Vec<[f32, 2]>)
  let interleaved: Vec<f32> = self.take(to_take)
     .flat_map(|x| x.to_vec()).collect();

  // to complex DspVec data struct
  let mut complex = interleaved.to_complex_time_vec();

  // interpolation is HERE !
  complex.interpolatef(
    &mut self.interp_buffer,
    &self.sinc_function,
    (1.0 / self.playback_rate) as f32,
    0.0,
    16,
  );

  // re-frame the signal
  let chunked = complex.to_slice().chunks(2);
  let mut out = Vec::new();
  for chunk in chunked {
    let c = [chunk[0], chunk[1]];
    out.push(c);
  }

  // send full buffer
  return out;

}


texts:
I am trying to understand the bode plot result I have obtained.

The blue, green, and red plots are $a$, $ b$, and $a+b$ respectively.

What I'm having a hard time understanding is how a+b is formed.  Looking at the magnitude of $a$, it is in the range of $500 dB$ or $$10^{25}$$.  Whereas, looking at $b$, it is a negative dB gain, so a number of small magnitude.  Although I do expect to see the gain of $a+b$ to be lower than that of $a$, I can't explain why the drop is so significant.

Am I missing something fundamentally important here? or is this fishy to anyone else?

If it makes any difference, I directly got this plots in matlab using

bode(a,b,a+b)

TIA


texts: Correlation of two signals at different sample rate

texts:
I have a signal $X$ and I want to find signal $Y$ in $X$, I can achieve this by cross-correlation. Signal $X$ is sampled at $8000\textrm{ Hz}$, and $Y$ signal is sampled at $44100\textrm{ Hz}$.

My question is "should both signals be sampled at same rate to correlate to find $Y$ in $X$ " Or is it okay to correlate two signals sampled at different rate to find  $Y$ in $X$?


texts: Can channel capcity be explained in simple terms?

texts: Can channel capacity be explained in simple terms - without mutual information and such probabilistic concepts? In its most general form, what are all the parameters that it depends on?  What is its dependence on physical parameters such as the electrical parameters if we talk of information conveying through electrical signals?  Can we measure/calculate the channel capacity of any channel given its physical parameters?


texts:
Given plots of magnitude and phase of a Fourier transform, how can I obtain the analytical formula for the signal?

I know that $$X_1(j\omega) = 3j\omega[H(\omega+3\pi)-H(\omega-3\pi)]$$ but I can't figure out a way from this point to the answer, which is $$x_1(t)=\frac{3}{\pi t^2}[3\pi t\cos(3\pi t)-\sin(3\pi t)]$$