texts: How to understand energy function in thin plate splines?

texts:
guys! I have carried out some experiments using matlab warp tools from here. But the results are not good, so I decide to read its original paper Principal Warps: Thin-Plate Splines and the Decomposition of Deformations from here to find reasons. However, I am stuck to part II.....

It says the energy function $\mathbf{I_f}$  is proportional to $\mathbf{WKW^T = V(L^{-1}_nKL^{-1}_n)V^T}$, then, I can not understand the whole first paragraph of part II:

Why does the pattern of y-displacements represent an expansion of lengths along y-axis? Actually, I know energy function, but I can not understand its relations to $\mathbf{WKW^T}$, then the pattern of displacements, and then the actions around x or y axis. So, I ask this question about energy function. Thanks a lot!


texts: Does there exist a simple regression method to fit a single sinusoid period to a data set?

texts:
I'm currently working on a project where I want to fit a single sinusoid period to a data set. Essentially I have very good control/knowledge of the signal's dominant frequency, so I'm only sampling for a single period. What this leaves me with is the problem of determining the amplitude and phase offset in the presence of the sensor's noise.

I'm curious if anyone is aware of a easy to implement method for estimating these two parameters that is optimal in some sense. Most of my research results have been complex algorithms implemented in MATLAB which is not suitable for my application (I'm trying to measure this in real-time and don't have a license).


texts:
I have been asked  to study an experiment which says "what is the effect of variation of duty cycle of sampling frequency on  a reconstructed signal"

I am not an electronics/electrical engineer I am CS guy.This experiment probably is a part of Analog to Digital Converion techniques, by now I have been able to understand Nyquist theorem and reconstruction is based on this.But for above experiment I am not clear with what to do or study so far what ever I have googled has not given me fruitful results let me know if some one has an answer to above problem.Following is the link which I have studied and could follow http://mhakhan.tripod.com/07sampling.pdf  a bit for the above question,what is exactly the thing which I should search?


texts: Filter design - Rules of thumb regarding ratio of sampling and stop/pass frequency

texts:
The issue: I have a signal in the time domain sampled at 51.2KHz. I want to apply a highpass filter with a pass frequency of 20 Hz.

The problem: Using Matlab to do that, the design of the filter takes forever to design. Having made a bit of research on this, it seems that this happens because of the high ratio of $\frac{sampling frequency}{Fpass}$ .

This got me wondering, is there a rule of thumb regarding this quantity?


texts:
Under the simulation of a communication system suffering from noise (AWGN) and co-channel interference, what is the appropriate way of calculating SINR?
Could one calculate it at a packet level as:

$$
SINR_{packet} = \frac{\sum{\left \| x \right \|^{2}}}{\sum \left \| n \right \|^{2}+\sum \left \| z \right \|^{2}}
$$

where $x$ is the total useful waveform for that packet, $n$ is the corresponding noise waveform and $z$ is the waveform of the interferer, or is it necessary to calculate it on a per bit basis as

$$
SINR_{bit} = \frac{\sum{\left \| x_{b} \right \|^{2}}}{\sum \left \| n_{b} \right \|^{2}+\sum \left \| z_{b} \right \|^{2}}
$$

where $x_{b}$ is the useful waveform for that particular bit, $n_{b}$ is the corresponding noise waveform for that bit and $z_{b}$ is the waveform of the interferer for that bit and then take the average across all bits of the same packet like so:

$$
SINR_{packet} = \frac{1}{k}\sum_{bit=1}^{bit=k}{SINR_{bit} = \frac{1}{k}\sum_{bit=1}^{bit=k}{\frac{\sum{\left \| x_{bit} \right \|^{2}}}{\sum \left \| n_{bit} \right \|^{2}+\sum \left \| z_{bit} \right \|^{2}}}}
$$

Which one is supposed to converge more quickly and give more realistic results? Which one would be better related to bit error rate (BER) and which to packet error rate (PER)? Finally could Central Limit Theorem be of any use proving that the first and last expression converge to the same number for large number of bits in a packet?

I believe real systems calculate SNR as in the first formula, except of course for the $z$ term.


texts: Output power of a signal passing through an IIR filter

texts:
suppose I have the following discrete time signal
$x\left(n\right)=0.7\sin\left(2\pi \times 623\:nT\right)$
where $N=4000$ and $T=\frac{1}{f_s}$, where $f_s=40$ kHz
I pass this signal through a simple IIR filter $ H\left(z\right)\:=\:\frac{0.2}{1-0.3z^{-1}}$
I calculate the power of the ouput signal $P_y\:=\:P_x\:\left|H\left(z\right)\right|^2_2$
basically am using the filter 2-norm in this case.
for that my calculations would be as follows:
$P_y\:=\left(0.245\right)\left(\frac{0.2^2}{1-0.3^2}\right)\:=\:0.0108$
However when I develop this in Matlab, I get the output power as 0.0199
Here's my code in Matlab
fs = 40000;
f1 = 623;
N = 4000;
n = 1:N;
b = 0.2;
a = [1 -0.3];
x = 0.7*sin(2*pi*f1*n/fs);
P_x = var(x);
y = filter(b,a,x);
P_y = var(y);

I don't understand what I was doing wrong, would someone please clarify for me? I would appreciate it. Thank you in advance.


texts: Directly compare subpixel shifts between two spectra — and get believable errors

texts:
I have two spectra of the same astronomical object. The essential question is this: How can I calculate the relative shift between these spectra and get an accurate error on that shift?

Some more details if you are still with me. Each spectrum will be an array with an x value (wavelength), y value (flux), and error. The wavelength shift is going to be sub-pixel. Assume that the pixels are regularly spaced and that there is only going to be a single wavelength shift applied to the entire spectrum. So the end answer will be something like: 0.35 +/- 0.25 pixels.

The two spectra are going to be a lot of featureless continuum punctuated by some rather complicated absorption features (dips) that do not model easily (and are not periodic). I'd like to find a method that directly compares the two spectra.

Everyone's first instinct is to do a cross-correlation, but with subpixel shifts, you're going to have to interpolate between the spectra (by smoothing first?) -- also, errors seem nasty to get right.

My current approach is to smooth the data by convolving with a gaussian kernel, then to spline the smoothed result, and compare the two splined spectra -- but I don't trust it (especially the errors).

Does anyone know of a way to do this properly?

Here is a short python program that will produce two toy spectra that are shifted by 0.4 pixels (written out in toy1.ascii and toy2.ascii) that you can play with. Even though this toy model uses a simple gaussian feature, assume that the actual data cannot be fit with a simple model.

import numpy as np
import random as ra
import scipy.signal as ss
arraysize = 1000
fluxlevel = 100.0
noise = 2.0
signal_std = 15.0
signal_depth = 40.0
gaussian = lambda x: np.exp(-(mu-x)**2/ (2 * signal_std))
mu = 500.1
np.savetxt('toy1.ascii', zip(np.arange(arraysize), np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth), np.ones(arraysize) * noise))
mu = 500.5
np.savetxt('toy2.ascii', zip(np.arange(arraysize), np.array([ra.normalvariate(fluxlevel, noise) for x in range(arraysize)] - gaussian(np.arange(arraysize)) * signal_depth), np.ones(arraysize) * noise))


texts:
Gentlemen! I have to reproduce MATLAB random normal generator (randn) with some fixed seed.

Also, I have to reproduce it in an external script, without MATLAB itself. Has anyone done this? What random normal generator implementation is better to use?

P.S. Preferred languages are Python and C/C++, but other languages may be useful too.


texts: bandwidth limited pulse,what window function is it?

texts:
I want to know what kind of shape a multiple cycle sine wave pulse ( or maybe the correct term is burst ) should have to have the narrowest bandwidth possible at given pulse width. In the laser science field, modelocked lasers producing ultrashort pulses are said to make transform limited pulse. Transform limited and bandwidth limited are the same thing: it's a pulse with minimum spectral width possible at the specific length of that pulse. When I was looking at the amplitude envelope shape of these pulses, it reminded me of the window function like it's used for example in fft, and it looked kind of like gaussian (the sort of gaussian that starts and ends with zero):

My question is what amplitude envelope shape should that sinewave multiple cycle pulse have to have minimum spectral bandwidth for its duration? At first I thought, that's easy, it must be gaussian (the "confined gaussian" in wikipedia, zero at edges), but when I thought about it more I am not so sure anymore. Wouldn't for example a pulse with Nuttall shape be narrower spectraly? What about Hanning and Blackman? I know all these look almost the same but hey, the true transform/bandwidth limited pulse can be only one shape, so which is it? https://en.m.wikipedia.org/wiki/Bandwidth-limited_pulse


texts: Compute real signal from a discrete analytic signal

texts:
I have a 128MHz-wide down converted signal that gets processed in a FPGA via a polyphase filter bank to give 8x16 MHz baseband analytic signals. How do I convert this analytic signal to a real-valued signal? The Wikipedia article on analytic signal mentions that interpolation/up sampling will be required if the original signal was down converted. Here is how I plan to proceed:

Extract the real part of the analytic signal
interpolate (up-sample by a factor of 2, followed by a low pass filter)

Would the above give me a real-valued signal?

EDIT 1: I believe, I have to up-sample (i.e add zeros to every alternate point) and then use a low pass filter, e.g. a FIR to interpolate.


texts: Linear convolution of discrete signals with defined lengths

texts:
What is the maximum count of non-zero elements, that can a linear convolution of discrete signals of "lengths" 5 and 7 have?

When I label the length of signal $x[n]$ as $M$, and the length of
signal $h[n]$ as $N$, then the length of their convolution is the signal $y[n]=M+N-1$. However this is not that thing I was looking for. So how can I find the maximum count of non-zero elements?


texts: 3dB and 6dB Per Octave Rolloff Digital Filters / Generating Pink and Brown Noise

texts:
So I would like to have a design method for filters with 3dB and 6dB per octave roll off -- for generating pink and brown noise respectively.

I know the following 'pinking' filters exist:
Filter to add 3dB per octave?

But the poster says he doesn't remember the value for A, nor does he give any idea of how he did it in the first place.

So my main question is:

How do we pick the poles and zeros of such a filter?


texts:
What is the purpose of DFT & IDFT blocks which are placed after the Analysis filter bank & before the Synthesis filter bank.
The reason i am asking this question is because I tried to design Polyphase DFT filter bank & got the same output with & without DFT-IDFT block.


texts: How to do image reconstruction from array data?

texts: i have 15 set of data which is the s-parameter, consist of 2 column (1001x2) each set. The first column is frequency and the second column is magnitude in dB. I have transform the magnitude from (dB) into value and frequency into time. But i have no idea how to transform these 15 sets of data into an image by using Matlab. Can anybody help me?


texts: Deconvolve using divsion of FFT of shifted signal in time domain

texts:
I have FFT of two signals. Y=120 , Y1=80 Hz. to convolve in time domain I can convolve them using their ffts. as Y3= (Y).*(Y1);
I can get Y by dividing as Y3./(Y1)

Y=Y3./(Y1);

But if take ifft of Y3. and shift it to 5 sample and at the end of the signal is padded with some zeros. after shifting I take FFT(shifted_Y3). and try to get Y using FFT(shifted_Y3)./(Y1). But I get distorted signal... Why that so? although frequency  spectrum is same of Y3 and shifted_Y3.
last left plot is result of fft(shifted_Y3)./Y1 and last right plot is result of Y3./(Y1)

Why I am getting recovered signal with added components. Is there any way to avoid this in shifted signal?

Fs = 1000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 1000;             % Length of signal
t = (0:L-1)*T;        % Time vector
S = 0.7*sin(2*pi*120*t);

bb=zeros(1,200);

figure
%subplot(4,2,1);
%plot(S)
Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
subplot(4,2,1);
plot(f,P1)
title('Single-Sided Amplitude Spectrum of Y(t)')
xlabel('f (Hz)')

S = 0.7*sin(2*pi*80*t);
%subplot(4,2,3);
%plot(S)
Y1 = fft(S);
P2 = abs(Y1/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
subplot(4,2,2);
plot(f,P1)
title('Single-Sided Amplitude Spectrum of Y1(t)')
xlabel('f (Hz)')

Y3=Y1.*Y;
P2 = abs(Y3/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
subplot(4,2,3);
plot(f,P1)
title('conv by multiplying fft of both above signals Y3 =Y1.*Y;')
xlabel('f (Hz)')

out_ifft = ifft(Y3);
subplot(4,2,4);
plot(out_ifft)
title('conv-result in time domain')

 div_fft = Y3./Y1;
  P2 = abs(div_fft/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
  subplot(4,2,8);
plot(f,P1)
title('in fft not-timeshifted./Y1')

c=cat(2,out_ifft,bb);
out_ifft1=c(5:1:1004);
subplot(4,2,5);
plot(out_ifft1)
title('in-timeshifted= 5 sample shifted conv-result in time domain')

 in_fft = fft(out_ifft1);
 P2 = abs(in_fft/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
 subplot(4,2,6);
plot(f,P1)
title('in-fft-timeshifted = Single-Sided Amplitude Spectrum of 5 sample shifted conv-result')
xlabel('f (Hz)')

  div_fft = in_fft./Y1;
  P2 = abs(div_fft/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
  subplot(4,2,7);
plot(f,P1)
title('in-fft-timeshifted./Y1')


texts: Dimension of channel matrix elements in MIMO systems

texts:
In a MIMO system we have:
$y= \mathbf{H} \star \mathbf{x} +n$ with $x \in \mathbb{C}^{N_t \times 1}$ ($N_t$ the number of transmission antennas) and $\star$ the convolution product.

For the matrix, $\mathbf{H}$ we have: $\mathbf {H}=(h_{i,j}(t,\tau))$ with $i \in \{1,...,N_r \}$, $j \in \{1,...,N_t \}$, $t$ the time and $\tau$ the delay for the impulse response

Here is my question:
What is the dimension of $h_{i,j}(t,\tau)$?It's a number or a vector because in many book, there is $size(\mathbf{H})=N_r \times N_t$. But if $h_{i,j} (t,\tau)$ is a vector, it's false.

An example is article Survey of channel and radio propagation models for wireless MIMO systems.


texts:
The values that one can find in a wave file, e.g.

    0.0036
    0.0026
    0.0174
    0.0050
    0.0026
    0.0108
    0.0154
   -0.0114
   -0.0006
    0.0115

are - as I assume - relative values between -1 and 1, the range you normally see when watching a waveform in some audio editing program. Now, how are these relative values mapped when you playback a file? E.g. if there is a value of 0.5, does it mean playback at half of the maximal possible loudness? Is there some kind of reference value stored in the WAVE file to which the relative values refer to?

Another thing: In the example above there are 7 positive values follow by 2 negative values. So the membran vibrating in a speaker to create the sound is pressed longer in one direction than in the other direction, or how does one has to interpret this? Positive values are one direction of the vibration and negative ones are vibrations into another direction or am I wrong about this?

Thank you very much :-)


texts: How can I sample a signal at 4 samples per cycle?

texts:
I am in Chapter 5 of "Understanding Digital signal processing" 3rd edition by Richard G. Lyons.
I am stuck with Section 5.4.

In order to get the bandpass filter centred at $F_s/4 \mathrm{Hz}$, I need to multiply the low pass filter coefficients by a sinusoid of $F_s/4 \mathrm{Hz}$ samples at 4 samples per cycle. I have tried to do this in MATLAB but cannot figure out how to sample at four samples per cycle. Since samplerate was not given, I created a sinusoid:
$$
n=0:30;\\
x=\sin(n\pi/2)%
$$
Is this correct for $F_s/4$ signal?

How do I sample it at four samples per cycle in MATLAB?

I feel like I missing something very fundamental.


texts: Rayleigh fading with frequency selective fading channel

texts: I Read in one reference that, the Rayleigh distribution is commonly used to describe the statistical time varying nature of the received envelope of a flat fading signal, or the envelope of an individual multipath component. Can Rayleigh and Rician distributions be used with frequency selective fading channels? I need some explanation please


texts:
The Polynomial Bandlimited Step (PolyBLEP) methodology involves approximating the sinc() function with a polynomial. The most basic polynomial to choose is a uni-polar triangular pulse as it is the linear approximation of the windowed sinc() pulse. Integrating it produces the bipolar integrated triangle which can be expressed as:
\begin{cases}0\,\textrm,&t<-1\\t^{2}/2+t+1/2,\,&-1\leq t\leq 0\\t-t^{2}/2+1/2,&0<t\leq 1\\1,\,&t>1\end{cases}
The PolyBLEP two part residual is then obtained by subtracting the perfect step from the bipolar integrated triangle and converting again to bipolar.

Residual is then applied to the naive wave shape to offset the samples around the discontinuities, shifting them up or down to producing a smoothed waveform.
How are those polynomials and two part residual actually, step-by-step, 'derived'?
I tried to solve this two part residual through those bipolar integrated triangles and ended up to different result:

which results two part residual as:
-((|t|-1.0)(-|t|+1.0))
((|t|-1.0)(-|t|+1.0))

See: https://www.desmos.com/calculator/xwqyghn8pb
Are my results OK/usable?


texts: DFT of signal with frequency-dependent time shift

texts:
I have a signal $s[n]$, and I want to compute the DFT of it with a frequency selective shift (computing all $k$ here but I really only care about $k > 0$)

$$S[k] = \sum_{n=0}^{N-1} s[n-\dfrac{N}{2k}] e^{(-j2\pi nk/N)}\space\space\space  -N/2 \le k \lt N/2 $$

Following the derivation of the general shift theorem here, I can make the change of variables $m = n-\dfrac{N}{2k}$:

$$S[k] = \sum_{m=-N/2k}^{N-1-N/2k} s[m] e^{(-j2\pi (m+N/2k)k/N)}$$

Carrying through the product in the exponential term:

$$S[k] = \sum_{m=-N/2k}^{N-1-N/2k} s[m] e^{(-j2\pi mk/N)} e^{-j\pi}$$

And realizing $e^{-j\pi} = -1$:

$$S[k] = -\sum_{m=-N/2k}^{N-1-N/2k} s[m] e^{(-j2\pi mk/N)}$$

Which, to my surprise, indicates that the shift drops out, though now I'm summing over shifted intervals.  I believe that's just another phase adjustment that has to be made though, given $s[n]$ is periodic.

I have a hard time believing this, I would have expected something like a quadratic phase shift in frequency to result.  Can anyone point to where I've made a mistake?


texts: A cubic interpolation function: folkloric copypasta or clever trade-off?

texts:
I've been reading on interpolation methods recently and I have come across an implementation of cubic interpolation that is leaving my head scratching. Every other variant and example of cubic interpolation I've come across has included lots of fractions and exponents. This implementation I found on Paul Bourke's website uses only additions and subtractions to figure the coefficients. I've seen this implementation used in several codebases but unfortunately without any documentation.
double Cubic(double x, double a, double b, double c, double d) {
    double A = d - c - a + b;
    double B = a - b - A;
    double C = c - a;
    double D = b;
    return A * (x * x * x) +
           B * (x * x) +
           C * x +
           D;
}

If I may, let me first walk through how this implementation might have come to be just to make sure I at least understand the algebra going on here. But I'm curious as to why it can be made so simple.
As I understand it, cubic interpolation takes four sample points, say $a$, $b$, $c$, and $d$, and interpolates between b and c using the function:
$$
f(x) = Ax^3 + Bx^2 + Cx + D \tag{1}\label{1}
$$
Where $0 \le x \le 1$.
A major condition I often see (one that's obviously important for me) is that $f(x)$ must run through the sample points. So, in this case, samples $b$ and $c$ correspond to $x = 0$ and $x = 1$ respectively.
$$
f(0) = b \tag{2}\label{2}
$$
$$
f(1) = c \tag{3}\label{3}
$$
There is another condition that says that the derivative at samples $b$ and $c$ are the same as the slope between their surrounding points. I believe this is called a Catmull-Rom spline?
$$
f'(x) = 3Ax^2 + 2Bx + C \tag{4}\label{4}
$$
$$
f'(0) = \frac{c - a}{2} \tag{5}\label{5}
$$
$$
f'(1) = \frac{d - b}{2} \tag{6}\label{6}
$$
Now, if I follow that path, I get Paul Breeuwsma's solution here and that all makes sense to me.
However, in order to get Bourke's implementation, I have to multiply the slopes at $b$ and $c$ by 2! That is:
$$
\begin{align}
f'(0) = c - a \tag{5a}\label{5b}\\
f'(1) = d - b \tag{6a}\label{6b}
\end{align}
$$
If I follow from that, I get the below for the coefficients $A$, $B$, $C$, and $D$ in order to match Bourke's.
First, $D$:
$$
\require{cancel}
$$
$$
\begin{align}
f(0) &= \cancel{A(0)^3} + \cancel{B(0)^2} + \cancel{C(0)} + D \\
f(0) &= D \\
D = b \tag{7}\label{7}
\end{align}
$$
Now, $C$:
$$
\begin{align}
f'(0) &= \cancel{3A(0)^2} + \cancel{2B(0)} + C \\
f'(0) &= C \\
C = c - a \tag{8}\label{8}
\end{align}
$$
Now, $B$:
$$
\begin{align}
f(1) &= A(1)^3 + B(1)^2 + C(1) + D \\
f(1) &= A + B + C + D \\
c &= A + B + c - a + b \\
B = a - b - A \tag{9}\label{9}
\end{align}
$$
Finally, $A$:
$$
\begin{align}
f'(1) &= 3A(1)^2 + 2B(1) + C \\
d - b &= 3A + 2(a - b - A) + c - a \\
d - b &= 3A + 2a - 2b - 2A + c - a \\
d - b &= A + a - 2b + c \\
A = d - c + -a + b \tag{10}\label{10}
\end{align}
$$
Ok, so... why the heck? This looks like nothing I've seen while reading up on this.
By doubling the slope at $b$ and $c$, besides having a fast and definitely appealing implementation, is this a trade-off or are there additional conditions that were able to be made that I missed? It seems like this could cause the interpolant to shoot around much more but I haven't yet coded up a comparison to look at and listen to. I've been scouring articles, blog posts, and papers and I can't seem to map anything to this. If anything, other implementations end up quite hairy and difficult for me to understand. Is there a name for this kind of interpolation? Where did it come from?
I know this was a long question. Thanks for reading!


texts: Phase plot of $H(e^{j\omega})=(1-re^{-j\omega})\left(1-\dfrac{1}{r}e^{-j\omega}\right)$, $0<r<1$

texts:
I want to find the phase plot of $H(e^{j\omega})=(1-re^{-j\omega})\left(1-\dfrac{1}{r}e^{-j\omega}\right)$, $0<r<1$ for the interval $0\leq \omega \leq \pi$.

Method 1:

$H(e^{j\omega})=1-\left(r+\dfrac{1}{r}\right)e^{-j\omega}+e^{-2j\omega}=e^{-j\omega}\left[2\cos \omega-\left(r+\dfrac{1}{r}\right)\right]$

$\angle H(e^{j\omega})=-\omega+\angle \left[2\cos \omega-(r+\dfrac{1}{r})\right]=-\omega+\pi$

Method 2:

$\angle H(e^{j\omega})=\angle(1-re^{-j\omega})+\angle(1-\dfrac{1}{r}e^{-j\omega})$

$$\angle H(e^{j\omega})=\begin{cases}
\tan^{-1}\bigg(\dfrac{r\sin \omega}{1-r\cos\omega}\bigg)+\pi-\tan^{-1}\bigg(\dfrac{\sin \omega}{\cos \omega-r}\bigg),0\leq \omega\leq \cos^{-1}(r)\\
\tan^{-1}\bigg(\dfrac{r\sin\omega}{1-r\cos\omega}\bigg)+\tan^{-1}\bigg(\dfrac{\sin\omega}{r-\cos \omega}\bigg),\cos^{-1}(r)\leq \omega\leq \pi\\
\end{cases}$$

I tried applying the formula $\tan^{-1}x+\tan^{-1}y=\tan^{-1}\bigg(\dfrac{x+y}{1-xy}\bigg)$, but I am confused when to use this formula and unable to proceed forward. Please help in solving this.


texts:
I have generated colored noise from the white gaussian noise using auto-regressive model (by solving Yule-walker equation).
 The covariance matrix that I used to generate the colored noise is [1 0.2;0.2 1].

It shows that variance(power) of white noise and colored noise are same. So if I don't whiten the colored noise, the BER performance of the white noise and colored noise should be same because power of both noise are same (1 in this case).

But the simulation shows that the BER performance of colored noise is nearly 1.4 dB degraded as compared to the white noise. It seems the covariance term (0.2) is causing this BER degradation.

Can anyone please explain why this happened? I used the matched filter as the original pulse shape in both cases. Does this affect the performance?

I am considering my threshold as: if the matched filter output at signal interval >0 then bit 1 is transmitted else bit 0 is transmitted.


texts: Nyquist Rate (Sampling Frequency) for $ {f}^{2} \left( x, y \right) $

texts:
We are given that $f(x,y)$ is highest frequency is $\omega$ what will be the frequency sample rate if we want to restore the function of the form $g(x,y)=f^2(x,y)$

Would it be correct to say that because $\sin^2x=1-\cos 2x$, we will have to sample it at $4\pi$ ?


texts: Step response of a differentiating system

texts:
Consider I can find the step response of a system with the following method:

We can easily find the step input of a system from its transfer function.  Given a system with input $x(t)$, output $y(t)$ and transfer function $H(s)$

$$H(s) = \frac{Y(s)}{X(s)}$$

the output with zero initial conditions (i.e., the zero state output) is simply given by

$$Y(s) = H(s)X(s)$$

so the unit step response, $Y_\gamma(s)$, is given by

$$ Y_\gamma(s) = \frac{1}{s} H(s)
$$

If my input is a "differentiating system" then what will be the step response ?


texts:
So I have an image.
I take it into HSV format and v has its intensity values.

now the intensity values of v vary from 0 to 1 (Black to White)

I want to divide v into 3 parts:

1] v1= with intensity values 0<v<0.3

2] v2= with intensity values 0.3<v<0.75

3] v3= with intensity values 0.75<v<0.1

How do I do this?


texts:
While studying about Wavelet Transform, I'm confused as why this is called a [1 -1] Filter , Here is what i'm talking about.

  In this hypothetical example the student does fairly well the first
  half of the     term then neglects his or her studies for the last
  half. Thus the exam scores     for the term were 80%, 80%, 80%, 80%,
  0%, 0%, 0%, and 0%*     We can tell the average of all the scores
  (40%) and when the scores “tanked”     after the 4th exam just by
  looking. Knowing the answer in advance, however,     is a good way to
  learn and to verify the wavelet transforms. Then we     can use them
  with confidence on real-world data where we can’t simply     “eyeball”
  the final values.

  We will now walk through the CWT process step by step using the
  simplest     of the wavelet filters on this example. We begin by
  comparing the humble Haar wavelet filter, [1 –1], with the data as
  shown

80 80 80 80 0 0 0 0

If We now keep subtracting each value from the next value we will eventually get,

[0, 0, 0, 80, 0, 0, 0]

Now If I want to stretch the Filter to three points i.e If the filter is stretched from [1 –1] to [1 0 –1].

How will I calculate these points as I did in case of two point Filter where i simply subtract each value from the next value ?


texts: Measuring frequency response range of piezoelectric disc

texts:
How can I measure the frequency response range of a piezoelectric disc below.

I don't have an oscilloscope but I do know how to program in matlab / octave I was thinking of creating a swept signal in matlab and playing it through my laptop or tablet and playing it through the disc.  But how do I measure the frequency response?


texts: What is the difference between wide sense and strict sense stationary processes?

texts:
What is the difference between wide sense and strict sense stationary processes (SP) ?

According to the definition (by Heinrich Meyr, Marc Moeneclaey, Stefan A. Fechtel in ￼"Synchronization, Channel Estimation, and Signal Processing") :

strict-sense SP = not time dependent
wide-sense SP = not dependent on variable t (time)

Could someone please illuminate what this really means ?

Perhaps pointing out a good book in which the difference is explained ?

If an SP is not dependent on variable t (time) then how can it depend on time ( a process which is not strict sense stationary but wide sense stationary)?

Could you give an example of such a process ?

Why and when is making this difference useful/neccessary?


texts:
"Rectangular window" and "Triangular window" have gibbs phenomenon
So which window should be used in fir filter design to have minimum or zero gibbs phenomenon


texts: The Wiener Filter for the Time Derivative of the Signal

texts:
Assume the observation of the original signal $s(t)$ is $x(t) = s(t) +n(t)$, where the signal and noise are independent.
Then we need to design the Wiener filter $g(t)$ to estimate $\frac{d s(t)}{dt}$ which is the derivative of the original signal $s(t)$.
The design of the Wiener filter is from the MMSE: $E[e]$. But how to formulate the error function here? $e(t) = \frac{d s(t)}{dt} - x(t)$ or $e(t) = \frac{d s(t)}{dt} - \frac{d x(t)}{dt}$?
I think this can be generalized to any linear operators like estimating $y = Hs$ (deconv, deblur etc.). But I have no idea the begining step. Could anyone give me some hints? Thanks!


texts:
I'm having troubles trying to use one relatively simple noise estimation method, here it is:A Fast Method For Image Noise Estimation Using
Laplacian Operator and Adaptive Edge Detection

My problem is with their proposed noise deviation formula (page $2$, ($3$)), I can't seem to comprehend it. What I want to know am I right that this $N$ is amount of pixels in image?

And also I'm having troubles implementing first part of the formula since its result is float undeflow. Maybe I can just use regular pixel_value - pixel_mean formula? Will it produce same results?

Altrough images with massive noise will have mean that actually pretty much depends on the noise amount and thus estimation will be considerably lower than the real value.

Where can I read about this formula and how to use, or how it's even named?

Maybe I've implemented it wrongly, here's my code that does it:
noise deviation implementation (C++ CImg)
It gives me very big numbers, like for noisy image of bears I got something like 1830152.

First of all, what I'm supposed to get from that formula? Isn't it supposed to be some kind of 0.xx coefficient or maybe percent of noise in image?


texts: How does a low-pass filter programmatically work?

texts:
I have been working on a simple low pass filter for < 100 Hz metering in my application. But so far, I am struggling with the theory behind it all. It's cool that I got it working, but I'd really enjoy it if I knew how/why it is working.

I found the following code:

void getLPCoefficientsButterworth2Pole(const int samplerate, const double cutoff, double* const ax, double* const by)
{
    double PI = M_PI;
    double sqrt2 = sqrt(2);

    double QcRaw  = (2 * PI * cutoff) / samplerate; // Find cutoff frequency in [0..PI]
    double QcWarp = tan(QcRaw); // Warp cutoff frequency

    double gain = 1 / ( 1 + sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) );

    by[2] = ( 1 - sqrt2 / QcWarp + 2 / ( QcWarp * QcWarp ) ) * gain;
    by[1] = ( 2 - 2 * 2 / ( QcWarp * QcWarp ) ) * gain;
    by[0] = 1;

    ax[0] = 1 * gain;
    ax[1] = 2 * gain;
    ax[2] = 1 * gain;
}

To calculate the Coefficients. Then, in the audio samples, I 'low pass' them this way:

        xv[2] = xv[1];
        xv[1] = xv[0];

        xv[0] = pData[j];
        yv[2] = yv[1];
        yv[1] = yv[0];

        yv[0] = (ax[0] * xv[0] + ax[1] * xv[1] + ax[2] * xv[2]
                   - by[1] * yv[0]
                   - by[2] * yv[1]);

        pData[j] = yv[0];

To get a lowpass design.

I am wondering a few things:

I receive the audio samples in a simple float* array. What is that float number? The only thing I see is a number, how is that a piece of sound?
The code is using past calculations (three of them) in the new calculation per sample. Does that mean that the first 2 data samples are not filtered correctly? (not that it would matter because it's just 2 samples, but just wondering)
Trying to learn everything, I found a couple of formulas for the Butterworth (2nd Pole) filter. How are those formulas reflected in this code? None of the formulas I found have these calculations that you can see in the 'getLPCoefficientsButterworth2Pole()' function.


texts: Convert sine wave signal to square wave signal

texts:
I have a signal containing a variation of digital ASK type on-off keyed sine wave among which some stochastic irrelevant harmonic waves occur, like you can see on this figure:

In its current form each pulse duration is the same, but it's not always the case.

I'd like to calculate each pulse duration. To do that, I was thinking of converting this signal to a square wave signal like this:

I'm not sure of the different steps needed for this conversion but I have some leads:

Remove the noise => Low pass filter I presume?
To have the signal between 0 and the amplitude, I multiply the signal by itself
I apply a magic filter to convert the signal to square signal.

So my questions are:

Do I need a low pass filter to remove the noise?
What kind of filter do I need to convert the signal to square signal? Is it possible?

For information, I'm using Python and Numpy.

If there is any other method to calculate the duration of each pulse in the sequence, I'm open to all ideas.

Thanks!

Edit: Here's my signal:

   freq = 1000000
   fs = int(2*freq) # sample rate
   x = np.linspace(0, 10, fs)
   noise1 = 0.01*np.sin(2*np.pi * 1*x) + 0.05*np.sin(2*np.pi * 1.5*x) +     0.04*np.sin(2*np.pi * 0.2*x)
   noise2 = 0.1*np.sin(2*np.pi * 0.1*x) + 0.2*np.sin(2*np.pi * 1.8*x) + 0.1*np.sin(2*np.pi * 0.4*x)

  realSignal = 2*np.sin(2*np.pi * freq*x)
  sig =  noise1

  sig = np.append(sig, [realSignal, noise2, realSignal, noise1*noise2])


texts:
I'm quite new to this; I'm try to classify textures as defective or non-defective. I've used a Gabor filter bank with Matlab which outputs a column vector of the Gabor features of an image. I have a data set of non-defective images and defective images.

My question is, what can I now do with this (or these) feature vectors to classify the texture? I've read about many types of classification (SVMs, Neural networts, nearest neighbour etc.), but couldn't find any similar types of implementation to help me get an idea of what I'm doing - perhaps someone could point me towards a practical example to help my understanding. Many thanks.


texts:
First of all, I have to stress that I am not a professional of coding, no more than a professional of signal processing. I am a chemist that happen to be working on a project involving both. So in advance, please forgive me if my questions are not well worded, or seem trivial to you (in addition, I am not a native english speaker). I have several questions that are not directly related to each other, so I might open several threads as I got answers.

So, I am working on a project involving the interference patern of a laser, and more importantly, shifting of this pattern alongside the horizontal (x) axis. Here is the typical diffraction pattern I observe through my camera :

My goal is to precisely measure the shift of that pattern that happens during my experiment. For a better understanding of what I mean, I attach an exemple with a signal shifted from the other.Keep in mind that the shift presented here, that is of about 40 pixels, is considered HUGE, and will never really happen during my work, I am typically observing shift of about 0-5 pixels.

To achieve this goal, I wrote a C# code that computes the cross-correlation between two arrays (I reduced my image to a single 1024 pixels line, by calculating the mean of the signal alongside each 768 pixels row). The cross-correlation has been computed using this forumula (using the Accord Framework):

corr(a, b) = ifft(fft(a_and_zeros) * conj(fft(b_and_zeros))).

The code seems to work, here is the kind of cross-correlation figure I get :

This particular figure is obtained when cross-correlating two identical images (I guess I could call it auto-correlation, but I am not sure, and, as expected, I find a shift of 0 pixels). Keep in mind that this figure is the product of a modification of the cross-correlation results array, as applying Fast Fourier Transform results in the peak being "cut in half" at the edge of the array.

My problem arise when I zoom on the peak :

As you see, the peak is sharp, and I think it should be smooth. This is a problem when I want to observe low shift, as the "sharpness" of the peak hides the real maximum :

The problem disappears when the shift becomes bigger, but then the presence of a really sharp peak I can't explain is exposed :

Note that this peak doesn't come from me modifying the array to center the peak, it is present at the edge if I don't modify the results.

This ends up with my results not being realistic for low shifts (AKA the shifts I care about) as shown here :

As you can see, things are not exactly on point (this is the very start of the project, so the experimental setup is far from perfect, as is the code), but the tendancy is there, and I would really be able to solve this peak issue before going further.

So I come to you in hope of finding out why my experiment are behaving this way, and what I could do about it ! Is there a reason linked to the way cross-correlation works, or something else ? I couldn't find any issue in the code, as of now.

Thanks a lot !


texts:
I have written Matlab script for BPSK modulation. The way my code works is as follow:

series of 0 and 1----> NRZ data (continuous rectangular pulses with amplitude -1 and 1)--->multiply by carrier----> add noise-->  multiply by carrier---> integration--->detection---> Bit Error collection?

I would like to know can anyone tell me where I should apply DAC and ADC ?
I am wondering should it be after NRZ data or before that for DAC ?
:!:


texts: Where is the mistake in the convolution?

texts:
Let :
$$
x[n]=\begin{cases}
1&\text{if $0\leq n\leq 4$}\\ \\
0&\text{if otherwise}
\end{cases}
\qquad
\text{and}
\qquad
h[n]=\begin{cases}
\alpha^{n}&\text{if $0\leq n\leq 6$}\\ \\
0&\text{if otherwise}
\end{cases}
$$
We must compute the convolution :
$$
(x*h)[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]
$$
Observe that $x[n]=u[n]$ for $n\in[0,4]$ and $0$ otherwise. Furthermore, $h[n]=\alpha^{n}u[n]$ for $n\in[0,6]$. Thus :
\begin{array}{c | c c c c c c c }
\mathbf{x} & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
\hline
1 & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
1 & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
1 & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
1 & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
1 & 1 & \alpha & \alpha^{2} & \alpha^{3} &\alpha^{4} &\alpha^{5} &\alpha^{6}\\
0 & 0 & 0 & 0& 0& 0 & 0 & 0 \\
0 & 0 & 0 & 0& 0& 0 & 0 & 0
\end{array}
\begin{align*}
\displaystyle (x*h)[n]&=\left[\sum_{k=0}^{0}\alpha^{k}\;,\;\sum_{k=0}^{1}\alpha^{k}\;,\;\cdots\;,\;\sum_{k=0}^{4}\alpha^{k}\;,\;\sum_{k=1}^{5}\alpha^{k}\;,\;\sum_{k=2}^{6}\alpha^{k}\;,\;\cdots\;,\;\sum_{k=5}^{6}\alpha^{k}\;,\;\sum_{k=6}^{6}\alpha^{k}\right]\\ \\
&=\begin{cases}
\displaystyle\sum_{k=0}^{n}\alpha^{k}&\text{if $0\leq n\leq 4$}\\
\displaystyle\sum_{k=1}^{n}\alpha^{k}&\text{if $n=5$}\\
\displaystyle\sum_{k=n+2}^{6}\alpha^{k}&\text{if $0\leq n\leq4$}
\end{cases}
\end{align*}
I wish to know where my mistake is and how to fix it. Thank you :)


texts: Design of Parks-McClellan Lowpass Filter

texts: Can we design Parks-McClellan Lowpass Filter with only sampling frequency of 12KHz and cutoff frequency of 4KHz and order of 30? I went through Matlab help, where we need parameter like passband cutoff frquency, stopband cutoff frequency, passband ripple e.t.c.,


texts:
Simple math question. The convolution theorem states that multiplication in time domain is equal to convolution in frequency domain and vice versa. There is a condition that the signal has to be properly zero padded as to not cause aliasing.

This question concerns convolution in the frequency domain. The difficulty arises from the fact that we are dealing with a complex signal which has a positive and a negative side. Let's consider two frequency domain signals, with negative frequencies included, presented in polar coordinates []:

a=[0,1(3π/4),0,1(π/2),0)]
b=[1,1,1,1,1]

So we have a sine wave and a zero-phase dirac delta. It's apparent that the time domain product is zero. But coming to this conclusion via using convolution in frequency domain doesn't seem to be so simple. I get (rectangular):

[-j,-j,-j,-j, 0,j,j,j,j,j]
=> [-2j,-2j,0,2j,2j]

Which is not correct. Even for a cosine I get a row of twos, which is not correct either (should be row of ones). Where am I going wrong?


texts: What is the moment when all oscillators aligned to make a jump called?

texts:
Say we have a square wave and its Fourier series. When the wave doesn't jump, its oscillators aren't aligned:

But if they are aligned, the wave will jump:

What is this moment called? They might not be necessary in the same angle, and I don't think "resonating together" is correct, since they have different frequencies. What is it called when all oscillators aligned to make a jump?

Meta discussion: Why are questions asking for terms off-topic?


texts: Change pitch of mp3 audio on a microcontroller

texts:
I have an MP3 decoder running on a ESP32, and I'd like to change the pitch of the streaming audio in real time.

I plan to build a little toy that streams Amazon Alexa and changes her voice so that it is higher pitched. So no pro audio requirements at all, it is just supposed to be a funny effect.

I was thinking maybe I can operate on the raw DCT data? Or is there another way that is fast enough for that chip? I cannot use time compression/dilation because the audio is streamed.

I'm using the MAD fixed point decoder.

A description of (or a link to) an algorithm working on DCT or PCM data would be fine, but code or pseudo-code in whatever language would be very welcome.


texts: Synchrosqueezing Wavelet Transform explanation?

texts: How does Synchrosqueezing Wavelet Transform work, intuitively? What does the "synchrosqueezed" part do, and how is it different from simply the (continuous) Wavelet Transform?


texts: How can poles and zeros exist at infinity?Can anybody explain using a system function?


texts: What is the absolute phase and the relative phase, of a signal?

texts:
I need to know what is the absolute phase and the relative phase, of a signal?
and why these phases are important in the signal processing?


texts: What is the difference between object detection and target detection? Are they the same thing?

texts: I was reading some articles regarding object detection and target detection and have a hard time understanding the difference between object detection and target detection. i want to know whether they are different from each other? and if so, how?


texts: how sinc pulses can be used to distinguish different symbols

texts:
Let's say we want to transmit a sequence of numbers Using PAM, this sequence is transmitted with the signal $s(t)$
The pulses in $s(t)$ are transmitted at a rate $R_p=\frac{1}{T_p}$. This is called the "pulse rate" or the "symbol rate". We want $T_p$ to be short, to transmit information faster. However, this increases the bandwidth. The maximum spectral efficiency can be achieved using sinc pulses. I have the following questions please:

How can sinc pulses be used to distinguish the different numbers? If we are using 16-QAM, then we have 16 different symbols. Is it the amplitude of the sinc pulse that will change to distinguish the 16 symbols??
The bandwidth of a single sinc pulse is $\frac{R_p}{2}$, what is the bandwidth of the whole signal $s(t)$ consisting of several pulses, each of width $\frac{R_p}{2}$?


texts: About the post "How to use “Filter” matlab function to apply channel to signal?"

texts:
I would like to add my question to the post How to use “Filter” matlab function to apply channel to signal?, but I don't have permission to  do that.

The user @Anna Koroleva discussed with @Enginner the Matlab function "rayleighchan". Usually, I am using outputs of the function as filter(h,s) but there they discussed that it is possible to fill outputs in the matrix.

How to fill PathGains in the matrix? Could you please give me an example of how to do it?

EDIT 1:

I have wrote a script to simulate user @Engineer suggestion :

N=4; Nr=2; Nt=2
h = rayleighchan(Ts, Fd, tau, pdb);
chan = zeros(N,Nr, Nt);
ofdm_chan = zeros(N,Nr, Nt);
for k=1:N
    chan(:,k) = filter(h,ones(N,1));
end

But it doesn't work. Should I implement in another way?


texts:
How effective is Nash-equilibrium for big data cases in scenarios like crowd and traffic monitoring?
Crowd and traffic scenarios are competitive environments which suit well to Nash-Equilibrium. Are there algorithms using Nash-Equilibrium which are successfully deployed?


texts: Which one is the correct graph for $u[2-n]$?

texts:
$u[n]$ is the unit step function in discrete form. I want to know which one of the following two is the correct waveshape for $u[2-n]$, where $k$ is a constant.

What's the answer? The top one is $u[n]$, the unit step function. So, what's the graph for $u[2-n]$, middle or bottom?


texts: Plotting the spectrum of a continuous sinusoid function

texts:
I am just learning about spectrums and I am a little confused. I have given a simple example and where I get stuck.

  Plot the spectrum of  $x(t) = 5\cos\left(800\pi t + \frac{3 \pi}{7}\right)$.

So I use reverse-Euler's to write this as
$$x(t) = \frac{5 e^{i \frac{3 \pi}{7}}}{2} e^{i 800\pi t} + \frac{5 e^{-i \frac{3 \pi}{7}}}{2} e^{-i 800\pi t}$$

What is the spectral amplitude in each case? I am unclear if the frequency of $800 \pi t$ has amplitude $5$ or $\frac{5 e^{-i \frac{3 \pi}{7}}}{2}$.
If it is the second do I make two seperate plots for the complex and real parts, or do I just plot $\bigg\lvert \frac{5 e^{-i \frac{3 \pi}{7}}}{2}\bigg\rvert $?


texts: Fixing drift in a recursive fixed point filter

texts:
I'm implementing a 80-72-64-48 multi pass moving average filter for a embedded system in C and in fixed point. The implementation is a circular buffer where i'm keeping a running sum and calculating

$$ y[n] = y[n-1] + x[n] - x[n-M] $$ where $M$ is the length of a filter. This is done for each subfilter with the output from one serving as the input for another.

I'm scaling my coefficients by $2^{15}$ which gives me coefficients of length $2^{9}$ or $2^{10}$ depending on the filter length. Then the result is scaled down again by $2^{-15}$ to get the correct output.

Now, everything looks good on short time scales but over long times i get a drift. The reason for the recursive implementation is to save calculations on a embedded processor.

I have included image of some of the internals of my filter, this is when a step response is applied and we can see the transfer functions of the filters taking shape, square, triangle, then approximating a gaussian so the filter is working as expected.

Now here is my problem, the  drift

Is there any way to fix this, and where is the most likely source of this?, is this drift due to a bit getting lost in the shifting or something else. The drift is not present for DC inputs, but for AC signals it slowly drifts.

SOLUTION: The problem was in the accumulator as robert suggested in the comments. The issue was that one element of the calculation had gone through an extra up and down shift compared to the rest, which created a round offset that accumulated.


texts: As of 2019, which discrete nonlinear, time-invariant systems with memory are considered "easy" to model and identify?

texts:
There are several types of discrete nonlinear time-invariant systems with memory ("NTIM") which are considered "easy" to model and identify. Any such system can be represented using a Volterra series, but Volterra kernels become very difficult to model for any more than quadratic or perhaps cubic nonlinearities.

However, for certain types of NTIM systems, there exist clear methods to completely measure and replicate some arbitrary unknown system, typically via measurement of certain test signals.

One good example are those systems for which each Volterra kernel has nonzero coefficients only on the diagonal. These are sometimes called "Farina," "diagonal-Volterra," or "Generalized Hammerstein" systems (Fig. 1). Given some suitable smoothness assumptions, these can be perfectly modeled using dynamic convolution, which is a variant of regular convolution that has a different impulse response for each scaled impulse. As a result, an arbitrary system of this type can be perfectly characterized by simply measuring the impulse response for each scaled impulse. Some info on that in this related answer: Dynamic convolution vs Volterra series. (There are also techniques involving synchronized sine sweeping)

Figure 1. A diagonal Volterra series system, which has nonzero Volterra coefficients only on the diagonal of each Volterra kernel, is equivalent (given suitable smoothness assumptions) to a sum of cascades of a memoryless waveshaper (illustrated as a sigmoid) and a linear time-invariant (LTI) system.

I am wondering if there exists any larger set of nonlinear systems that is also considered easy to model. There also exist the Wiener (Fig. 2) and Wiener-Hammerstein systems (Fig. 3) , which are not diagonal-Volterra; are those easy to model? And all of these systems seem to be subsets of "NARMAX" systems; how easy are they to model? Is there any easily modelable type of system which generalizes both the Weiner-Hammerstein and diagonal-Volterra systems?

Figure 2. A Wiener system is a cascade of an LTI system and a memoryless waveshaper.

Figure 3. A Wiener-Hammerstein system is a cascade of an LTI system, a memoryless waveshaper, and another LTI system.

Questions

As of 2019, which nonlinear time-invariant systems with memory are clearly easy to model and profile?
Is there any type of nonlinear system, which includes all the Volterra-diagonal ones as a subset, which is easy to model?
Is there some type of system which includes both the Volterra-diagonal and the Weiner-Hammerstein systems, and is easy to model?
All of these systems are a subset of the "NARMAX" nonlinear models - are those easy to model?


texts: How to take the inverse FFT of Welch's power spectral density in Python?

texts:
I have calculated the power spectral density of a weakly-stationary signal using Welch's method (as implemented by scipy.signal.welch in Python), and am looking to take the inverse Fourier transform of this PSD to find the autocorrelation of the function.

(The reason why I'm not calculating the autocorrelation from the signal directly is because I'm applying transformations to the original PSD, which will yield autocorrelations different than that of the signal.)

signal.welch returns two arrays: one of positive frequencies and the other of the corresponding amplitudes. I figured from this that one could apply numpy.fft.irfft to the amplitude array to find the autocorrelation of the signal behind the PSD.

I am having trouble accomplishing this, however. With an amplitude array a from signal.welch, my method has been to apply irfft(a) and plot over the time-domain of the original signal. This is not yielding what I'm expecting (as there is appearing an unexpected symmetry in the time-domain, among other artifacts).

My question boils down to: how does the inverse FFT (calculated using fft.irfft) of an array a of positive frequency amplitudes calculated from signal.welch map to the time-domain of the original signal?


texts: How to design a FIR low pass filter and force a subset of filter coefficients to 0?

texts:
I'm looking into interpolating a signal by 2x and came into this matlab function interp:

  interp

  Interpolation — increase sampling rate by integer factor

  Algorithms

  interp uses the lowpass interpolation algorithm 8.1 described in 1.

  1.It expands the input vector to the correct length by inserting 0s between the original data values.

  2.It designs a special symmetric FIR filter that allows the original data to pass through unchanged and interpolates to minimize the
  mean-square error between the interpolated points and their ideal
  values. The filter used by interp is the same as the filter returned
  by intfilt.

  3.It applies the filter to the expanded input vector to produce the output.

Then about intfilt:

  intfilt

  Interpolation FIR filter design

  Algorithms

  The bandlimited method uses firls to design an interpolation FIR
  filter. The polynomial method uses Lagrange's polynomial interpolation
  formula on equally spaced samples to construct the appropriate filter.

Finally the documentation of firls is here

In short, the algorithm of interp is to insert l-1 zeros between each pair of neighboring input samples (up-sampling) and FIR filter it, where l is the multiplier of interpolation. The FIR filter is designed by intfilt.

E.g. intfilt(2,4,0.5) returns a 15-tap FIR filter with 0.5*fs/2 cut-off for 2x interpolation, in this case:
    [-0.0068, 0.0000, 0.0395, -0.0000, -0.1427, 0.0000, 0.6098, 1.0000, 0.6098, 0.0000, -0.1427, -0.0000, 0.0395, 0.0000,-0.0068]

Observation:

In the above example: all the even coefficients are 0 and the middle coefficient is 1. Recall that this FIR filter is to be used with a 2x up-sampled signal, this means the "original" sample in the input will be unchanged (only multiplied with the center 1), while all the inserted 0s will the calculated by the FIR filter.

This offers two useful properties: 1) computational efficiency: "original" samples are unchanged; 2) lossless-ness: interpolation then down-sample leaves the original signal unchanged.

Question:

I found the source code for intfilt and it calls h=firls(n-1,F*2,M); with n=15, F = [0, 0.1250, 0.3750 , 0.5000] and M=[2, 2, 0, 0], then returned h is the same as h = intfilt(2,4,0.5);

Why this n/F/M combination would set all even coefficients of h to 0?


texts: Describing a time delay at a maximum for maxima and minima but tends to 0 elsewhere

texts:
I have two force-time signals obtained from a couple of force plates, corresponding to input/output signals for some system. Plotting the two together one of the first things I noticed was a time delay. Initially I used cross correlation to gain some approximation for this feature however now I look a little closer I realize things are little more complex.

Basically the output trails behind the input with a time delay that is maximum at the maxima/minima on the two signals and in between the delay tends towards 0.

I literally have no idea how to express this. Can anyone suggest anything?

EDIT: Rough plot below. Genuinely apologize for the MSPaint... I don't have access to anything better or my plots right now! Also the absolute integral of both signals is approximately equal.

I'm thinking of making the time delay dependent on amplitude but sadly my mathematical skills at the moment restrict me to something as basic as pulling off maxima/minima with basic differentiation and recording the separation of corresponding features as a function of amplitude. Far too crude!


texts: Sinc Interpolation implementation in Matlab for reconstructing or interpolating?

texts: I am working on my project and I need to have a Sinc Interpolation function for a real world image. But, it seems that I have not caught the concept of Sinc interpolation. Could some one help me with this issue? this function is for adding some value among samples or it is used for reconstructing? How can I use this function for reconstructing signal.


texts:
I am learning Wiener filtering for 2D images by myself. From one book, it reads

It should be noted that the study of the more general problem of signal denoising dates
back to at least Norbert Wiener in the 1940s. The celebrated Wiener filter provides the optimal solution to the recovery of Gaussian signals contaminated by AWGN. The derivation
of Wiener filtering, based on the so-called orthogonality principle, represents an elegant solution and the only known situation where constraining to linear solutions does not render
any sacrifice on the performance.

It mentions Gaussian signals. I know Gaussian noise. What is Gaussian signal? I cannot find the definition of Gaussian signal via google. Does it mean one signal is contaminated by Gaussian noise?


texts:
I am working on plotting range-velocity heatmap of the receiving FMCW signal. I actually have already plotted. However, I do not really understand the mapping relationship between Doppler FFT bin and velocity. I learned to map from 1d FFT bin to range with the following codes as the x-axis.

x = (0:ADCsample-1) - ADCsample/2;
x = x*Fs/ADCsample;
% freq to dist;
x = x*3e8/2/slope;

Some notes, I plot the 2D FFT with nFFT = No of chirp in one frame =64.


texts: How to Calculate Local Spatial Variance of an Image?

texts:
For those who work often on image processing, I'm facing an algorithm of Non-Uniformity Correction (NUC) that requires the calculation of a parameter called "Local Spatial Variance". I can't find a definition for it.
Does anyone have any idea about it ?
The article is called A Novel Infrared Focal Plane Non Uniformity Correction Method Based on Co Occurrence Filter and Adaptive Learning Rate, IEEE Access, 2019.
The "Local Spatial Variance" is on the 3rd page if you please.


texts: Why is an RC circuit a first order system?

texts:
I am going through the following tutorial: Time-Response Identification of a Resistor–Capacitor (RC) Circuit

Why do we call the RC circuit a first order system?

I understand the order of a system is based on the highest derivative involved in its governing equation. The input governing equation of an RC circuit is described by equation 1 in the tutorial:

$$e_i-iR-\frac{1}{C}\int i \;dt =0$$

I understand the equation. But, I'm confused because there is no derivative term in it and yet it is called a first order system.


texts:
I want to know more about the latest attempts to create synthesizers for singing. I would ideally like to know more about both the most widely respected algorithms/models as well as mainstream software or programs that implement them (whether free or proprietary). Can someone offer me resources that would give me a decent grasp of the above?

I chose to post this question on Signal Processing Stack Exchange since I know such software relies to some degree on signal processing type methods for analyzing or creating waves.

I'm very new to this area of research so any suggestions for where to start would be welcomed.

Thanks.


texts: What is 'Normalized frequency in the range [0,1)', à la DTMF & Goertzel algorithm

texts:
I have a functioning DTMF program. It takes a larger signal and breaks it up into many small components, then analyzes each segment with an FFT to determine if the magnitude of the response for a particular "bin", corresponding to a specific frequency, is present, thereby indicating that the number which is respresented by that frequency has been pressed, straight forward enough, but I want to implement the Goertzel Algorithm of the library jmathstudio and I can't quite figure out how to go about searching for the frequency I want. the documentation says this:

I construct the bins like this, the term '200' is dictated by the fact that my FFT is currently of size 400, and 2000 is determined by the fact that the sample rate is 4000.

int sixNineSeven = (int)(697.0*200/2000);
int sevenSevenZero = (int)(770.0*200/2000);
int eightFiveTwo = (int)(852.0*200/2000);
int nineFourOne = (int)(941.0*200/2000);

int twelveZeroNine = (int)(1209.0*200/2000);
int thirteenThirtySix = (int)(1336.0*200/2000);
int fourteenSeventySeven = (int)(1477.0*200/2000);

I guess the main thing I need to do is this step:

Normalized frequency in the range [0,1)

But does that mean that I can totally ignore the sample frequency?

I suppose I would be doing this Geortzel operation in lieu of an FFT, so should I just search for this normalized frequency in each of the small segments?

But when I do this:

   Complex coeff = su.goertzelFrequencyAnalysis(new Vector(c), (float)(1/697));

   System.out.println("coeff "+coeff);

My output is just a bunch of garbage like this:

   coeff org.JMathStudio.DataStructure.Complex@5e228a02
   coeff org.JMathStudio.DataStructure.Complex@7bd63e39
   coeff org.JMathStudio.DataStructure.Complex@2e8f4fb3
   coeff org.JMathStudio.DataStructure.Complex@42b988a6


texts:
I am measuring ECG signals on subjects and I am having issues filtering a very small signal from the inherent noise of the collecting device. I am using an 18 bit ADC to digitize the signal
Here is an example of the noisy signal I get:

I would like it to look like this:

I have tried different filters and the best signal yet has been from the 8-point mean averaging of the signal:

However, I am trying to find the p-wave, which is a small upward deflection just before the large QRS complex of the pulse. This p-wave is buried in the noise and I can't seem to get it out. I have tried creating a median beat of 300 heartbeats, and by doing this I can show that the p-wave does exist, but I need to show that it exists for every heartbeat.
Median beat:

Is there a more elaborate way to tease out the p-wave using the median beat, or could you point me to other strategies for improving the signal to noise ratio?


texts:
There exist many, many analog and digital data transmission protocols. It is often difficult to tell which mode one listens to and tedious to do for many signals.

In Ham Radio, a Reed-Solomon Identifier may precede the signal, but this is not the case for every signal in existence.

Are there algorithms that automate the detection and mode differentiation from raw signal data that don't include metadata such as RSIDs?


texts: Extended Kalman Filter (EKF) Linearization of Non Linear Functions

texts:
I have more general question about the extended Kalman filter usage. What is not clear to me is why the EKF uses non-linear functions $f$ and $h$ for state prediction and estimate, while in other places the Jacobian of these functions is used.

Why the following is never used?

First calculate the linearized state and measurement models at previous estimate point using Jacobian. Use the linearized state transition and measurements matrix everywhere instead of non-linear in this specific iteration.


texts: Spatial spectrum of EEG data: non-uniform DFT?

texts:
I want to compute the spatial spectrum of EEG data collected using a non-uniformly sampled grid  of sensors (as  in the figure below). One way to do this would be to interpolate the data on a uniformly spaced rectangular grid. However, I would like to try a non-uniform DFT approach.

QUESTION 1: Do you think an nonuniform DFT makes sense for such an electrode grid?

QUESTION 2: The reason I would like to used the nonuniform DFT is that I would like to extract the most parsimonious frequency representation. I will try to explain it better. With the uniform DFT, one can compute ~N/2 Fourier coefficient (N = number of points) at evenly spaced frequency intervals. One can estimate the DFT at further frequency point, theses, however, will only be an interpolated version of the set of N/2 coefficients. For the non-uniform case it is unclear to me, how many non-interpolated coefficients one can extract and what these coefficients are. Any idea?

QUESTION 3: I was checking the NFFT implementation, however I am confused by the fact that the grid points must be constrained to the [-1/2 1/2] grid. Is it enough if I scale down my grid to this range?

QUESTION 4: The electrodes are actually placed on the scalp surface. The grid below is a 2-D representation of the 3D scalp surface. Is there a way to compute a nonuniform DFT directly on the 3D grid?


texts:
All these three terms are often used in communication texts and apparently they are some how similar
But what is exact difference between each of them?


texts: Commenting frequency spectrum of a sampled data

texts:
the question below asked in my latest signal processing exam. I have no idea about it and i could not find any solution in literature. Could you please help me finding out the answer. Actually i believe that i am quite good at signal processing until this question.

A voltage signal is sampled digitally in the laboratory at three different sampling frequencies, and FFTs are performed on the sampled data. In all cases, 512 data points are taken. Here is a record of the observations:
‒ Data sampled at 115 Hz: A peak appears in the frequency spectrum at 2 Hz, with amplitude close to 3.0 V. Another peak appears in the frequency spectrum at 45 Hz, with amplitude close to 1.0 V.
‒ Data sampled at 280 Hz: A peak appears in the frequency spectrum at 48 Hz, with amplitude close to 3.0 V. Another peak appears in the frequency spectrum at 95 Hz, with amplitude close to 1.0 V.
-Data sampled at 500 Hz: A peak appears in the frequency spectrum at 185 Hz, with
amplitude close to 1.0 V. Another peak appears in the frequency spectrum at 232 Hz,
with amplitude close to 3.0 V.

Question: What are most likely the true frequencies (and their corresponding amplitudes) in the voltage signal.


texts: Why is distorted / quantized white noise still white noise?

texts:
White noise can be created by setting every sample in the time domain to a random number with a gaussian distribution. Every sample is random, so the frequency response is obviously random. Easy enough.

Next step: Set every sample to a random number between 0 and 1, but the numbers are equally distributed instead of using a gaussian distribution (please ignore the DC offsets, they don't really matter). Samples are random enough, so it still makes sense that it would sound like white noise.

Now, set every sample to either 0 or 1 with a 50/50 chance. This feels like it shouldn't sound like white noise anymore, but it still does.

Finally, set every sample to either 0 or 1 with a 90/10 chance. How does this still sound like white noise when it looks nothing like it?

EDIT
I did some experimenting (with sample rate 48000 Hz).
1/2 chance is basically white noise.
1/4 chance sounds almost like white noise, but if you know what you're looking for it's easy to tell the difference.
1/16 chance sounds like rain smattering against something (you can hear individual impulses), with white noise in the background.

Varying the amplitudes of the impulses did not seem to have any effect on the sound, other than volume.


texts: Physical significance of Power spectral density of sum of correlated random processes

texts:
For two processes $X(t)$ and $Y(t)$, PSD of $Z(t) = X(t) + Y(t)$ is
$S_Z(\omega) =S_x(\omega)+ S_Y(\omega) + S_{XY}(\omega)+S_{YX}(\omega) $.

If $X(t)$ and $Y(t)$ are orthogonal, PSD of the sum process is sum of the PSDs of individual processes. However for correlated processes, an additional PSD term comes into picture.  What is the physical significance of this behavior?
To be more specific, the resulting PSD is greater than the sum of individual of PSDs.How can we explain this behavior qualitatively?


texts:
In the paper, The averaging method for finding exactly periodic dispersion-managed solitons
, there is an expression for the chirp $C$ of a signal given by-
$C = \frac{\text{Im}\int_{-\infty}^{\infty} u^2u_t^2 dt}{|u|^4dt}$
I am writing a program to calculate the chirp on MATLAB. The chirp function I have written so far is-
function c = chirp(signal)
u = signal;
ut = gradient(signal);
utc = conj(ut);
numval = imag(sum(u.^2 .* utc.^2));
denval = sum(abs(u).^4);
c = numval/denval;
end

However, this is giving incorrect answers. Any help pointing me in the right direction would be much appreciated.


texts: MATLAB produces too short output for frequency modulation

texts:
I have a voice data contained in a vector. Vector length is 118796 and sampling frequency is 8000 Hz. So it's around 15 seconds of voice data.

When I doing frequency modulation with fmmod (Carrier 138 MHz, out sampling frequency 5*Carrier, freq. deviation 12.5 KHz) I getting the vector of the same size.

[voice fs] = audioread('female.wav');

Fc = 138 * 10^6;
Fs = 5*Fc;
frequency_deviation = 12.5 * 10^3;
modulated = fmmod(voice, Fc, Fs, frequency_deviation);

I thing output vector should be significantly longer because it represents much higher frequency but have the same "real time" duration of 15 seconds.
Can anyone explain why matlab produce such output and how to get correct realtime output vector? I only assume to make input signal interpolation to make input vector as long as output vector should be.


texts: Are there amplitude and phase characteristics of harmonic signal components?

texts:
Nonlinear amplifiers are characterized by their AM/AM and AM/PM. That's what I thought. But now I have a question: The amplitude and phase distortion characteristics, do they only belong to the fundamental signal component?

I mean, to fully characterize my amplifier, are there AM/AM and AM/PM characteristics for the harmonic components as well?


texts: Determine surface homogeneity (luminance, colour, shape) by fourier transform

texts:
Let's say I have a computer screen, or a TV, or a flat panel ceiling lamp, i.e. a certain luminance spread over an area.

If I wanted to determine, how homogeneous this source of light is, I'd take a picture of it and put it through a fourier transformation, e.g. a DFT.

If the panel were perfectly homogeneous, there'd only be one single peak at the zeroth frequency and nothing else. Of course, they never are, and the calculations are done in a discrete system, so it will resemble a 1/x function.
Depending on the amplitudes at frequencies > 0, the panel will be homogeneous to a certain degree. High frequencies might indicate sharp edges, circular blobs in the actual image might be represented by circles in the amplitudes spectrum as well.

From such an amplitude spectrum alone, it would be very difficult to get a sharp definition for "good" or "bad".

As far as I know, the phase spectrum of an image contains most of the information for reconstructing the actual image. But looking at the phase spectrum, I see just "noise". Would it be possible to use the phase information as well?

How would you determine, from a fourier spectrum alone, if the source of light is homogeneous or not?

I have been reading about fourier transforms for the last few months, but I feel like I am still missing a huge point when it comes to connecting fourier transformations and image processing. It seems everyone is just using the amplitude information of an image while being in the fourier space, and ignoring phase.


texts:
I have read about spectral leakages in other posts here. From what i understand, it occurs because you dont have an integer number of time periods in your sampled data. By default a rectangular window (=no window) is used and it corresponds to a sinc function in the frequency domain. If we have an integer number of time periods in our data the zeros of the sinc cancel out every other frequency component giving us the correct value. If however we dont have an integer number of time periods in our sample, that corresponds to a shifted sinc and this is convolved with the frequency spectrum of the original resulting in leakage. In essence there is always leakage but it is masked by carefully choosing the number of samples to be an integer number of time-periods.

I wrote a simple program in scilab to see if i can remove this leakage with windowing. f(x) = Asin(20*pi*X). Maximum frequency = 10, sampling rate Fs = 30. I use a sinc window. A finite sinc window leaks in the frequency domain. For n samples the FFT(sinc) = (m1, m2, m3,...mn), m2 to mn are mirrored about the center. I changed the FFT(sinc) to FFTnew = (-m1, m2/abs(m2)-m2,-m3,-m4,...,mn/abs(mn)-mn) so that abs(FFT + FFTnew) = (0,1,0,...,1). I now create a new window called modifiedSincWindow = sinc + InvFFT(FFTnew). (in scilab invFFT(FFT(f(x)) = f(x), no scaling issues). the abs(FFT(modifiedSincWindow)) has all 0s except for 2 1s as expected. This is obvious from the linearity property of the DFT.

a) NumSamples = 30. FFT(f(x) x modifiedSincWindow) gives me the correct amplitude A.

b) NumSamples = 31. FFT(f(x) x modifiedSincWindow) does not give me the correct amplitude A. it doesnt even look anything like in (a). why is this? if we are convolving the frequency spectrum of the signal with the spectrum of the modifiedSincWindow and since the component falls in some bin shouldnt we get the same result as in (a). the spectrum for the modifiedSincWindow is exactly 1 for a particular bin and zero for the rest. And convolving this with the spectrum of f(x) should get me the same result as (a) but it isnt. Can anyone help me understand what exactly is happening here? why is it not working as expected? thanks for your comments.

EDIT: Added scilab code below.

My question basically boils down to if we are getting leakage because of the convolution by the sinc function, why not multiply the signal with a window whose magnitude response is 1 for one bin and zero for the rest. This should give us the right amplitude value irrespective of whether the number of samples is an integer number of time periods or not. right?

maxFreq = 10;
T = 3*maxFreq;

numSamples = T; %%(b) numSamples = T+1
timeRes = 1/T;

%%sampling the sinc function

sincVals = zeros(numSamples);

for i = 1:1:numSamples
    xval = (i - numSamples/2) * timeRes;
    if(xval == 0) %%then
        sincVals(i) = 1;
    else
        sincVals(i) = sinc(xval);
    end
end

fftVals = fft(sincVals);

modifiedFFTArray = zeros(numSamples);

for i = 1:1:numSamples
    modifiedFFTArray(i) = -fftVals(i);
end

modifiedFFTArray(2) = fftVals(2)/abs(fftVals(2)) - fftVals(2);
modifiedFFTArray(numSamples) = fftVals(numSamples)/abs(fftVals(numSamples)) - fftVals(numSamples);

invModified = ifft(modifiedFFTArray);

%%sincDiffs magnitude response is 0 for all bins except for 2 bins where it is 1

sincDiff = sincVals + (invModified);
signal = zeros(numSamples);

for i = 1:1:numSamples
    xval = i * timeRes;
    signal(i) = 1.5*cos(2*pi*maxFreq*xval) * (sincDiff(i));
end

binRes = T/numSamples;

%%disp(binRes)
binID = 1+floor(maxFreq/binRes); %%indexing is from 1 so add 1
%%disp(binID)
fftVals = fft(signal);

plot(abs(fftVals))


texts: How to do sliding window multiplication in frequency domain

texts:
I want to perform the following computations with the impulse response of length L of a LTI system:

Multiply the impulse response with a window w(n) of length M << L.
Take the FFT of the result in 1 and save as a vector Y2.
"Move" the window to the "right" by M/2 samples, i.e. take w(n-M/2), and multiply the impulse response with this new window.
Take the FFT of the result in 3 and save as a vector Y2.
Move the window to the right again by M/2 samples, i.e. take w(n-2*M/2), and multiply the impulse response with this new window.
Take the FFT of the result in 5 and save as a vector Y3.
Etc... (keep doing this process until the end of the impulse response, I hope you get the idea).

I think this is easy to do with a for loop in Matlab. However, I already have immediate access to the frequency response (the FFT of the impulse response). I was therefore wondering if there is a way to get my vectors Y1, Y2, Y3, etc... by computations in the frequency domain instead? I know that multiplication in the time domain corresponds to convolution in the frequency domain. So to get for instance the vector Y5 with the above notation, would it be possible to do as follows?
Take the FFT of the window w(n-5*M/2), and simply convolve it with the frequency response? And how do I know if it's faster/slower than just taking the inverse FFT of the frequency response and doing the steps described above?


texts: How do I derive complicated robotic motion models easily?

texts:
I have a filter that tracks a robot. I want it to use a 2D coordinated turn polar velocity motion model (from page 15 here):

But I want to expand on this motion model: I want an additional velocity state that is orthogonal to the heading, and I also want acceleration states for the two velocities.

How do I find the resulting discrete-time linearized model easily? I guess I could someone smarter than me could find the find the continuous time differential equation by hand (by thinking really hard?), and then discretize/linearlize (or linearlize/discretize) it. But that seems like an awful amount of work. I think I can manage the "discretize/linearlize"-step, but I don't know how to start on the "find the continuous time model" step.

It would be nice if there were some program that takes my desired states as input and outputs a model. But my googling turns up nothing. Is it really hard to find motion models? Do I fail to see the gravity of the work if I feel like this could be done in an afternoon? What should my workflow look like? Where do I start?


texts: Efficient Method for Interpolating between bins in FFTW

texts:
I'm working on some oscillator classes right now and perform FFT around 100 times per second. The issue I'm running into is interpolating between the bins so there is not a noticeable change from each execution of FFT. I'm fairly new to using FFTW right now and I'm still getting to know its functions.
For some more context, I'm working on some spectral morph functions and other fun like that. As of right now, I'm not using SIMD for my interpolation and it's just a basic for loop that writes to a fftw_complex pointer based on the two fft_complex bins I'm interpolating on. The idea is that I only do FFT when needed, at a max of roughly 100 times per second, and between every n number of frames or so, calculate the interpolated bin between the current one and the last one. My FFT length is currently 2048 samples because I've found that's the sweet spot for CPU and low harmonics. This makes the FFT bin a size of [2048][2].
I'm also only a few months into SIMD, and plan on using it for the final interpolation algorithm. Am I missing an FFTW function that can do a lot of this for me, or is there some other method out their that is faster?
Edit: My bad - When I mean interpolating the bins, I am referring to them being in the time domain, not frequency. I am storing my waveform in the bin in time domain so I don't need to do memory stuff when I want to do forward FFT.
Second Edit: I suppose I was a little off, but I had the right idea. I was working with another FFT library for a while that did not allow easy conversions from complex objects to arrays. I was under the presumption that I had to be restrained to only using FFT bins. This clears a lot up


texts: If MFCC vectors are speaker independent, then how can we reconstruct unique speech signals from them?

texts:
Reconstructing a speech signal from a collection of MFCC vectors seems to work pretty well, but I've heard that one advantage of MFCCs is speaker-independence i.e. they are more-or-less the same across different speakers for a given phoneme. How then can a speech signal, with all its speaker-dependent idiosyncrasies (accent, etc.) be reconstructed from a supposedly speaker-independent MFCC vector? Are MFCCs not in fact speaker-independent then? If not, what does determine speaker-independence vs. speaker dependence?

Thanks.


texts: Why does diagonal loading of a covariance matrix make an adaptive beamformer more robust in the case of a perturbed array?

texts:
It has been shown that 'diagonal loading' a covariance matrix derived for an adaptive beamformer can improve robustness of the beamformer when the antenna array is perturbed, albeit at the expense of background noise power.

The question is why does this work? Is the diagonal loading stopping the adaptive filter from reaching a local minima? Or is there some other mechanism at work here?


texts: What steps come after finding the gradients of the SIFT key points?

texts:
I am using a SIFT algorithm to extract features from an image.  I understand that the SIFT descriptor first finds the extrema and then finds the gradient and direction of each of these interest points.

What happens after this part?


texts: Are there frequency detection methods that are subject to less latency than the Goertzel algorithm?

texts:
I am trying to detect oscillations in a small range around 1.5Hz as early as possible. I am currently using the Goertzel algorithm with bins 1.0Hz, 1.2Hz,..., 1.8Hz, and 2.0Hz; with block size N = 500 and sampling rate Fs = 100 (the sampling rate, in my circumstance, cannot be modified).

I am applying the Goertzel algorithm on two signals, and checking for three conditions that, if met, result in the oscillation I'm concerned with being detected:

amplitude of signal 1 > threshold_1
amplitude of signal 2 > threshold_2
phase difference between signals > threshold_3

If all conditions are met, a flag is raised, and that is the purpose of this whole oscillation detection endeavor.

This all works fine, except that the Goertzel algorithm is inherently subject to latency in accordance with N. I am wondering if there are any other methods that do better than the Goertzel algorithm in terms of latency? Or if there's any way to reduce the latency of this implementation while still confidently detecting oscillations around 1.5Hz. I suppose I could reduce N but I'm not sure I want to be increasing the bin width.


texts: spectral subtraction for stationary noise

texts:
my recording was suffering from stationary background noise. I used spectral subtraction to remove it. the noise is reduced in loudness but not eliminated completely.
before:
http://vocaroo.com/i/s18FtoUwSah4
after:
http://vocaroo.com/i/s12U4feLs6iu

is spectral subtraction working?What can I do to further improve my results.


texts: Fractional spaced equalizer + timing (clock) recovery

texts:
I am trying to figure out how timing recovery is used after equalization (or, being more correct, when the equalizer is in the timing recovery loop), having a T/2 fractional spaced equalizer and using the Gardner algorithm as timing error detector.
The usual situation is to have the equalizer output at 1sample/symbol, by having downsampled, or combining the outputs of the even and odd filters. The latter is my case (case (a) in figure). But Gardner algorithm works at 2samples/symbol, so I wonder what I must do with the samples coming out from the even and odd filters. Must I multiplex the outputs of the two filters so I have a 2samples/symbol, as depicted in (b)?

Regards,


texts: DSP, Related Math, Statistics, and Jitter Primers

texts:
I'm a CS who took an internship at an EE shop, and I'm trying to do some work developing in house and client applications for a chip designer, and I'm having some fundamental communication problems. Mainly I shore up any applications they need developed from their specs, including at the moment a tool that demonstrates Jitter on a given PLL. I've had some high level stats, linear algebra, and DIP courses that I thought would have prepared me for this, but right now I've got a 50 pages on how to calculate Peak to Peak Jitter, Period Jitter, Random Jitter, and Long Term Jitter, and I can't seem to get a solid equation/algorithm out of any of it. It seems to approach RMS like it's a constant number.

I'm not sure what I need to tackle this. Does any one have a(some) primer(s) on the subject(s)?


texts: Why is being jointly WSS important in signal estimation with LMMSE estimator?

texts:
Assume that we have $y_k = s_k + n_k $.

We have observed $y_k$ and want to estimate $s_k$. The goal is to use LMMSE (Linear MMSE) estimator to find $y_k$ and on the other hand, we know that our filter is causal, it means we have

$$\hat{s_k} = \sum\limits_{m=-\infty}^k  h_{k-m}y_m\text,$$ but why should $s_k$ and $y_k$ be jointly WSS?


texts:
What are some applications of virtual microphones (as in here http://sail.usc.edu/aigaion2/index.php/attachments/single/119.)?
Are there good references (books or comprehensive tutorial papers) for this technology?


texts: BIBO Stability and the convergence of the frequency response of a system

texts:
It is my understanding that an LTI system is BIBO stable if and only if its impulse response $h(t)$ is absolutely integrable. This also happens to be one of the Dirichlet conditions for the convergence of the Fourier Transform, from which we can conclude that $H(j\omega)$, the frequency response of the system (i.e. the Fourier transform of the impulse response) converges if $h(t)$ is absolutely integrable.

Now think of an ideal low-pass filter. Its impulse response $h(t)$ is a $\mathrm{sinc}$ function which is not absolutely integrable. In other words, it is not BIBO stable. However, it does have a Fourier transform despite violating the condition of BIBO stability (or equivalently one of the Dirichlet conditions). Why is that? Is it because the Dirichlet conditions are a set of sufficient but not necessary conditions? Or is it because the Fourier transform of a $\mathrm{sinc}$ function doesn't "exist" in the usual sense unless we allow step functions?


texts: Envelope following using FFT bin power -> Low pass effect on the modulation transfer function

texts:
I am deriving the temporal envelope of a signal (in a particular frequency band) by summing some FFT bin powers and taking the average, I do the FFT frame by frame, windowing and some overlap.

Please help me understand the following: when the envelope of the signal changes slowly (say 5 Hz), I can track the envelope fine. However when the envelope is modulated at a higher frequency, say at 100 Hz, the dynamic range and amplitude of the detected envelope are very reduced.

So the Modulation Transfer Function (MTF) of this method has a low-pass look.

I am trying to understand the reasons for this low-pass roll-off effect and the maths behind it. I see what's happening but I can't clearly conceptualise the underlying cause . It should also be possible to come up with a number regarding the slope of the MTF.

Thanks in advance!

GP,

Edit 1 : Thanks for your answer. I can't post images cause I'm a new user apparently.. so here is a link: Link . you see on the left the input signal is 4 kHz modulated by 5 Hz and I take some fft bin powers to track the 5 Hz envelope and it works fine. But at 75 Hz modulation (right), the amplitude and dynamic range of the envelope are reduced, and it gets worse as you increase the modulation freq (low pass effect). This is what I am trying to get my head around: where is this effect coming from ? Is it because I destroy phase information when I take the power of the bins ? I apply a hann window on the time domain input frames, before computing the FFT.


texts: Is it right to talk about dynamic range on frequency domain?

texts:
I'm comparing three frequency response's plots from simulations and experiments and I observed that two of them look like clamped to a DC component, where the max dB value is i.e. 100 and the min is 40 dB, while on the third one the max is i.e. 120 dB and the min is 10 dB. I want to use the correct term to talk about these intervals.

I've used the term dynamic range to express an amplitude variation interval on time domain, but is it right to use it on frequency domain too?


texts: Are the HOG features of the images good enough to do image classification alone? Why are their variances low?

texts:
I've also asked this question here:https://cs.stackexchange.com/questions/77884/are-the-hog-features-of-the-images-good-enough-to-do-image-classification-alone

I've been working on classifying human poses like the pictures below. For that, I'm computing HOG features, and while in MATLAB, after using the command extractHOGfeatures, it's giving me 1764 HOG features, I also notice that maximum variance of all the 1764 HOG features are as low as 0.04 for a sample size of 2660, meaning the HOG features aren't probably changing much from sample-to-sample. My questions are:

1) why do you think we're getting low variances for HOG features? 2) what could we intuitively say about HOG features of these pose images? Are they good or interesting enough to classify images alone?


texts:
Feature detection is an essential task in low-level vision.

Good features are those that resist to different perturbations such as

noise addition,
blur,
geometric transforms (3D rotation with perspective, scaling),
radiometric transforms (monotonic or non-monotonic grayscale remapping).

By resistance, I mean that despite those alterations, the same features will continue to be detected.

My question is how can we mathematically express these types of resistance in general terms ? What would be the form of the equations, assuming grayscale images and features defined from the pixel values in a region of arbitrary shape ? (I am not asking for an assessment of the stability of known features described in the literature.)

Answering the question for a 1D signal would already be a good start.

Can anyone help ?


texts: What is the frequency response function for the Bode Plot pictured here?

texts: I've forgotten how to get the frequency response function out of a Bode Plot. The phase component of $H(j\omega) = e^{j\pi/2}$, right? But I am having trouble finding $|H(j\omega)|$. How would I get it from decibels?


texts: Correlation between a signal segment and its stretched replica

texts:
In my application the measurements are affected by temperature and the signal is stretched over time, though preserving relatively similar structure.

I want to find incremental stretching between pairs of consecutive signals. Having a pair of signals, I divide the first into N segments (windows) and for each window I want to find it's shift and new length in the next signal.

I tried to estimate these incremental shifts with peak detection, but got low accuracy estimates, probably due to the stretching effects on the waveform structure.

In addition, I was advised to try Dynamic Time Warping algorithm, but it manipulates both signals to find the minimal Euclidean distance. In my case I would want only the first segment to be warped until it correlates perfectly with some waveform in the 2nd signal. Couldn't find the way to do so with DTW.


texts: Why does knowing the impulse response allow you to determine the output for any LTI system?

texts: Going over some notes on LTI systems and I’m quite confused about how knowing the impulse response of a system allows us to recover the systems’ response to any signal. In the textbook screenshot below, it looks like they just replaced the Dirac deltas with $f$ when defining $y(n)$.


texts: Retrieve pixel values from the output jpeg image

texts:
I have this jpeg image (abc.jpg) and I change some pixel values of this image and store this on hard disk, to later see the changed the pixel values. The name of the output jpeg image is say, xyz.jpg. The problem here is in the image xyz.jpg. The changed pixels does not give me the exact values which I changed them to.

Example:

In abc.jpg at the pixel position (x,y)=(1,10) the pixel value is (231, 120, 120). Now, I want to change it to (230, 120, 120). Now, I save this output image as xyz.jpg. When I try to retrieve the values from position (1,10) I get the pixel value as (232, 120, 120). Now you see the problem, I can't get the value to which I have changed. I am doing all this in MATLAB. And yes, I have tried to save the output image as .tiff, .bmp, .gif. But the size of these formats is very large as compared to .jpg. Please suggest some solution to this problem.


texts: Fourier Transform of a PSD and response of a PSD input

texts:
Does a Fourier Transform of a white noise exist?
If so, what is its general form?

It is possible to compute the fourier transform of a Power Spectral Density?

My problem in detail. I have to compute the response of a dynamic system given the PSD of an input ($ \Phi_{vv} $).

I thought to solve the problem in three step:

Build a shape filter (system which, given a white noise input ($ n $), gives a colored noise as output)
Express both the shape filter and the system dynamic with a fourier transform
If i can manage to express also the input (the white noise) in the fourier domain I should have solved the problem


texts: How to find out if a transform matrix is separable?

texts:
In image processing, when we have a series of basis images, how could we know if the transformation is separable or not?
For example, I know that following bases are separable and transformation can be done with matrix $A$, but I don't know how to find $A$ from $A_1$, $A_2$, $A_3$ or find separability.
\begin{equation}
A_1 = \frac{1}{2}\begin{bmatrix}
    1 & 1 \\
    1 & 1
\end{bmatrix}
A_2 = \frac{1}{2}\begin{bmatrix}
    1 & -1 \\
    1 & -1
\end{bmatrix}
A_3 = \frac{1}{2}\begin{bmatrix}
    1 & -1 \\
    -1 & 1
\end{bmatrix} \\
A = \frac{1}{\sqrt 2}\begin{bmatrix}
    1 & 1 \\
    1 & -1
\end{bmatrix}
\end{equation}


texts: What characterizies 'causality' for a finite FFT?

texts:
Causality of a LTI transfer function $G(\tau)$ in the continuous time domain, i.e. for
$$y(t)=\int G(\tau)x(t-\tau)d\tau$$
is characterized by
$$G(\tau < 0) = 0$$
By the way, in the frequency domain, causality is characterized by the Kramers-Kronig relations, which is why I am particularly interested in the Fourier transform $G(\omega)$.

Now I am analyzing the functions $x(t)$ and $y(t)$ in discrete time and for a finite set of equidistant samples (measurements). In order to work in the Fourier domain, I will have to assume that they are periodic functions (or rather series):
$$x(t+T)=x(t) \qquad y(t+T)=y(t) \qquad G(t+T)=G(t)$$
Even with such techniques as zero-padding, they are still periodic as seen from the view of FFT.

But what means causality for periodic functions? If $G(t)$ is the response of a periodic Dirac comb, it is not clear if it is the causal response of the Dirac impulse at $t=0$ or the non-causal 'response' at some of the later Dirac impulses at $t=nT$ where $T$ is the measurement duration. Strictly speaking, it is the response to all the Dirac impulses in the comb and it makes no sense to define the notion of 'before or after' the impulse.

Am I missing something? Is there a proper notion of causality for periodic functions and hence on finite domains?


texts: What is meant by 20 dB signal-to-noise ratio?

texts:
I'm reading a paper in which a discrete signal

$$x(n) = s(n) + w(n)$$

is considered. $s(n)$ is a known, deterministic series, and $w(n)$ is white noise with zero mean. The authors write that

  signals were generated with an SNR of 20 dB

What does this mean?

What does signal energy mean? There seems to be several ways of defining this, but no attempt to do so is done in this paper.
What does 20 dB SNR mean? I suppose that either $20=10\log{E_s/E_w}$, or $20=20\log{E_s/E_w}$, but am not familiar enough with these concepts to be sure.


texts: Noise and harmonic elements w/in a discrete time signal

texts:
So I have a discrete-time signal with around 300 samples. It is part of an environmental science application, so is not-strictly periodic although there clearly are periodic components. I have applied an FFT (using Python NumPy.fft) to calculate power spectral density plots and also then applied an IFFT to try and reconstruct the periodic components of the signal (shown in the attached plots).
I am looking to have a measure for what % of the signal is periodic, and what % is noise/residual. In addition, I would ideally be able to extract the equations harmonic functions with high power. I figured the IFFT could be compared to the original signal via a Person's correlation to find the share of the variance in the time signal explained by the periodic components captures by the FFT (using R^2). This seemed to make sense to me, but I have no idea how to proceed to find the equations for the harmonic components.
Is anyone familiar with a method that does this, or does anyone see some glaring issue with my IFFT-time signal correlation comparison?
Thanks!


texts: At what amplitude does spectral leakage start to be noticeable in the phase spectrum?

texts:
I understand that discontinuities at the boundaries of a finite signal may introduce spectral leakage, that alter both the amplitude and phase spectra of the FFT, specially at low amplitudes. That is the explanation to several questions like (baffled by fft phase spectrum!"). In such questions, the examples are normally based on sine waves, where there is a dominant frequency, and the spectral leakage is noticeable in all the other frequencies.

Signals are often more complex, with a more continuous amplitude spectrum, but even in that case - specially if they are low pass-filtered, there is a point where there is a noticeable discontinuity in the amplitude and phase diagrams.

Example (in Octave):

N = 101;
x = exp(-linspace(-5,5,N).^2);
f = fftshift(fft(x));
subplot(211)
semilogy(abs(f))
subplot(212)
plot(unwrap(angle(f)))

Here the amplitude spectrum is bell-shaped as I expected at low frequencies. But when the amplitude goes below approximately 1e-11, the amplitude spectrum flattens.

At the frequencies where the amplitude spectrum is bell-shaped, the phase spectrum is a straight line with slope $\pi (N-1)/N$, which coincides with the effect of a time shift of $(N-1)/2$ from an even signal. At higher frequencies, it goes more or less linearly towards $\pi/2$ (although in this case the slope is not strictly constant).

Why there is such a harsh change in the spectra at that frequency? Is it a matter of numeric error in the computation of the FFT? The amplitude at the discontinuity is small (about 1e-11), but greater than machine error. And for other signals with narrower spectra, the amplitude at the discontinuity is higher. For instance, the same example but using a "wider" time signal:

x = exp(-0.5*linspace(-5,5*(N-1)/N).^2);


texts: Noise sensitivity of the (classical) Empirical Mode Decomposition routine

texts:
I tried to apply a MATLAB Empirical Mode Decomposition routine to denoise a signal, basically retaining only the last IMFs, with a criterion based on the mode energy.

To validate the routine, I have built a synthetic signal with added Gaussian noise + a sinusoidal disturbance. I noticed that the EMD routine (at least mine) seems very sensitive to the noise. In fact, if I launch it twice, the generated noise is of course different and the IMFs are also quite different.

Do you have any suggestions on how to "stabilize" the routine?