Deconvolving a 1d Signal Using a Lookup Table of Kernels assuming I measure a signal that has different PSFs per position in time.
for example:
t = linspace(0,20);

% "ground truth" signal will be like:

x = @(t0) exp(-(t-t0).^2/0.1)  ;

% some made up impulse response (psf) that depends on t0 will be:

h = @(t0) diff(diff(  exp(-(fftshift(t)).^2./(0.1*t0) )));

% the convovled signal:
y = @(t0) conv(x(t0),h(t0) ,'same');

% now if I have signal from two positions, I get:
plot(t,y(3)+y(15))

Note that the two peaks now are distorted differently as function of their position.
What methods can I use here,  given that I have such a PSF lookup table, such as that h = @(t0)... above, to deconvolve my 1D signal, even though it will behave differently in different positions as seen in the plot? just a standard deconvolution wont work here.
EDIT : trying to clarify the question further. I am looking after a way to "deconvolve" the signal that is distorted by such position dependent PSF. So instead of these two features I will be able to trace back the original signal (this case just two peaks).  Using standard de-convolution schemes will not work well because they assume an effective single PSF, and here we have a "family" of PSFs. Is there a way to solve it?
I was hoping, for example, that extending the dimension of the PSF will allow to accommodate such effect, or maybe using other tools to "train" a system to understand it.
EDIT 2:
Here's a file that show an example of x - the ground truth signal, y - the signal convoluted by the positions dependent psfs (or kernels) and, psfs -  an array of kernels per position.

The way I understand the problem is each sample of the output is a linear combination of the samples of the input.
Hence it is modeled by:
$$ \boldsymbol{y} = H \boldsymbol{x} $$
Where the $ i $ -th row of $ H $ is basically the instantaneous kernel of the $ i $ -th sample of $ \boldsymbol{y} $.
The problem above is highly ill poised.
In the classic convolution case we know the operator matrix, $ H $, has a special form (Excluding the borders) - Circulant Matrix. With some other assumptions (Priors) one could solve this ill poised problem to some degree.
Even in the case of Spatially Variant Kernels in Image Processing, usually, some form is assumed (Usually being block circulant matrix, and the number of samples of each kernel is larger than the number of samples in the support of the kernel).
Unless you add some assumptions and knowledge into your model the solution will be Garbage In & Garbage Out:
numInputSamples = 12;
numOutputSamples = 10;

mH = rand(numOutputSamples, numInputSamples);
mH = mH ./ sum(mH, 2); %<! Assuming LPF with no DC change

vX = randn(numInputSamples, 1);

vY = mH * vX;

mHEst = vY / vX;

See the above code. You will always have a perfect solution yet it will have nothing to do with mH.
Now, if I get it right, you say I don't know $ H $ perfectly, but what I have is a pre defined options.
So let's say we have a matrix $ P \in \mathbb{R}^{k \times n} $ which in each row has a pre defined combination:
$$ H = R P $$
Where $ R $ is basically a row selector matrix, namely it has a single element with value $ 1 $ in each row and the rest is zero.
Something like:
mP = [1, 2, 3; 4, 5, 6; 7, 8, 9];
mH = [1, 2, 3; 7, 8, 9; 7, 8, 9; 4, 5, 6; 4, 5, 6];

% mH = mR * mP;

mR = mH / mP;

So our model is:
$$\begin{aligned}
\arg \min_{R, \boldsymbol{x}} \quad & \frac{1}{2} {\left\| R P \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} \\
\text{subject to} \quad & R \boldsymbol{1} = \boldsymbol{1} \\
& {R}_{i, j} \geq 0 \quad \forall i, j \\
\end{aligned}$$
It is still exceptionally hard (Non convex) problem but with some more knowledge it can be solved by utilizing alternating methods where we break the optimization problem as:

Set $ \hat{\boldsymbol{x}}^{0} $.
Solve $ \hat{R}^{k + 1} = \arg \min_{R} \frac{1}{2} {\left\| R P \hat{\boldsymbol{x}}^{k} - \boldsymbol{y} \right\|}_{2}^{2} $ subject to $ R \boldsymbol{1} = \boldsymbol{1}, \; {R}_{i, j} \geq 0 \; \forall i, j $.
Solve $ \hat{\boldsymbol{x}}^{k + 1} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| \hat{R}^{k + 1} P x - \boldsymbol{y} \right\|}_{2}^{2} $.
Check for convergence, if not go to (2).

Now each sub problem is convex and easy to solve.
Yet I still recommend you add better assumptions / priors.
Such as the minimal number of contiguous samples that must have the same PSF (Similar to 2D in images where we say each smooth area is smoothed by a single PSF).
Remark
We didn't employ the fact each element in $ R $ is either 0 or 1 as the straight forward use of that will create a Non Convex sub problem.
In case the number of PSF's is small we can use MIP solvers. But the model above assumed each row is a PSF so for large number of samples even in case we have small number of PSF'w the matrix is actually built by shifting each PSF as well. So we'll have large number in any case.
Another trick might be something like Solving Unconstrained 0-1 Polynomial Programs Through Quadratic Convex Reformulation.
Yet the simplest method would be "projecting" $ R $ into the space (Which is not convex, hence projection is not well defined). One method might be setting the largest value to 1 and zero the rest.
Update
In comments you made it clear you know the kernel per output sample.
Hence the model is simpler:
$$ \boldsymbol{y} = A \boldsymbol{x} + \boldsymbol{n} $$
The least squares solution is simply $ \boldsymbol{x} = {H}^{-1} \boldsymbol{y} $.
For better conditioning and noise regularization (Actually prior about the data, but that's for another day) you can solve:
$$ \hat{\boldsymbol{x}} = {\left( {A}^{T} A + \lambda I \right)}^{-1} {A}^{T} \boldsymbol{y} $$
This is a MATLAB code for proof of concept:
load('psfs.mat');

mA = psfs;
vY = y;
vX = x;

vParamLambda    = [1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1];
numParams       = length(vParamLambda);

vValMse = zeros(numParams, 1);

mAA = mA.' * mA;
vAy = mA.' * vY;
mI  = eye(size(mA));

for ii = 1:numParams
    paramLambda = vParamLambda(ii);
    vEstX = (mAA + paramLambda * mI) \ vAy;

    vValMse(ii) = mean((vEstX(:) - vX(:)) .^ 2);
end

figure();
hL = plot(vParamLambda, 10 * log10(vValMse));
xlabel('Value of \lambda');
ylabel('MSE [dB]');

This is the result:


Image processing: How to find center of biggest blob in an image? Let's say I have an image like this one, and have masked it so all yellow pixels are 1 and all black pixels are 0. I would like to (1) fit the largest possible circles containing solely yellow pixels inside the yellow blobs; multiple circles may fit in irregular blobs, (2) from all such circles, identify which has the biggest radius, and (3) find the location of its center. What image processing techniques can I use for this? I am familiar with some morphological operations (like these), but can't get to where I want to go yet.

do a distance transform. you'll see why that's a good idea: for every pixel you get the shortest distance to a border. that's exactly the radius of an inscribed circle.

from this, just find the pixel with the largest value.

if you're curious, throw a "non-maximum suppression" on it. that is a kind of "morphological" kernel operation where you set a pixel to 0 if any neighbor has a strictly larger value.


Q number format vs IEEE 754 format? I'm reading up on DSP implementations and Q number format for representing decimal and decimal fractional numbers. Why would you use the Q format and not a standard like IEEE 754? When would you use one instead of the other? They both encode numbers into a sequence of bits, right?

Q format numbers are fixed-point, which means they can be manipulated by integer ALUs, rather than needing to use a floating-point unit.  In a DSP setting, this is useful for greater speed and lower power consumption, since fixed-point arithmetic is much simpler than floating-point.  On the other hand, the main advantage of using a floating-point representation is greater dynamic range.


What is the acceptable range for bandwidth utilization factor in OFDM systems? The bandwidth utilization factor is defined by (alpha=bite rate/bandwidth), for OFDM system, what is the acceptable range for alpha?

It depends how good your system's anti-aliasing filter is. Your transmitter/receiver will have an antialiasing filter after the DAC (transmitter) or before the ADC (receiver). This filter will be designed to pass the band of frequencies that contain the signal of interest, while cutting off any content at higher frequencies that would alias into your signal's passband.

The higher the bandwidth utilization factor, the closer the passband edge is to the Nyquist rate, and thus the narrower your filter's transition region. A narrower transition region requires a more complicated filter (i.e. it will have a higher order, which can require more discrete components). Since such a filter would be implemented in the analog domain, this causes an increase in:

Design complexity
Physical size
Cost

OFDM is often used in widely-adopted wireless standards. In these cases, the final bullet above is a prime concern. Take WiFi as an example. WiFi transceivers are made by the zillions; even a small cost savings that can be realized by using simpler analog filtering can add up to significant money when multiplied by the numbers of units that mainstream electronic devices are sold in nowadays.


Have you ever seen logarithmic burn-in on a Kalman filter before? I'm reviewing a paper that uses a multiplicative factor $\ln i$ where $i$ is the number of steps since the Kalman filter was initialized.  The idea is to slowly build confidence as the filter sees more data.

The authors replace the typical Kalman gain $K$ with one weighted by this multiplicative factor $K_{new}$:

  $K_{new} = ( \ln i) K$

Have you ever seen such an approach used before? I'd like to give the authors a reference to cite to tie them into the literature on this subject if it exists.

The authors cite this work in their response to the review:

  Zhao H, Lu L. Adaptive recursive algorithm with
  logarithmic transformation for nonlinear system identification in
  α-stable noise [J]. Digital Signal Processing, 2015:S1051200415002535.

They take their inspiration from equation (14):

  $J_p(n)=\sum_{i=1}^n \lambda^{n-1} \cdot \log^p (1+|e(i)|)$

It is notable that $e(i)$ is (as one might expect) the error in the signal, and not linearly increasing with the number of steps since the filter was initialized.


Should the order of a filter be proportional to the sampling rate? I have designed a fir filter with my desired frequency response at a specific sampling rate and filter order. If I were to redesign my filter for a new higher sampling rate, should I change the filter order? I am having trouble understanding how the same filter order influences the desired frequency response at different sampling rates.

Keeping the FIR filter the same order and doubling the sample rate will effectively halve the accuracy of frequency resolution which will be noticable at low frequencies (this is based on the assumption you are using something like FFT to create your FIR, otherwise the order of the FIR and it's resolution is a bit more complicated).


Double Integrating Gaussian Noise My question has to do with integrating gaussian noise.

Let us assume we have samples of discrete gaussian white noise  with mean $\mu = 0$ and variance $\sigma_{th}^2$. These noise samples are passed through the system shown in the Figure (a cascade of two integrators with outputs $y_1[n]$ and $y_2[n]$, respectively).

What will be the mean and variance of $y_1$ and $y_2$ (let us say after $N$ cycles)?

Thanks to Dilip, I managed to figure out the answer to the question.
$y_1[n]$ will  be a normal random variable, with its variance given by the following equation:

\begin{equation}
\sigma_{y_1}^2=N\sigma_{th}^2
\end{equation}

This can be explained by using the fact that the variances will just add up. The mean will stay zero.$\\$
A similar approach can be used to find the variance of $y_2[n]$ ($\sigma_{y_2}$), which begins by writing the equation for $y_2[N-1]$ ($y_2[n]$ after $N$ cycles):

\begin{equation}
y_2[N-1]=(N-2)x[0] + (N-3)x[1] + (N-4)x[2]...+x[N-3]
\end{equation}

Taking variances on both sides:

\begin{equation}
\sigma_{y_2}^2=\sigma_{th}^2(1+2^2+3^2+...(N-2)^2)=\sigma_{th}^2\frac{(N-2)(N-1)(2N-3)}{6}\approx \sigma_{th}^2\frac{N^3}{3} (N>>1)
\end{equation}

Again, the mean will stay zero.


How to apply low-pass filter to an I/Q signal to get baseband signal? Can anyone help me with the demodulation technique of a modulated signal? I mixed it with I/Q signal. But I am very confused about low-pass filtering. I am trying to apply a low-pass filter, but I am not getting the baseband signal. What kind of filter can I use to get demodulated baseband signal?

My carrier frequency is $f_c=1 \ \mathrm{GHz}$ and my sampling frequency is $f_s=10 \ \mathrm{GHz}$.

You did not specify the bandwidth of your modulated signal, but for discussion let us assume it is 50 MHz. (Meaning all the signal of interest is in 50 MHz centered at your carrier at 1 GHz).

You may be downconverting to a lower frequency IF or doing a direct conversion to baseband in which case you will need to do quadrature downconversion to not have image issues. I will assume this is what you are doing which means you will be multiplying your real signal at 1 GHz with a sine and cosine using a local oscillator at 1 GHz (all digitally since you are in a sampled system, so an NCO).
After you downconvert with this multiplication each I and Q will have your signal of interest at baseband plus an image at 2GHz (the sum of the two frequencies at the multiplier inputs and this is what you want to filter. ("I" representing the real component of your complex baseband signal, and Q representing the imaginary component of your baseband signal). You need to eliminate the image at 2 GHz, but given your signal is only 50 MHz wide (+/-25 MHz at baseband), you would ultimately filter down to that bandwidth to eliminate interference and reduce noise, but to "see" your signal you could simply use a moving average to implement a very simple low pass filter given this significant separation between your baseband signal and the image.  A moving average filter would for a Sinc function in frequency with its first null at 1/T where T is the length of your average in time. So for example, if you wanted to place the first null at 500 MHz, you would average over 20 samples:

$T=\frac{1}{500e6}= 2$ ns. With a sampling rate of 10GHz, this is equivalent to 20 samples. What is convenient about this choice is that it will also have a null directly at 2 GHz which is the image you are trying to filter out!  So simply performing a moving average over 20 samples at the output of your mixer will significantly supress your image.  See the frequency response below where 1 on the frequency axis is half your sampling rate or 5 GHz.  In the example I gave with 50 MHz BW, the passband droop is negligible and the image supression (at .4 of the plot) would likely exceed 30 dB (I am eyeballing that) over your BW.  Point is that is it a quick filter to confirm your downconversion without having to get into an explanation on more detailed filter designs (which you would ultimately likely want).

Here are more details on the down-conversion process in case it helps:
To demodulate this signal directly to baseband, you would multiply it by $e^{-j\pi f t}$ where f = 1 GHz.  This is implemented with a sine and cosine at 1 GHz, and for best noise performance if you were concerned with that, you would do a full complex down-conversion as : $$(I_1+jQ_1)(I_2-jQ_2) = (I_1I_2+Q_1Q_2)+j(I_2Q_1-I_1Q_2)$$
Where

$(I_1+jQ_1)$ represents your modulated signal at 1GHz and

$(I_2+jQ_2)$ is your local oscillator for performing the down-conversion.

Since your 1 GHz signal is likely a real signal, to implement the lowest noise approach I suggest, you would need to split in quadrature before sampling, or alternatively implement the Hilbert transform digitally.  This also would not require a filter (although recommended of course to eliminate noise and interference out of band) as this is a full complex multiplier and therefore provides image rejection. Otherwise if you are ok with the 3 dB noise penalty and can filter out the image, $(I_1+jQ_1)$ would be real so $Q_1 = 0$, and the downconversion simplifies to

$$(I_1)(I2−jQ2)=(I_1I_2)-j(I_1Q_2)$$

Once you do that multiplication, you should see your signal of interest centered on baseband as $$I_{out} + jQ_{out}$$

Where
$$I_{out} = I_1I_2$$
and
$$Q_{out} = -I_1Q_2$$

Describing the I output ($I_{out}$) (Q would be similar) to see the image that would need to be filtered:

$I_1= cos(\omega_ct+\phi)$ where $\omega_c$ is the carrier frequency of the received signal with some arbitrary phase offset $\phi$ .

$I_2 = cos(\omega_ct)$ is the local oscillator I signal, tuned to the carrier frequency.

Using the cosine product relationship: $2cos(\alpha)cos(\beta)= cos(\alpha-\beta)+cos(\alpha+\beta)$

$I_{out}= I_1I_2 = cos(\omega_ct+\phi)cos(\omega_ct)= \frac{1}{2}(cos(\phi)+cos(2\omega_ct+\phi))$

Showing the high frequency component to be filtered out (in your case 2 GHz) and the signal at baseband with the phase offset as received.

Following the same process for Q results in

$Q_{out} = -I_1Q_2 = \frac{1}{2}(sin(\phi)+ sin(2\omega+\phi))$


Why more poles than zeroes? I read that an "improper system" "has more zeros than poles; it is not causal, cannot be implemented, has a strictly proper inverse and has infinite high-frequency gain."

Does causality fail due to instant signal change? I try to make up a recurrence, where $x_n$ depends on both $x_{n-1}$ as well as $x_{n+1}$ but they all can be interpreted as $x_{n+1}$ dependence on $x_{n}$ and $x_{n-1}$ and, thus, are causal.

Is it related with solving recurrent relations saying that P(x) in generating function $f(x) = f_0 + f_1 x + f_2 x^2 + \cdots = P(x)/R(x) \ \ \ \ \ $ must have degree less than R(x) as Miguel A. Lerma says in Generating function of Linear Homogenous Recurrence. Here, $R(x) = 1 + r_1 x^1 + r_2 x^2 + \cdots$ is a reciprocal of the characteristic polynomial of our recurrence $x_n + c_1 x_{n+1} + \cdots = 0.$ This indeed seems to happen whenever which R(x) I try.

On the other hand, https://ccrma.stanford.edu/~jos/fp/Existence_Z_Transform.html seem to identify causality with pole values (if your pole > 1 then series does not converge, which means anticausality) rather than their number and others say that causality is $x_n = 0$ for all $n < 0.$ I do not see how this is related with poles/zeros ratio and my guess that poles > zeroes = causality may be wrong. I have asked to relate various kinds of causality here.

I think the system is always non-causal if it has more zeros than poles. @Hilmar, your implementation might be possible digitally but if you see a transfer function like this:
$H(z) = \frac{(z-1)(z-2)}{(z-3)}$, taking the inverse Z-Transform,
$y(n-1)-3y(n-2) = x(n)-3x(n-1)+2x(n-2)$ is clearly non-causal and cannot be implemented real-time(of course you can store the signals and do processing offline, but real time implementation is impossible). So I guess Non-Causality is the answer, apart from the infinite high-frequency gain.

Update by Val
Ok, I understand now. In wikipedia, we see that $$a_0 y_n + a_1 y_{n-1} + a_2 y_{n-2} + \cdots = b_0 + b_1 x_{n-1} + b_2 x_{n-2} + \cdots$$ is identical to $$H(z) = {b_0 + b_1 z^{-1} + b_2 z^{-2} + \cdots \over a_0 z^{0} + a_1 z^{-1} +\cdots}.$$ I have finally noted that there are no positive powers. That is, $a_0$ corresponds to the current output, $z^0 = y[n]$, and $\cdots = a_{-2} = a_{-1} = \mathbf{0} = b_{-1} = b_{-2} = \cdots$. If we had more zeroes then we have more power in the nominator and, after normalizing $H(z)$ so that $a_0$ has $z^0$, we still have a positive power of z in the nominator, which means that output depends on future input.

In the simplest case, $Y(z) = zX(z)-x_0$ corresponds to $y[i] = x[i+1]$, that is output depends on future input (we can also drop $x_0$ assuming our sequences are causal, i.e. $0 = y_{-1} = x_0$ aka $z^{-1}Y = X$ or $Y = zX$). It is curious that causal is an example of non-causal. Could anybody clarify it further?

It also worth noting that $(1-az^{-1})(1-bz^{-1})$ is not a zero of order two. It has equal number of zeroes and poles, since it is equal to $(z-a)(z-b)\over z^2.$


First Differentatior and Integrator For discrete time series $Y$, first differentiation ($D_i=Y_i-Y_{i+1}$) and integrator ($S_i=Y_i+Y_{i+1}$) can be defined as two highpass and lowpass LTI digital filters. Where, transfer function for the first difference is:
$$
H(z)=1-z^{-1}
$$
and for the first integrator is:
$$
H(z)=1+z^{-1}
$$
These two transfer function corresponds to $\text{Conv}(Y,[1,-1])$ and $\text{Conv}(Y,[1,1])$ for differentiator and integrator, respectively.

I was wondering what the power transfer function of these two filters should be.

I have to add, in terms of implementation in Matlab, I can get $S_i$ and $D_i$ of a white noise by:

Y = randn(1,100);
D = diff(Y);
S = Y(1:end-1)+Y(2:end);

By filter in Matlab:

S_Filt = filter([1  1],1,Y);
D_Filt = filter([1 -1],1,Y);

Or by convolving conv the corresponding windows:

S_Conv = conv(Y,[1  1]);
D_Conv  = conv(Y,[1 -1]);

Note, all D* and S* variables should be identical (regardless of head and tail of time series, of course).

Thanks in advance.

I think I found an answer to this question. So, I will post it here.

For power transfer function of $1-z^{-1}$ we know:
\begin{eqnarray*}
1-z^{-1}&=&1-e^{-jw}\\
&=& 1-\cos(\omega)+j\sin(w)
\end{eqnarray*}
Then by using complex conjugate properties we have:
\begin{eqnarray*}
|1-z^{-1}|^2&=&1-2\cos(\omega)+\cos(\omega)^2+\sin(\omega)^2\\
&=&2(1-\cos(\omega)).
\end{eqnarray*}
And simillarly, for $1+z^{-1}$, we have:
\begin{eqnarray*}
|1+z^{-1}|^2&=& 1+2\cos(\omega)+\cos(\omega)^2+\sin(\omega)^2\\
&=& 2(1+\cos(\omega)).
\end{eqnarray*}

-

Note that the integrator in this question is a special case (see the expression for $S_i$) and is different than the conventional integrator $y[n]=y[n-1]+x[n]$.


How do I deepen my voice (I'm a female) to sound like a man in a video for a school project? Some extra info. that I think might help or make things clearer:

The video has several other people in it also talking (I don't want
to change their voices).
I have some experience in editing videos(nothing major I've only ever used
Flimora and Imovie)but I've never had to edit someone's voice.

I know some audio software that does a good job:

Audacity (free software) - effect -> change pitch (set the pitch to down, lower down the semi-tones and the percent change)
Adobe Audition (paid software) - there is a plugin called pitch shifter (lower down a little bit the semi-tones and cents)

Both are audio software which means you will need to export your audio from the video clip and because there are other people talking in the video you will need to select the parts only were you are talking and then apply the effect.


Autocorrelation matrix derivation Hi I am trying to derive the autocorrelation matrix and I am unsure about how exactly it works. I can see that the $4\times 1$ matrices result in the Hermitiain and Toeplitz matrix? Surely the only non-zero values in each of these is at $n = 0$ which is the variance?

Any help would be greatly appreciated!

You get a diagonal matrix because the values of $\eta(n)$ at different times (sample indices $n$) are uncorrelated, as shown in your first equation. This means that all off-diagonal elements of the auto-correlation matrix must be zero (they equal $E[\eta(n)\eta(n+m]$ for $m\neq 0$). Only the diagonal elements are non-zero, and they are equal to the variance.


Finding where signal is showing some (unknown) periodicity I have readings of the vertical displacement of someone over time, as in the picture below:

At some points in time, the person will be exercising or jumping up & down, creating more or less obvious periodic motions (yellow regions of the graph). I am now trying to find a way to automatically spot these yellow areas, by looking for periodicity in the recorded signal.

My initial attempt has been to use a spectrogram, using the specgram built-in function in Octave:
segmentLength=150;
window=hanning(150);
overlap = 80% * segmentLength;
specgram(x, segmentLength, 1, window, overlap)

Which results in the following output:

This seems to be hinting in an ok direction since I see for example that something gets picked up just before frame 2000, however it is not very conclusive. Given my limited understanding of the spectrogram, I am not sure how to proceed from here.

Is the spectrogram the right approach? If yes, how can I improve the output? If no, what approach would you recommend for this task? I am aware that some of the windows that I am trying to identify will be too difficult, but was expecting to at least pick up windows 1, 2 and 5 (starting from the left). Thanks for any pointers!

EDIT

I have tried to implement a simple ASDF in octave with the following code:

#Parameters for asdf analysis
  N=400; %Window size
  kmin= 12;
  kmax= 40;
  step = 50;

  n0min = floor((N+kmax)/2);
  n0max= size(x,2) - floor((N+kmax)/2) -1;
  Q=zeros(kmax-kmin+1, floor((n0max-n0min)/step)+2);

  i=1;
    for k= kmin : kmax %i

    j=1;
    for n0= n0min : step : n0max %j
      for n= 1 : N-1
        Q(i, j) = Q(i,j) + ( (x(n+n0-n0min) - x(n+n0-n0min+k))^2 * hanning(N)(n+1) );
      end
      Q(i, j) = (2/N) * Q(i,j);
      j++;
    end

    fprintf('Calculated %d out of %d periods \n', k-kmin+1, kmax-kmin+1);
    fflush(stdout);

    i++;
  end

  imagesc(Q)

Tweaking around with the parameters, I ended up with the ones shown in the code, that yield the following result:

This is already getting closer to what I expected. Other suggestions welcome :)

If you're gonna use ASDF to measure periodicity and you want to window the data with something other than rectangular window, do it after subtracting and squaring:

$$ Q_x[k, n_0] \triangleq \frac{2}{N} \sum\limits_{n=0}^{N-1} \left(x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor] \ - \ x[n+n_0-\left\lfloor \tfrac{N+k}{2}\right\rfloor + k] \right)^2 w\left(\tfrac{n}{N}\right) $$

where

$\left\lfloor \cdot \right\rfloor$ is the floor() function and, if $k$ is even then $ \left\lfloor \frac{k}{2}\right\rfloor = \left\lfloor \frac{k+1}{2}\right\rfloor = \frac{k}{2} $

and $w\left(\tfrac{n}{N}\right)$ is a window function of non-zero width $N$ samples centered at $n=\tfrac{N}{2}$.  A Hann window would be

$$ w(u) \triangleq \begin{cases}
\tfrac12 - \tfrac12 \cos(2\pi u) \qquad & 0 \le u < 1 \\
\\
0 & \text{otherwise} \\
\end{cases}$$

To make this ASDF into "autocorrelation" (in the neighborhood of sample $x[n_0]$) defined from the ASDF:

$$ R_x[k,n_0] = R_x[0,n_0] - \tfrac12 Q_x[k, n_0]  $$

where

$$ R_x[0, n_0] \triangleq \frac{2}{N} \sum\limits_{n=0}^{N-1} \Big(x[n+n_0-\left\lfloor \tfrac{N}{2}\right\rfloor]\Big)^2 w\left(\tfrac{n}{N}\right) $$

Since $Q_x[0, n_0] = 0$ and $Q_x[k, n_0] \ge 0$ for all lags $k$, that means that $ R_x[k, n_0] \le R_x[0, n_0] $ for all lags $k$.

Suppose for a minute that $x[n]$ is periodic with period $P$ (and $P$ happens to be an integer), then

$$ x[n+P] = x[n] \quad \forall n $$

and  $Q_x[mP, n_0] = 0$ and $R_x[mP, n_0] = R_x[0, n_0] \ge R_x[k, n_0]$ for any integer number of periods ($m$ is an integer).  So you get a peak at $k=0$ and at $k$ equal to any other multiple of $P$ if $x[n]$ is periodic.  If $x[n]$ is not perfectly periodic, what we might expect is the biggest peak at $k=0$, another peak (but slightly smaller) at $k=P$ (the period we are looking for) and progressively smaller peaks for larger multiples of $P$.

The measure of periodicity in the neighborhood of $n_0$ would be

$$\frac{R_x[P, n_0]}{R_x[0, n_0]}$$

Hmmm, I just thought of something.  Because of my audio/music-centric POV, I have been assuming no DC bias (i.e. $x[n]$ should maybe come out of a DC blocking high-pass filter).  But I don't think that is the case with your application.  So maybe just stick with ASDF, but somewhere you will need a threshold for how low $Q_x[k, n_0]$ can go (for $k \ne 0$) for you to consider that $x[n]$ is periodic in the neighborhood of $n_0$.


Circle detection with incomplete edges I want to detect circles in the image below but this approach, of using cv2.approxPolyDP on contours and counting the sides of the polygon, doesn't work here since some of the circle is missing. Is there a faster way of detecting shapes in this case than using the Hough Transform?

this is John BG. Found a way to detect the traffic signal circle avoiding imfindcircles that requires the guessing of its parameter Sensitivity.

1.- Acquiring image

clear all;clc;close all

A=imread('001.jpg');
figure(1);imshow(A);
% 001
A1=A(:,:,1);
[sz1 sz2]=size(A1)     % sz1: Y axis  sz2: X axis

th1=120                     % binarize
A2=A1;
A2(A2>th1)=255;
A2(A2<=th1)=0;
% figure(2);imshow(A2)

A3=bwmorph(A2,'skel',Inf)             % clean and thin
figure(3);imshow(A3)
% 002

A4=bwlabel(A3,8);                         % number objects
hf4=figure(4);imshow(A4)

range_A4=unique(A4)';                      %  count objects
amount_segments=max(range_A4)

PA={}

PA: each element of cell PA contains all pixels of 1 segment

CM= [1,1,0;  % yellow
         1,0,1;  % magenta
         0,1,1;  % cyan
         1,0,0;   % red
         0,1,0;   % green
         0,0,1;   % blue
         1,1,1]   % white
nCM=[mod([1:1:amount_segments],size(CM,1))];
nCM(nCM==0)=size(CM,1)

for k=1:1:amount_segments % put pixels same object in same cell and colour objects
    [col1,row1,v]=find(A4==k);
    P_segment=[col1 row1];
    PA{k}=P_segment
    hold all                                                           % check all found points belong to same segment
    for m=1:1:length(P_segment)
        plot(hf4.CurrentAxes,P_segment(m,2),P_segment(m,1),'*','LineStyle',':','Color',CM(nCM(k),:))
    end

end
% 003

nPA=0;

confirming that majority objects are small and of no interest

for k=1:1:size(PA,2)
    nPA=[nPA size(PA{k},1)]
end
max(nPA)
figure(5);histogram(nPA,max(nPA))
nPA(1)=[]
% 004

sizeCelln=@(n,cell1) size(cell1{n},1)

n2 format: [amount pixels in given segment  , unique numeral identifying each segment in A4]

n2=[0]
for k=1:1:amount_segments
    if sizeCelln(k,PA)>floor(max(nPA)/3)
        n2=[n2 k ];                % k is same as A4(PA{k}(1,1),PA{k}(1,2))
    end
end
n2(1)=[];
PA2=PA(n2)

keep only the ones large enough to potentially look like circle

close(hf4);
hf4=figure(4);imshow(A4)
hold all
for k=n2                                                    % check all found points belong to same segment
    L3=PA{k}
    for m=1:1:size(L3,1)
        plot(hf4.CurrentAxes,L3(m,2),L3(m,1),'b*')
    end
end

Following, the best candidates (by amount pixels only) to circle

s=[0 0 0 0 0 0 0];

s format: [numeral mean(R) var(R) range(x) range(y) midpoint(x) midpoint(y)]

for k=n2
    L2=PA{k};
    x_midpoint=.5*(min(L2(:,1))+max(L2(:,1)))
    y_midpoint=.5*(min(L2(:,2))+max(L2(:,2)))
    R=((L2(:,1)-x_midpoint).^2+(L2(:,2)-y_midpoint).^2).^.5
    s=[s; k mean(R) var(R) range(L2(:,1)) range(L2(:,2)) x_midpoint y_midpoint]
end
s(1,:)=[];

The segment(s) that best resemble(s) circle(s) has/have characteristics clearly distinctive against other geometries like straight lines, squares, rectangles or ellipses.

nid=s(:,1);
R0=floor(s(:,2));
Rvar=.1*floor(10*s(:,3));
x_range=s(:,4);y_range=s(:,5);
x_mpoint=s(:,6);y_mpoint=s(:,7)

T1=table(nid,R0,Rvar,x_range,y_range,x_mpoint,y_mpoint)

T1 =
  10×7 table
     nid      R0      Rvar       x_range     y_range     x_mpoint     y_mpoint
     ___    ___    _____    _______    _______    ________    ________
     24     105     31.2         226           199          122             217.5
     28      67    260.6          92           166          300              218
     29      70    245.8          95           165         407.5            217.5
     32      58    416.1        142           121          137              205.5
     35      52    280.1        126           107          139              204.5
     43      57    329.4        139           117         106.5            234.5
     46      51    288.3        124           105          106              235.5
     50      28    140.9         44             93          409               244.5
     54      85    733.5        241           108         121.5             280
     57      69    628.9        196            67          344               272.5

Circles have properties like:

single point centers
constant radius
range(x) ~ range(y)

Only segment points distributed along a circle meet these 3 conditions.

n=[1:1:numel(nid)]

R2half_x_range_ratio=R0-abs(x_range)/2
R2half_y_range_ratio=R0-abs(y_range)/2
dxdy=abs(abs(x_range)-abs(y_range))

hf8=figure(6);plot(n,R0,n,Rvar,n,dxdy); grid on
xticks(hf8.CurrentAxes,[1:1:numel(nid)])
hf8.CurrentAxes.XTickLabelMode='manual'
hf8.CurrentAxes.XTickLabel=string(nid)'
xlabel('A4 segment numerals')
legend(hf8.CurrentAxes,'R','var(R)','\Deltax-\Deltay')

In this context the key parameters to catch circles are lowest radius dispersion and at the same time among segments with largest amount of pixels.

Even with the horn marked with a red arrow in the previous figure, the circle found is fairly centred and the radius is correct.

lowest Radius dispersion

n10=find(Rvar==min(Rvar))
nid(n10)

C1=[T1.x_mpoint(n10) T1.y_mpoint(n10)]
R1=T1.R0(n10)
da=2*pi/100;
a=[0:da:2*pi];
Px1=R1*cos(a)+C1(2);Py1=R1*sin(a)+C1(1);
hold all;plot(hf4.CurrentAxes,Px1,Py1,'r','LineWidth',10)

So, this is it, the traffic signal circle has been correctly detected.

Following, additional criteria that have been checked, but are not as clear as the lowest radius centred variance:

Similar x and y ranges

n11=find(dxdy==min(dxdy))
nid(n11)

C2=[T1.x_mpoint(n11) T1.y_mpoint(n11)]
R2=T1.R0(n11)
Px2=R2*cos(a)+C2(:,2)
Py2=R2*sin(a)+C2(:,1)
hold all;
plot(hf4.CurrentAxes,Px2(1,:),Py2(1,:),'m','LineWidth',2)
plot(hf4.CurrentAxes,Px2(2,:),Py2(2,:),'m','LineWidth',2)

Segments with R Most similar to half ranges

RxRy=abs(R2half_x_range_ratio)-abs(R2half_y_range_ratio)
n12=find(RxRy==min(RxRy))
nid(n12)

C3=[T1.x_mpoint(n12) T1.y_mpoint(n12)]
R3=T1.R0(n12)

discarding those segments with Radius too small

Without removing shortest segments

this is, without removing any segment with radius 1/3rd or shorter than peak radius.

% Run lines 1 to 68 of this script, then hop here:

PA2=PA
n2=[1:1:numel(PA)]

T2=table(nid,R0,Rvar,x_range,y_range,x_mpoint,y_mpoint)
 T2 =
   85×7 table
     nid     R0      Rvar        x_range    y_range      x_mpoint     y_mpoint
     ___    ___    _____    _______    _______    ________    ________
      1       4        3.2            14            4                57               3
      2      10      37.4           37           16               85.5             9
      3      10      36.1           36           23               94              12.5
  ..

     83       2       2.2            7             6              431.5           422
     84       2       0.9            7             2              437.5           428
     85       3       0.5            9             7              448.5           430.5

The compact script is available by email, feel free to contact, email address  in Profile. Also welcome any correction comment or suggestion for further development of this answer.

Appreciating time and attention

John BG

%%%%%%%%%%%%%%%%%%%%%%

Comments:

1.- When too many small objects the graphs get too cluttered.

2.- The start picture does not have perspective correction.

This causing circles to look like ellipses.

3.- The above steps also allow correction of perspective angle, and even calculate dimensions of nearby objects.

4.- The traffic signal dish has to be relatively large to all other circle-like shapes which is the same as the photo has to be taken relatively close to the signal traffic, for these steps to.

5.- The command imfindcircles requires the parameter Sensitivity to go up to .95 in order to start detecting circles, but it oddly shows a really span of Sensitivity where just a valid circle is caught. Just a bit off 0.95 and as shown in the 2nd figure below a multitude of non-relevant circles are added.

Rmin=20;Rmax=200;

close(hf4);
hf4=figure(4);imshow(A4)
[centers, radii] = imfindcircles(A4,[Rmin Rmax],'ObjectPolarity','dark','Sensitivity',.98);
viscircles(centers, radii,'Color','b');

close(hf4);
hf4=figure(4);imshow(A4)
[centers, radii] = imfindcircles(A4,[Rmin Rmax],'ObjectPolarity','dark','Sensitivity',.98);
viscircles(centers, radii,'Color','b');
% 012


DFT for something similar to convolution I have the following problem:
Let $x$, $y$ be finite real valued sequences defined on $0, \ldots, N-1$ and let $g$ be a non negative integer.

Define:

$$z_n=\sum_{k=0}^n x_{n-k}y_k, \quad \text{also on $0, \ldots, N-1$.}$$

Note that the sum is from $0$ to $n$, not $N-1$. In addition, the DFT of $y$ is known in closed form.

Is there a way to write $z$ as some cyclic convolution, so that with the help of the convolution theorem $z$ can be calculated in $N \log N$ instead of $N^2$?

I tried following the convolution theorem proof but I get stuck:

\begin{align}
z_n&=\sum_{k=0}^nx_{n-k}y_k= \sum_{k=0}^{N-1}x_{n-k}\cdot \theta (n-k)\cdot y_k\\
\Rightarrow \tilde{z_q}&=\sum_{n=0}^{N-1}\sum_{k=0}^{N-1} x_{n-k}\cdot \theta (n-k)\cdot y_k e^{-2\pi iqn/N}\\
&=\sum_{n=0}^{N-1}\sum_{k=0}^{N-1} x_{n-k}\cdot \theta (n-k)e^{-2\pi iq(n-k)/N}\cdot y_k e^{-2\pi iqk/N}\\
&=\sum_{k=0}^{N-1} y_k e^{-2\pi iqk/N} \sum_{n=0}^{N-1} x_{n-k}\cdot \theta (n-k)e^{-2\pi iq(n-k)/N}\\
&=\sum_{k=0}^{N-1} y_k e^{-2\pi iqk/N} \sum_{n'=-k}^{N-1-k} x_{n'}\cdot \theta (n')e^{-2\pi iqn'/N}\\
&=\sum_{k=0}^{N-1} y_k e^{-2\pi iqk/N} \sum_{n'=0}^{N-1-k} x_{n'}e^{-2\pi iqn'/N}=?????\\
\end{align}

The problem is that the second sum depends on $k$ so the double sum doesn't factor to the product of DFTs.

What am I missing?

You might be looking for the convolution property of the DFT. If $z[n] = x[n] * y[n]$, where $*$ represents circular convolution, then $Z[k] = X[k] Y[k]$, where $Z[k]$ is the DFT of $z[n]$ and so on.

This property is used in fast convolution techniques like overlap-save to implement convolution in $O(N \log N)$ time instead of the $O(n^2)$ time required by direct implementation of the convolution sum.


Measuring amplitude of a spectral component I'm trying to measure the amplitude of a spectral component in a signal. The component's frequency is always known beforehand.

I take a given number of samples, apply the flattop window to reduce scalloping loss, then use FFT and look at the bin of interest.

The earlier method used adaptable sampling frequency ($f_s$ was an integer multiple of the signal frequency) -> no windowing was needed.

Compared to the earlier method, measurements with the new method have a lot more variance. I suspect that it is the window that picks up the noise (convolution over the whole spectrum, where there are other strong components).

The $f_s$ is either 312 KS/s or 625 KS/s, the signal frequency is 2-15 KHz.

How could I improve my measurement?

Is there an other method which exploits the fact that the frequency of the signal is always exactly known?

A variation on @geometrikal's answer would be to use the Goertzel Algorithm.

$$
y[n] = x[n] + e^{+j\omega_0} y[n-1]
$$
where $x[n]$ is your input signal, $\omega_0$ is the frequency you know, and $y[n]$ is the (complex-valued) estimate of the discrete Fourier coefficient.


Is the definition of recursion and iteration in signal processing different from computer science? I was talking to my advisor today about this. He is a very experienced professor in this field, so I assume he is very knowledgable about it.

He pointed out that in signal processing, the word iteration means the repeated process that you repeat to gain a more accurate result (not the exact word, but generally this meaning). On the other hand, the word recursion means the actions that go on with every time instant. Basically, iterations happen within the same time instant, and recursions happen along the time instants.

However, as far as I can remember from the algorithm course that I took in undergrad, recursion means a function that can call itself repeatedly until there's no data left. The wikipedia page of recursion is consistent with my understanding.

What I guess is the situation now is, probably the definitions of recursion and iteration are different in SP than CS. I wanted to know, how exactly they are defined in signal processing. It is also possible that there is a different layer of the meaning of recursion that I didn't know, which made the two ways of understanding consistent. If there was a textbook or material that has the formal definition of them, it would be great.

I wanted to bring up the wikipedia page today during the meeting, but we ran out of time. I tried to google but couldn't find the answer. Guess I will ask him about it next time when we meet. :)

In my understanding, the term iteration in DSP is usually used in the mathematical sense, i.e., for describing iterative methods solving certain numerical problems. The power method for computing the largest eigenvalue of a matrix and its corresponding eigenvector is a typical example of such an iterative method used in signal processing.

The term recursion is normally used in the sense of recurrence relation, where a value $y_n$ is generally a function of $n$ and of its values $y_k$, $k<n$:

$$y_n=f(n,y_{n-1},y_{n-2},\ldots)$$

In signal processing applications, such as filtering, the index $n$ is interpreted as a time index. E.g., when implementing a (causal) filter, the current filter output can be computed from the current input and previous input and output values. This is a specific form of a recurrence relation. On the other hand, when talking about iteration, the iteration index is usually not related to time. I think it is this latter point that your professor referred to.


Is there any tool to anlyze audio specturm of a song and give the average visually ? I have a lot of songs I made over the 20 years and was thinking it would be interesting to see how the overall EQ changed on the various songs I made.

Any tool out there to help me with something like this ?

Audacity will let you look at a song's spectrum/spectrogram, and it works with most audio formats.

It's also free!


Swift vDSP rFFT not same as Python np.fft.rfft() I am trying to implement real FFT in iOS, for that i am using Accelerate Framework. Here is my Code for Swift.
    class FFT {

private var fftSetup: FFTSetup?
private var log2n: Float?
private var length: Int?

func initialize(count: Int){
    length = count
    log2n = log2(Float(length!))
    self.fftSetup = vDSP_create_fftsetup(vDSP_Length(log2n!), FFTRadix(kFFTRadix2))!
}

func computeFFT(input: [Float]) -> ([Float], [Float]) {

    var real = input
    var imag = [Float](repeating: 0.0, count: input.count)
    var splitComplexBuffer = DSPSplitComplex(realp: &real, imagp: &imag)

    let halfLength = (input.count/2) + 1
    real = [Float](repeating: 0.0, count: halfLength)
    imag = [Float](repeating: 0.0, count: halfLength)

    // input is alternated across the real and imaginary arrays of the DSPSplitComplex structure
    splitComplexBuffer = DSPSplitComplex(fromInputArray: input, realParts: &real, imaginaryParts: &imag)

    // even though there are 2 real and 2 imaginary output elements, we still need to ask the fft to process 4 input samples
    vDSP_fft_zrip(fftSetup!, &splitComplexBuffer, 1, vDSP_Length(log2n!), FFTDirection(FFT_FORWARD))

    // zrip results are 2x the standard FFT and need to be scaled
    var scaleFactor = Float(1.0/2.0)
    vDSP_vsmul(splitComplexBuffer.realp, 1, &scaleFactor, splitComplexBuffer.realp, 1, vDSP_Length(halfLength))
    vDSP_vsmul(splitComplexBuffer.imagp, 1, &scaleFactor, splitComplexBuffer.imagp, 1, vDSP_Length(halfLength))

    return (real, imag)

}

func computeIFFT(real: [Float], imag: [Float]) -> [Float]{

    var real = [Float](real)
    var imag = [Float](imag)

    var result : [Float] = [Float](repeating: 0.0, count: length!)
    var resultAsComplex : UnsafeMutablePointer<DSPComplex>? = nil

    result.withUnsafeMutableBytes {
        resultAsComplex = $0.baseAddress?.bindMemory(to: DSPComplex.self, capacity: 512)
    }

    var splitComplexBuffer = DSPSplitComplex(realp: &real, imagp: &imag)

    vDSP_fft_zrip(fftSetup!, &splitComplexBuffer, 1, vDSP_Length(log2n!), FFTDirection(FFT_INVERSE));

    vDSP_ztoc(&splitComplexBuffer, 1, resultAsComplex!, 2, vDSP_Length(length! / 2));
    //
    //// Neither the forward nor inverse FFT does any scaling. Here we compensate for that.
    var scale : Float = 1.0/Float(length!);
    var copyOfResult = result;
    vDSP_vsmul(&result, 1, &scale, &copyOfResult, 1, vDSP_Length(length!));
    result = copyOfResult

    return result
}

func deinitialize(){
    vDSP_destroy_fftsetup(fftSetup)
}

Here is my code of Python for computing of rFFT and irFFT
    # calculate fft of input block
    in_block_fft = np.fft.rfft(np.squeeze(in_buffer)).astype("complex64")

# apply mask and calculate the ifft
    estimated_block = np.fft.irfft(in_block_fft * out_mask)

Question:

Swift: If I compute rFFT for a 512 frame and apply irFFT to the resultant of rFFT, I get the same original array.

Python: Same goes for Python, if I take rFFT and irFFT I get the original array in return.

Problem: The problem occurs if I compare the results of Swift rFFT and Python rFFT. Their results are different in decimal values. Sometimes the real part is same but the imaginary is totally different.

I tried different Framework in Python like Numpy, SciPy and TensorFlow and their results are exactly the same (Slight different in decimal part). But when I calculate rfft on same input for iOS using the Swift code above, results are different.
If anybody who has experience with Accelerate Framework and holds the knowledge of FFT help me out with this one would be very helpful. I have limited knowledge of FFT.

Sanity check your data layouts.
It is very common for FFT libraries to pack the (real-valued) Nyquist-bin into the unused imaginary component of the DC-bin and it happens almost naturally when you implement a real-transform as a half-size complex-transform.
As far as I can tell, vDSP wants such "packed" layout, where as Python apparently stores it as a separate "unpacked" real-component in a larger array.


Problem with applying elliptic filter to ecg signal I'm currently working on an article for ECG classification which it says that it has used elliptic bandpass filter with 0.5Hz and 50Hz critical frequencies to eliminate base-line and interference etc. First problem is that I think it should be stopband instead of bandpass, am I wrong?
Second problem is when I try to do the filtering in matlab. The original signal shape is something like this:

the sampling rate of the signal is 300Hz.
[b, a] = ellip(10, 1, 100,[0.5 50]/150, 'bandpass')
fvtool(b, a, 'Fs',300)

This is the frequence response of the filter:

using this line, I apply the filter to my signal:
s_filtered = filtfilt(b, a, s)

and what it returns is all NaN! Am I doing something wrong?

Numerical problems. 64-bit double precision floating point is not nearly enough to implement  this filter in the way you have done it.
Your filter is extremely steep: The poles are way too close to the unit circle and the order is way too high to implement it in transfer function form (which is a bad anyway).
Things to do

Implement the filter as cascaded second order section
Review the requirements for your filter. It seems way steeper than it needs to be. Such an aggressive filter will destroy a lot of time domain details since the filter has massive amounts of time domain ringing. The settling time alone is multiple seconds and by using filtfilt you double it again & create an enormous non-causality.

Filter design requires a lot of trade offs, make sure you understand them and optimize to the specific requirements of your application.


Help choosing coefficients for a simple FIR filter to extract DC component I am a beginner in design of digital filters... I have just completed Coursera's Linear Circuits course!

I have a signal which is sampled at 50 Hz where the DC component has a value of 800000 and there is some noise with a standard deviation of around 200. The DC component should remain constant but it may change very slowly over time.

I want to design a simple FIR filter to obtain the DC component: would a simple weighted average algorithm be most appropriate? How do I go about choosing suitable coefficients?

MA is almost never a "most appropriate" algorithm (it's only advantage is simplicity). My advice is to choose Gaussian filter which has least possible group delay. Calculate its coefficients keeping in mind that cut-off frequency of a Gaussian filter is $$f_0 = \frac{1}{2 \pi \sigma}$$
Estimate which variances of low frequency component you'd like to preserve in your application, which are acceptable to be attenuated to 50%, this will give you desirable cut-off frequency.


BandPass Signal Vs PassBand Signal Am reading Rick Lyons - Understanding DSP book about Sampling. I have few doubts related to this chapter.

what is the difference between Bandpass Signal and a Passband Signal ?
Can somebody give an example for a Bandpass system which produces Bandpass signal as output !
When it comes to Baseband Sampling ( Low pass Sampling ) sampling theorem states that " Fs > twice the Highest frequency Component" and for Bandpass Sampling, sampling theorem says Fs > twice the Bandwidth of the signal. why two different sampling statements for the same signal of B hertz.

I understand in Bandpass sampling if I sample at the rate of twice the highest frequency component, then it'll be oversampling.

A PassBand is the Band that can be passed through a BandPass filter, or for that matter, any filter. In other words, the filter is a Bandpass, what signal passes through is a Passband.
The classic is a high pass and low pass together. For instance, a capacitor, inductor, and resistor in series is an example of a bandpass filter. Digitally, you combine a high pass and a low pass filter together.
Essentially, if you are limited to a small region, you can infer what frequency you are limited to. The classic example is Nf->2*Nf. The frequencies reported at 0 is 2*Nf, the frequency reported at Nf is actually NF, and continue the pattern for the rest. This only works when you are guaranteed to only have frequencies in the particular band.


Is the value of a frequency bin of a DFT-output the average of the 'real' frequency values within that bin's range? I am wondering whether the value of a frequency bin with a certain resolution is the average of the fourier transform values of the 'real' frequencies within that bin's range.

I'll start with a counter-example:

Assume your continuous-time signal is a sinusoid at freqeuncy $f_0$. You then sample that signal, and use a rectangular window to select $N$ samples. Feed those $N$ samples to your DFT, and you will see non-zero values in the whole spectrum (unless in the special cases where $f_0$ is a multiple of $f_s/N$).

So you start with a signal whose energy is concentrated in a single frequency, and you end up with non-zero values in all DFT bins.

To understand what is happening:

Start with a continuous-time signal $x_c(t)$.
(Optionally) filter the signal to make it band-limited to $f_s/2$.
Sample the signal to get $x[n] = x_c(n/f_s)$.
At this point, the spectrum is given by the sampling theorem. You will have spectrum aliasing if the sampled signal was not band-limited in step 2.
Now apply a window of size M. The window type depends on your application (Rectangular (ugly), Hann, Hamming, Blackman, Flat top, etc.).
At this point, your spectrum is that of the sampled signal, but convolved with the window (depending on window size and type, you will get some energy spread to other frequencies). This spectrum is $X_w(e^{j\theta})$.
The DFT of the windowed signal will return samples of its spectrum, at frequencies $Y[k] = X_w( e^{j 2\pi k/N})$.


Estimate transfer function from Bode curve I have measured the magnitude response Y_mag(f) and phase response Y_phase(f) of an unknown physical system.

Is it possible to estimate the s-domain transfer function without a priori information about the order of the polynomials in the numerator or denominator?

I am more or less looking for an out-of-the-box approach that estimates the transfer function in a sufficient way. Performance and computing time are not an issue.

Python and MATLAB (Signal Processing Toolbox, Control System Toolbox) are available.

So this page suggests that you might be able to use this package in Python.

The main issue is that that package seems to require a time series for estimation, rather than the frequency domain data you have.  To get around that, you may be able to use your frequency domain data to generate a time series.

Suppose your frequency domain data is
$$Y(\omega = 2\pi f) = Y_{\tt mag}(f)\exp(\imath Y_{\tt phase}(f)), \ \ f>0$$
and the complex conjugate for $f<0$.

Interpolate your data:

$$
y(t) = {\tt IFFT}(Y(\omega))\\
y_{\tt padded}(t) = \left \{
\begin{array}{ll}
y(t) & \mbox{for } 0 \le t \le T_y\\
0 & \mbox{for } T_y \le t \le 2T_y \mbox{  (say)}
\end{array}
\right .\\
Y_{\tt padded}(\omega) = {\tt FFT}(y_{\tt padded}(t))
$$
then just use this to filter a padded white noise sequence of length $T_y$ (padded to $2T_y$).


Most accurate method of Voice Activity Detection? What is currently the most accurate Voice Activity Detection algorithm with high levels of background noise? I would also appreciate if it does not require any manual calibration. Bonus points for examples.

EDIT: I am referring to accuracy as yielding the lowest amount of false positives and negatives.

If you have more than 1 detectors, you may try blind separation based on the power spectrum difference of voice and background noise


How do you get FFT for negative $k$ values? I am reading a book (am a programmer so I suck at math) and it states that for a number of $k$-values that are symmetric around 0 (for example $k = -3, -2, -1, 0, 1, 2, 3$), we need to calculate

$$\sum_{j = 0}^{N - 1} x_j e^{-2\pi i kj/N}$$

and the book claims this can be done with FFT.

I know nothing about FFT but Wikipedia tells me that FFT can calculate above sum for $k = 0, ..., N - 1$ in $O(N\log N )$-time.

How does the book then intend to apply FFT to get the values for the $k$ we are interested in here? We want the sum for $k = -3, -2, -1, 0, 1, 2, 3$, but the FFT will calculate them for $k = 0, 1, 2, 3, 4, ..., N - 1$.

There are N possible circular rotations but consistent with what FFTSHIFT in Matlab/Octave provides: Simply map the FFT according to:

For FFT k values before fftshift: 0 to N-1

The samples from 0 to floor((N-1)/2) would map to ceiling((N-1)/2) to N-1

The samples from ceil(N/2) to N-1 would map to -ceiling((N-1)/2) to -1

For example if you have an eight sample FFT  (N= 8)

floor((N-1)/2) is floor(7/2) = 3, so the samples from 0 to 3 map to 4 to 7

And the samples from 4 to 7 map to -ceiling(7/2) to -1 or -4 to -1

The samples at index k: 0 1 2 3 4 5 6 7

get mapped to: 4 5 6 7 0 1 2 3

corresponding to the new index -4 -3 -2 -1 0 1 2 3

Similarly when the FFT is odd, for example a nine sample FFT (N= 9)

the samples at index k: 0 1 2 3 4 5 6 7 8

get mapped to: 5 6 7 8 0 1 2 3 4

corresponding to the new index:  -4 -3 -2 -1 0 1 2 3 4


The Concepts Behind SVD Based Image Processing I am interested in information about the output of SVD of the matrix. Mathematical part of SVD process is clear to me but I don't understand how to read figures with graphs of SVD of a matrix. As I understood vertical axis represents singular values.

But I can't figure it out what horizontal one represents.
Also, what is the difference between singular spectrum and singular
vectors representations?
What singular spectrum means, is it singular values?

You may think on the SVD as a generalization of the Discrete Fourier Transform.
Namely, it is generates an orthogonal basis to represent the data.
The nice thing about it, it generates the basis according to data (Where the Discrete Fourier Transform basis is the same for any data).
Just like the Fourier Spectrum, you have the "Energy" - The eigenvalue.
The i-th eigen value represent the projection of the data on the i-th basis.
In the image above they show you some Eigen Vectors, the basis functions, of the data. As you can see, it is very similar to Fourier Basis which means the data is well represented by the harmonic discrete signals.

The horizontal axis is the sample number of the Eigen Vector - The basis.
The spectrum is the value of the projection. The vectors are the basis to project upon.
Yes. It is the set of the Singular Values.


FFT Phase interpretation of input signal with non-integer number of cycles in FFT window It is known that if we take the FFT of a signal at such a frequency which has non-integer number of cycles in the sample window we get different value of phase even if the original signal was at phase 0 (reference to cos wave). Looking at the time-domain signal can we predict what the phase output is going to be i.e. as there any intuitive sense behind the phase that is reported?

There is a circular waveform discontinuity at sample 0 of an FFT input (to sample N-1), if it isn’t exactly integer periodic in aperture.  However, if the waveform is continuous at sample N/2, then the phase can be measured at that point.

You can measure phase at the halfway point by doing an FFTshift, or by flipping the phase of every odd numbered FFT output bin.  If you know, or can estimate the frequency, then if you estimate the phase at sample N/2, you can use the frequency and phase at a known point to compute the phase at any other point in time.  Including at the beginning of the original sample data window.

This answer is for even N Lengths of FFT.  (Dawg may have a better answer for odd lengths.)


In Plain English, What Are the Within and Between Class Variances in Otsu's Thresholding? For the past week I've been through numerous research papers and tutorials on Otsu's thresholding. Even wrote to the authors, but every single one of the explanations are in terms of "Class" and "variance". I'm not very good at statistics, so it's hard for me to picture what these classes and variances mean in terms of the actual image. I have a few questions (marked in bold) that lead to the main questions.

What I've understood until now:
1. We need to separate the image into two parts. The dark pixels and the light pixels. These are what are called the two classes. For this explanation could we please just call them "dark" and "light" instead of "class".
2. A histogram array of 0 to 255 positions is created, where each position holds the total number of pixels of intensity 0, 1, 2,...255. Let's say for intensities 0 to 9, the number of pixels are 20,5,23,6,3,23,5,8,67.
3. Probability for any of these intensities (assume intensity 1) are calculated as = 1 / (180 * 200), where 180 is the length of the image and 200 is the width.
4. The sum of probabilities of the dark and light pixels are the weights of the dark and light pixels.
5. Question: When we calculate the mean of the pixels, what exactly does it imply? Is it like taking an average of the probabilities? If it is like an average, then if the threshold is at position 100 of the histogram, why can't we simply calculate it as: (sum of intensities from 0 to 99) / (99) ? What extra benefit do we have by doing a complex calculation of mean?
6. Then we have variances of distribution of pixels from the mean, which I understand, but Question: when we take the sum of variances of the light pixels and sum of variances of dark pixels, what does it imply? Does it speak of the amount of variation of intensity present in the dark pixels vs the light pixels?
7. Question: What is a "within class variance" and a "between class variance" in plain English? I can see the formula, but it makes no sense. Is the within class variance just a way of saying "how much the dark pixels vary with respect to mid-darkest pixel plus how much the light pixels vary with respect to the mid-lightest pixel"?
8. Question: What is a "between class variance" in plain English? The formula itself is odd, because there is no variance, but just a sum of the wights and means. But still, what does it mean in plain English, in terms of the pixels?

ps: Kindly refrain from statistical/mathematical terms and please try to explain it in plain English the way I've put some explanations in quotes above. I (and many people in this world) have a hard time understanding statistical jargon because it was taught in a way that never allowed us to fully understand it.

i'm avoiding labels light and dark because you could say that we choose 128 as threshold which perfectly separate light and dark pixels but consider you have a picture with only two gray levels which are 240 and 250 and you want to label them but both seems to be light.

now consider an example which you have a part of x-ray image which composed of soft tissue and bone and you want to separate them. we know bones generally have higher intensity than soft tissue in x-ray image. but some part of bones are very thin and some part of soft tissue is very thick so some pixels which  belong to bone have lower intensity than pixels belong to a soft tissue. in this problem any threshold that we choose lead to some miss labeled pixels.

lets consider more details. bones have thick and thin areas which lead to a variance around average intensity of bones, this variance is called within class variance and the average value represent the bone intensity value. the same is true for soft tissue. the difference between average of soft tissue and average of bone is called between class variance.

in ideal situation we have minimum within class variance and the averages must be as far as possible from each other (highest between class variance) to have the minimum mixing between these two groups.


Time scaling in DFT? Let $x[n]$ be a real signal whose length $N$ is even. Let $y[n]=x[2n]+jx[2n+1]$ a signal whose length is $M=\frac{N}{2}$. $X[k]$ for $k=0,...,N-1$ and $Y[k]$ for $k=0,...,M-1$ are the respective DFTs of those signals. Find $Y[k]$ as a function of $X[k]$.

I haven't been able to find a scaling time property for the DFT. I read all the Oppenheim-Schafer's discrete signals book and found nothing about the subject. I tried doing it by definition but I just can't realize what the relationship between both DFTs is.

This is a classic exercise and it shows you how to obtain the DFT of a real-valued length $N$ signal in terms of the DFT of a length $N/2$ complex-valued signal.

You have to write $X[k]$ in terms of the even and odd samples of $x[n]$. With $W_N=e^{-j2\pi /N}$ you get

$$\begin{align}X[k]=\sum_{n=0}^{N-1}x[n]W_N^{kn}&=\sum_{n=0}^{N/2-1}x[2n]W_{N/2}^{kn}+W_N^k\sum_{n=0}^{N/2-1}x[2n+1]W_{N/2}^{kn}\\&=\text{DFT}_{N/2}\{y_R[n]\}+W_N^k\;\text{DFT}_{N/2}\{y_I[n]\}\\&=Y_e[k]-jW_N^kY_o[k]\tag{1}\end{align}$$

where $y_R[n]$ and $y_I[n]$ are the real and imaginary parts of $y[n]$, respectively, and $Y_e[k]$ and $Y_o[k]$ are the even and odd parts of $Y[k]$. Eq. $(1)$ shows that after having computed the length $N/2$ DFT $Y[k]$, you need to split it into its even an odd parts, and from them you can obtain the length $N$ DFT $X[k]$.

Obtaining $Y[k]$ from $X[k]$ is also possible, even though that's usually the less practical problem, because often you have an FFT routine that takes complex input signals, so it's more efficient to compute a lenght $N$ FFT of a real-valued sequence by computing a length $N/2$ FFT of a complex-valued sequence. Anyway, $Y[k]$ is given by

$$Y[k]=\sum_{n=0}^{N/2-1}x[2n]W_{N/2}^{kn}+j\sum_{n=0}^{N/2-1}x[2n+1]W_{N/2}^{kn}\tag{2}$$

With $W_{N/2}^n=W_N^{2n}$, the first term on the right-hand side of $(2)$ can be written as

$$\begin{align}\sum_{n=0}^{N/2-1}x[2n]W_{N}^{k(2n)}&=\sum_{n=0}^{N-1}x[n]W_N^{kn}\frac{1+(-1)^n}{2}\\&=\frac12\sum_{n=0}^{N-1}x[n]W_N^{kn}+\frac12\sum_{n=0}^{N-1}x[n]W_N^{nN/2}W_N^{kn}\\&=\frac12\left(X[k]+X[k+N/2]\right)\tag{3}\end{align}$$

where I used $(-1)^n=W_N^{nN/2}$.

Eq. $(3)$ shows that decimation in the time domain corresponds to aliasing in the frequency domain, as expected.

In a completely analogous manner, the last term in Eq. $(2)$ can be written as

$$j\sum_{n=0}^{N/2-1}x[2n+1]W_{N/2}^{kn}=\frac{j}{2}\left(X[k]-X[k+N/2]\right)W_N^{-k}\tag{4}$$

Substituting $(3)$ and $(4)$ into $(2)$ gives the final solution:

$$Y[k]=\frac12X[k]\left(1+jW_N^{-k}\right)+\frac12X[k+N/2]\left(1-jW_N^{-k}\right)\tag{5}$$


Envelop detection using matlab hilbert function I am trying to detect the envelope of a high-frequency signal (similar to below)  using Matlab Hilbert function

This is my output

It generally follows the envelope, but there is a high-frequency ripple in the envelope, see zoomed version below

Here is my code:
t =out.tout';
x =out.Itxmod';
y = hilbert(x);
env = abs(y);
figure;
plot_param = {'Color', [0.6 0.1 0.2],'Linewidth',2};
plot(t,x)
hold on
plot(t,[-1;1]*env,plot_param{:})
hold off
xlim([0.001 0.002])
title('Hilbert Envelope')

Question: How can I get a smooth envelope the signal without this high-frequency ripple?

This must be an artifact of the implementation for the hilbert function, in that any implementation of finite length will have a finite ripple. It doesn't appear the underlying filter using in the hilbert function can be modified, so there are three suggestions to obtain a smoothed demodulated result:
Method 1: Implement the Hilbert using the FFT: take the FFT of the original signal, set all the negative frequencies to zero (the upper half of the FFT) and double the positive frequencies leaving bin 0 (DC) as is (and as Overlord pointed out in the comments, leave the Nyquist bin at $N/2$ as is when the FFT length $N$ is even). The IFFT will be the Hilbert Transform and taking the absolute value of this signal will be the desired envelope. (As @aconcernedcitizen pointed out in the comments, this IS the method specifically implemented by Octave using the hilbert function from the signal package.)
Method 2: Implement the Hilbert with quadrature phase tracking filters where the filter length versus ripple can be traded.
Method 3: Use a traditional AM demodulator by multiplying the modulated signal with the coherent carrier and then low pass filtering the result.
A demonstration of Method 1 is shown below, which is a suitable solution for the OP's case:

Zooming in on the start of the signal shows the distortion limited to the start-up condition and the smooth envelope after that.

The Matlab code for this is as follows:
N = length(sig);
sig_spectrum = fft(sig);
hilbert_spectrum = zeros(1, N);
hilbert_spectrum(1) = sig_spectrum(1);
hilbert_spectrum(2: ceil(N/2) - 1) = 2 * sig_spectrum(2:ceil(N/2) - 1);
if mod(N, 2) == 0
  hilbert_spectrum(N/2) = sig_spectrum(N/2);
endif
hilbert_time = ifft(hilbert_spectrum);
envelope = abs(hilbert_time);


What is the index modulation in OFDM and Subcarrier Index modulation? I'm new here, I'd like to ask about the OFDM-IM (OFDM - Index modulation), what is it and how it works?  Is it the same as SIM-OFDM?

I checked white papers, I see it's using the index of sub-carriers for transmitting bit too. But what does it mean by the index of sub-carriers?

Could you please explain it in clear way? or share any article where I can understand it well.

I'm currently trying to understand R. Abu-alhiga and H. Haas, "Subcarrier-index modulation OFDM," 2009 IEEE 20th International Symposium on Personal, Indoor and Mobile Radio Communications, Tokyo, 2009, pp. 177-181.

Yes OFDM-IM is somehow new and it's very hot topic of searching. That article you have shares is the basic of OFDM-IM. So if you are at the beginning to understand it, that's the right paper to to read and understand.

OK let me explain that paper for you in easier way,

consider we have block with length 32 bits you want to transmit it, in the traditional OFDM, using 4-$QAM$ modulation, we supposed to have $32/2 = 16$ FFT size (which represents the number of sub-carriers). According to that figure 01 in the article you've shared, all sub-carriers will be active to convey data to receiver following the common steps of OFDM.

Now, let's look at OFDM-IM, using the block whose length is $32-bits$ and 4-$QAM$. In OFDM-IM, we start by dividing our bits into two parts. each with length of $16-bits$ The first part will be handled similar to traditional OFDM, but the second part will be handled using OOK and using also that module added in Figure 01 called (sub-carrier index modulation).  How??

That second part, we start count how many zeros and ones we have in that block (it's called subblobk). suppose we have 10 ones and 6 zeros similar to that example in your article,That's ok, but the issue is how can we know the location of zeros and locations of ones? That's what the subcarrier index means!!!  It means that the index of sub-carrier will indicate that location.

For example, in that example, you've 16 FFT size representing the number of sub-carriers, so we use each index to point that location of zeros. In other words, we use the subblock in 4-$QAM$ to show the location of zeros in subblock OOK. It means, we look at first bit in OOK, if it's 0 we let the sub-carriers to be empty, if it's 1 we upload the bits of $QAM$ into sub-carrier and so on. So, it's called OOK fashion.

That's what the idea of paper you've shared, I will let the other parts of paper to you to analyze and understand. then you notice easily the gain of using OFDM-IM

Hope that will help


How is the complexity of applying Short term fourier transform over n samples with a moving window of size m amount to O(nmlogm) I have been trying to understand a publication here where Short time fourier transform is applied over n samples in steps of m samples each (m is the size of the moving window). I understand that the computational complexity of the Short time fourier transform for m samples is mlogm.

The operation also consists of extracting energy from the samples. Since there are m samples in the window, I believe it should be O(m) in the worst case.
So , as per my understanding, the complexity of the whole operation should be
O(n/m * (mlogm + m)) = O(nm). Since there n/m windows and each window has mlogm for fourier transform and m for signal strength extraction.

But the publication shows the complexity as O(nmlogm). Please help me understanding this.

evidently the author is performing an FFT $n$ times (or, if the window hop is not 1 sample, it would be proportional to $n$) and each FFT costs something proportional to $m \log(m)$).  all this is quite normal.


On the use of Permutation matrix to perform iFFT I have a question about using the permutation matrix for performing $iFFT$ for such matrix and then reshape it row-wise and column-wise way.

Let's say that we have a random matrix $x$ whose size is (32,8), so that matrix can be written as below:

where $x_k$ represents the column $k$ of matrix $x$ with size of 32 for each column.

$1$- Let's take column wise $iFFT$ for matrix $x$ and then reshape it row-wise resulting matrix $Y$ as below matlab code

x = randn(32,8);
y = ifft(x);
Y = reshape(y.',[],1);

$2$- let's manipulate that using another way of permutation matrix. first we can reshape the matrix $x$ column-wise into (256,1) as below:

Then performing $iFFT$ operation with size of $32$-point is corresponding to:

$F$ is the Fourier transformation matrix
In order to get the equivalent matrix of Y representing the reshape of $iFFT(x)$ in row-wise, we can multiply the matrix Z with permutation matrix P resulting vector $Y2$ which should be similar to vector $Y$ in case $1$. That can be done by matlab as:

z = reshape(x,[],1);   %
X = kron(eye(8),dftmtx(32));
Z = X*z;

P = zeros(m*n);       %Building the permutation matrix P
col = 1;
for i = 1:m
    for j = 1:n
        E = zeros(m,n);
        E(i,j) = 1;
        P(:,col) = E(:);
        col = col + 1;
    end
end

Y2 = P*Z

So, for $(1)$ and $(2)$,   $Y$ should be similar to $Y2$. Is that right? but what I get is different results!  My question is that calculation right? why do I get different results ?

Y and Z should show the same values

clear;
%Create random vector x
x = randn(32,8);

%take ifft along each column
y = ifft(x);
Y = reshape(y,[],1);

z = reshape(x,[],1);   %
d=conj(dftmtx(32))/32;
X = kron(eye(8),d);
Z = X*z;


How to determine the number of multipath? In a multipath propagation environment if we take the time index of reach channel tap by the integer multiples of sampling periods, assume the RMS delay spread $ \sigma_{tau} $ and maximum excess delay $ \tau _m $ is known, then how to determine the maximum number of paths?

You can't really tell that from anything. The maximum delay spread together with your sampling rate, however, gives you the number of paths that you can tell apart with a single observation of the channel. If you model every sample of your discrete-time channel impulse response as a path.
That's not really a sensible "channel model", however: the tapped delay line model works well enough for point reflectors and static channels. However, your channel changes, and your reflectors are moving: you get paths as "smeared" peaks in your ambiguity function.
Still, the tapped delay line model with each sample of the channel impulse response being one tap is immensely useful, as it allows you to equalize the channel by OFDM.
Many scenarious consider point reflectors (which sadly don't always sit "exactly" sampling periods apart) that are moving, thus "smearing" and leaking across multiple taps. As soon as you add anything that starts to track different path components, you might be able to sense more paths than you have taps!
So, honestly, there's no "maximum number of paths" that you can distinguish if you have enough time and if they are uncorrelated enough.


Discard circular/blob like structures from image I have a binary image that looks like this:

And I would like to get rid off all the small structures and just keep what's on red:

Image erosion is an option that works well to remove the very small groups of pixels without affecting the red regions. Is there any other method that can be applied to keep only the elongated structures on my image?

You can calculate the area of each connected area in your image. Then, you throw away the smaller areas and only keep the biggest. See this example Python code:

import skimage.io as io
import numpy as np

from skimage.measure import regionprops, label
from skimage.color import label2rgb

I = io.imread("http://i.stack.imgur.com/rtyFJ.png")
bw = I[:,:,0]

bigger = np.zeros_like(bw)
labels = label(bw)

for R in regionprops(labels):
    if R.area > 500:
        # draw the region (I'm sure there's a more efficient way of doing it)
        for c in R.coords:
            bigger[c[0], c[1]] = 1

plt.imshow(bigger)

You can play around with the area threshold, e.g. by calculating a histogram of all areas and choosing the 3 biggest etc. And, you need to define the neighborhood of a pixel to belong to a connected region. In my case, it uses 8-neighborhood (i.e. if two pixels are connected at a corner, they belong to the same region). This is up to your application.

Additionally, you might throw away regions based on other criteria. E.g. if you want to throw away "circular" structures, you might throw regions away, where the area and perimeter are close to pi for example. Or, you look at the extent or orientation or whatever.


Intuition behind the scaling property of Fourier Transforms The Fourier transform of $f(ax)$ is $\frac{1}{|a|}F(\frac{u}{|a|})$. So the frequencies are scaled horizontally but the magnitudes are also scaled when the graph of $f$ is scaled horizontally.

On the other hand the fourier series of $f(\theta)$ on $[-\pi,\pi]$ is $\sum_{-\infty}^{\infty} c_n e^{in\theta}$ whereas the fourier series of $f(a\theta)$ on $[\frac{-\pi}{a},\frac{\pi}{a}]$ is $\sum_{-\infty}^{\infty} c_n e^{ina\theta}$ so that the frequencies are  scaled horizontally but the magnitudes are left unchanged at $c_n$.

Now these are easy to show mathematically, but I'm left feeling uneasy conceptually. Is there same intuition geometrical view to leave my intuition satisfied that explains why there is this lack of symmetry?

Let $g(x) = f(ax)$.  Then, using a change of variables $y = ax$,
$$\begin{align*}
G(u) &= \int_{-\infty}^{\infty} g(x)\exp(-j u x)\,\mathrm dx\\
&= \int_{-\infty}^{\infty} f(ax)\exp(-j u x)\,\mathrm dx\\
&= \frac{1}{|a|}\int_{-\infty}^{\infty} f(y)\exp(- j (u/a) y)\,\mathrm dy\\
&= \frac{1}{|a|}F\left(\frac{u}{a}\right)
\end{align*}$$
not $\displaystyle \frac{1}{|a|}F\left(\frac{u}{|a|}\right)$ as you say (cf. Henry Gomersall's
comment on your question).
Now, if $f(x)$ is a periodic function with period $2\pi$, then its
Fourier transform exists only in a generalized sense (involves impulses).
The Fourier series
representation $\displaystyle \sum_n c_n \exp(j nx)$ of $f(x)$
has Fourier coefficients
$$\displaystyle c_n = \frac{1}{2\pi}\int_{-\pi}^{\pi} f(x)\exp(-j n x)\,\mathrm dx,
 ~ -\infty < n < \infty$$
and these also show up in its Fourier transform
which is given by
$\displaystyle F(u) = \sum_n c_n \delta(u -n)$. Note that
the inverse Fourier transform is
$$\int_{-\infty}^\infty F(u)\exp(jux)\,\mathrm du
= \sum_n \int_{-\infty}^\infty c_n \delta(u-n)\exp(jux)\,\mathrm du
= \sum_n c_n \exp(j nx) = f(x).$$
The function $g(x) = f(ax)$ of period $2\pi/a$ thus has Fourier transform
$$G(u) = \frac{1}{|a|}\sum_{n=-\infty}^{\infty} c_n \delta(u/a -n)
= \frac{1}{|a|}\sum_{n=-\infty}^{\infty} c_n \delta\left(\frac{u-an}{a}\right).$$
But, $\delta(\alpha x) = \frac{1}{|\alpha|}\delta(x)$, and so with
$\alpha = a^{-1}$, we have
$$G(u) = \sum_n c_n \delta(u - an)$$
whose inverse Fourier transform is
$$\int_{-\infty}^\infty G(u)\exp(jux)\,\mathrm du
= \sum_n \int_{-\infty}^\infty c_n \delta(u-an)\exp(jux)\,\mathrm du
= \sum_n c_n \exp(j nax) = g(x) = f(ax).$$

So, you are not missing anything. Scaling an ordinary
function scales both
the magnitude and the frequency axis of the Fourier transform,
but when impulses are involved, scaling the argument of the
impulse leads to a magnitude scaling that exactly compensates
for the magnitude scaling imposed on the Fourier transform,
and leaves the Fourier coefficients unchanged.
frequency scaling also scales the amplitude


How is this was derived? DTFT for symmetric pulse do not understand how the middle expression was derived from above expression,

after all summation\multiplication i get:

$$\frac{e^{-jwN} + e^{jwN} - e^{-jw(N+1)} + e^{jw(N+1)}}{2 - e^{-jw} - e^{jw}}$$

First, the following trick will be useful:
$$ e^{j a} - e^{j b} =  e^{j \frac{a+b}{2}} \left( e^{j \frac{a-b}{2}} - e^{j \frac{b-a}{2}} \right) $$

Then simply expand the sum as:
$$ X(e^{j\omega}) = \sum_{n=-N}^{N}e^{-j\omega n} = \frac{ e^{-j\omega (-N)} - e^{-j\omega (N+1)} }{1 - e^{-j\omega}} = \frac{ e^{j\omega N} - e^{-j\omega (N+1)} }{1 - e^{-j\omega}} $$

Apply the trick to the numerator and denominator by recognising $a$ and $b$ carrefully: for the numerator $a=\omega N$ and $b=-\omega (N+1)$ whereas for the denominator $a=0$ and $b=-\omega$ which yields:

$$ X(e^{j\omega}) = \frac{ e^{j\omega N} - e^{-j\omega (N+1)} }{1 - e^{-j\omega}} = \frac{ e^{j \frac{\omega N -\omega (N+1) }{2}} \left( e^{j\frac{\omega N +\omega (N+1) }{2}} - e^{-j\frac{\omega N +\omega (N+1) }{2}} \right) }{ e^{-j\omega/2} \left( e^{j\omega/2}  - e^{-j\omega/2} \right) }$$

which simplifies to:
$$ X(e^{j\omega}) = \frac{ e^{-j\omega/2}  \left( e^{j\frac{\omega N +\omega (N+1) }{2}} - e^{-j\frac{\omega N +\omega (N+1) }{2}} \right) }{ e^{-j\omega/2} \left( e^{j\omega/2}  - e^{-j\omega/2} \right) } = \frac{  e^{j \omega \frac{2N+1 }{2}} - e^{-j \omega \frac{2N+1}{2}} }{   e^{j\omega/2}  - e^{-j\omega/2}  }  $$

which is:
$$ X(e^{j\omega}) = \frac{  \sin( \omega (N+\frac{1}{2}) )}{   \sin( \omega \frac{1}{2}) }  $$


Formula for designing Lynn's low pass filter For the second order Lynn's low pass filter, the general form of the transfer function is:
$$
H(z)=\frac{(1−z^{-m})^2}{ (1−z^{-1})^2 }
$$
where $m$ is a positive integer. The gain for this is $m^2$, time delay is $m-1$ samples and the nominal frequency corresponding to the location of the lowest frequency zero is $f_l = f_s /m$. The frequency response is given by:
$$
|H(e^{j\omega T})|=\frac{\sin^2(m\omega /2f_s)}{\sin^2(\omega /2f_s)}
$$
At the 3 dB cutoff frequency, for unit gain we have:
$$
|H(e^{j\omega_c T})|=\frac{\sin^2(m\omega_c /2f_s)}{m^2\sin^2(\omega_c /2f_s)} \\ \therefore  \frac{1}{\sqrt2}=\frac{\sin^2(m\pi f_c /f_s)}{m^2\sin^2(\pi f_c /f_s)}
$$
Now, to design the filter for our required cutoff frequency $f_c$ and sampling frequency $f_s$, we need to find the value of $m$ corresponding to these inputs.
So how can we get this equation to the form $m = f(f_c, f_s)$, so that we can avoid getting the $m$ value by trial-and-error or numerical methods? Or is there some other quick approximation formula that we can use?

Just analyzing your equation, let $a = f_c / f_s$
Also we know that for $m = k/(2a)$ it reduces to $1/(m^2 sin(\pi a))$ for odd $k$ and $0$ for even $k$.
$$\lim_{m \to 0} \frac{sin^2( \pi m a)}{m^2 sin(\pi a)} = \frac{a^2 \pi^2}{sin(\pi a)} = \frac{1}{\sqrt{2}}$$
Even for this case, determining $a$ is not an easy task. I computed numerically the curve defined by this equation.


Processing a heterogeneous image-collection for CBIR I'm working with a grey-scale data-collection of Hebraic tombstones in content-based image retrieval CBIR). For example an image like this tombstone. I'm going to reduce the content of the image to the centered tombstone. And after some other processing steps I'm trying to find a suitable representation of this image to enable CBIR. So far SIFT and SURF were good detectors but others like MSER fail to produce good results. The problem is the heterogeneity of the data-collection (6.700 picture of different tombstone with different shapes, symbols, ornaments and Hebraic letters).

Does anyone can give me a hint which pattern-recognition approach might be suitable? OCR won't do.

Thank for your help! If you need more picture, take a look around at the institue of German-Jewish history.

From what you said in the comments, it seems that CBIR is not really what you need here.
  I am guessing that you would help the historians more by clustering the images, rather than by doing
retrieval.

SIFT and SURF are a good place to start.  You might also try the FREAK descriptor, which you can compute at the
points detected by SIFT or SURF.

Then the important question is what you do with all those features.  For a problem like this you should consider the bag-of-features approach or the pyramid match kernel.

Edit: you may be right to look for a different interest point detector. SIFT and SURF detectors give you centers of blobs (roughly circular patches of uniform intensity).  Since you are looking at man-made gravestones with Hebrew text on them, corners might work better.  A classic approach is the Harris-Laplacian corner detector.  Note that the original Harris detector is not multi-scale.  For this application you definitely need multiple scales.  Other corner detectors such as FAST and AGAST have been proposed recently.  I have heard good things about AGAST, but I have not tried it.


Autocorrelation function $R_{yy}(t_1,t_2)$? If $x(t)$ is a zero mean stationary Gaussian process and if $y(t)=x^2(t)$,then $\{y(t)\}$ is called a square law detector process. Now i want to find autocorrelation function $R_{yy}(t_1,t_2)$,that is
$$R_{yy}(t_1,t_2)=E({y(t_1)y(t_2)})$$
$$=E(x^2(t_1)x^2(t_2))$$
Now after this step every books give the result like this

  $$R_{yy}(t_1,t_2)=E({x^2(t_1)})E({x^2(t_2)})+2E^2({x(t_2)}X
(t_2))$$

But i don't know how this derivation came. Can anyone help me here please?

Let $\sigma^2$ denote the common variance of the random variables comprising the zero-mean stationary Gaussian process $\{X(t)\}$. We want to find
$$R_{Y}(t,s) = E[Y(t)Y(s)] = E\big[(X(t))^2(X(s))^2\big].$$
Now,
given that $X(t)=x$,
\begin{align}
E\big[(X(t))^2(X(s))^2\mid X(t)=x\big] &=
E\big[x^2(X(s))^2\mid X(t)=x\big]\\
&= x^2E\big[(X(s))^2\mid X(t)=x\big]
\end{align}
where we know (since $X(t)$ and $X(s)$ are jointly Gaussian random variables with covariance $R_X(t-s)$ and correlation coefficient $\rho = \sigma^{-2}R_X(t-s)$) that the conditional
distribution of $X(s)$ given that $X(t) = x$ is a Gaussian distribution
with mean $\rho x$ and variance $\sigma^2(1-\rho^2)$. It follows that
\begin{align}
E\big[(X(t))^2(X(s))^2\mid X(t)=x\big] &= x^2E\big[(X(s))^2\mid X(t)=x\big]\\
&= x^2 \left(\sigma^2(1-\rho^2) + \rho^2x^2\right).\end{align}
Hence, the random variable $E\big[(X(t))^2(X(s))^2\mid X(t)]$
equals $(X(t))^2 \left( \sigma^2(1-\rho^2) + \rho^2(X(t))^2\right)$ and the law of iterated expectation gives

\begin{align}
E\big[(X(t))^2(X(s))^2\big]&=E\bigr[E\big[(X(t))^2(X(s))^2\mid X(t)\big]\bigr]\\
&= E\big[(X(t))^2\left(\sigma^2(1-\rho^2) + \rho^2(X(t))^2\right)\big]\\
&= \sigma^2(1-\rho^2)E[(X(t))^2] + \rho^2 E[(X(t))^4]\\
&= \sigma^4\left((1 - \rho^2) + 3\rho^2\right)\\
&= \sigma^4\left(1 + 2\rho^2\right)\\
&= (\sigma^2)^2 + 2 (\rho\sigma^2)^2\\
&= \left(R_X(0)\right)^2 + 2 \left(R_X(t-s)\right)^2\\
&= E[(X(t))^2]E[(X(s))^2] + 2\left(E[X(t)X(s)]\right)^2.
\end{align}
This differs in subtle ways (not just the use of $t$ and $s$ in place of $t_1$ and $t_2$) from what the OP insists is the right answer as stated in several books that he claims to have read, which answer he has faithfully transcribed into his question above.


OFDM sensitivity to sampling time offset (with half of time interval between samples) In simple QAM modulation symbol synchronization is crucial: the samples need to be taken (or interpolated and resampled to ) exactly corresponding to the time when the symbol was emitted.

This is why symbol syncronization methods like Muller-Muller was developed.  I call this - for the sake of contrasting it with coarse OFDM symbol syncronization (achievable without interpolation and resampling) - fine-symbol-syncronization or smaller-than-time-interval-between-samples syncronizaiton.

Now, the question is, if such fine syncronization is important or not for OFDM. In OFDM there is the cyclic prefix which is much longer than the time interval between samples and there is also the FFT. It is not clear if in the case of OFDM the fine-symbol-syncronization is really neccessary. It might be that the guard interval+FFT makes it unneccessary.

One reason why fine-syncronization should not matter is that the purpose of the guard interval (cyclic prefix) is to make demodulation insensitive to channel response (multi path due to reflections).

Is it important with OFDM to carry out fine, i.e. smaller-than-time-interval-between-samples syncronization ?

EDIT: this book seems to discuss this question in Section "5.2.4 Synchronization Errors", it says that timing offset causes rotation - increasing in strength for sub-carriers further away from the center - need to read it more carefully.

Yes, fine timing synchronization is still very important for OFDM. As you know, an FFT is typically used to demodulate an OFDM signal once it's been properly synchronized and stripped of cyclic prefix. To see what would happen if you aren't finely synchronized in time, recall the time-shifting property of the DFT:

$$
x[n-l] \Leftrightarrow X[k] e^{j\frac{2\pi kl}{N}}
$$

That is, an arbitrary timing offset $l$ (which need not be an integer for this formula to be valid) causes a linear phase shift in each output bin that is proportional to its distance from zero frequency.

So, your subcarriers near the center of the band will not be affected that much, but as you move outward, your demodulated signals will experience an increasing amount of phase rotation. That phase rotation is related to the unknown timing error $l$, and this value will likely vary slightly from symbol period to symbol period (assuming you don't have perfect symbol rate synchronization either). This all adds up to degraded performance.

The solution? You need fine symbol synchronization in OFDM, just like with any other digital modulation format.


Normalization of a signal in MATLAB I have a signal in MATLAB

s = randn(1,1e6);

that is normalized in the following way:

s = s/sqrt(sum(abs(s .^2)) / 1e6);

I would be thankful if anyone could explain me why is the signal normalized in this way!

Assuming that $N$ is the length of your signal $s$, the normalized signal $s_n$ is given by:

$$s_n = \dfrac{s}{\sqrt{\dfrac{\sum_{i=1}^{N}\left|s_i^2\right|}{N}}} $$

The denominator is nothing else than Root Mean Square value of your signal. Thus the code is doing a simple RMS normalization.

You can think of it as a method of normalizing the average of your signal and still allowing for some peaks to clip it (instead of being set to 1). In other words the standard division by the maximum absolute value of your signal will always guarantee that the sample values will be within the range of $[-1;1]$, whereas RMS normalization doesn't.

This method is widely used for audio processing and speech processing.


Why does an anti-symmetric function has zero amplitude at the center of an even length window I am performing FFT on a real odd function and the resultant transform has zero amplitude in the last bin. Essentially if Y= rfft(X), then Y[-1] is always zero. I stumbled on this answer which says

  An anti-symmetric waveform has to be zero at the center of an even length window.

Now, this might be very easy to prove but I cannot think of why this is true. If I use the formula then to the first term should be zero, not $N/2$ th. It's clear to me why it should be real but why should the amplitude be zero for real-valued, odd function. Any help will be appreciated.

A real symmetric function has the following property $$f(t) = f(-t)$$

whereas a real anti-symmetric function has the following property $$f(t) = -f(-t)$$

So at $t=0$ you have $$f(0) = -f(0)$$ which implies that $$f(0) = 0$$

A complex conjugate-symmetric function has the following property $$f(t) = f^*(-t)$$

whereas a complex conjugate-antisymmetric function has the following property $$f(t) = -f^*(-t)$$

So at $t=0$ you have $$f(0) = -f^*(0)$$ which implies that $f(0)$ is real and further $$f(0) = 0$$


DFT of an audio signal stored as a multi-dimensional array I'm using the Python scipy library to get the data of a particular .wav file into array format. Now, I'd like to find the Discrete Fourier Transform of that signal. The formula for the DFT is,

$X_k = \sum_{n=0}^{N-1}x_{n}e^{-2\pi ink/N}$ for a signal $x(t)$ of length $N$.

The data returned from the wavfile.read function, however, appears as a multi-dimensional array; it begins: array([[0, 0],
       [0, 0],
       [0, 0],
       ...,
       [0, 0],
       [0, 0],
       [0, 0]], dtype=int16). I've seen the formula for the 2D Discrete Fourier Transform, but I was under the impression that this was used for image processing. How would I go about finding the DFT of an audio signal represented by a multi-dimensional array?

It looks like you are reading a stereo file, so you have a left and a right channel. You can do a DFT on the left or a DFT of the right or the mono or all of the above. It really depends on your specific application and requirements.


Low pass to low pass transformation coefficient? I am unable to solve this question, 10.10 from GATE IN 2004 (a previous year question paper for an exam targeted at engineering graduates in India.)

So I tried to solve the 10.10 by considering the real frequency shift transformation's general form,

I am not getting any answer close to the options! I wonder what I am doing wrong.

So 1) what is the answer to the question 10.10 shown above?

2) how is this formula for the coefficient derived?

This is an allpass transformation, i.e., the unit circle is mapped to itself. On the unit circle we have

$$e^{-j\omega_0}=\frac{1-\alpha e^{j\hat{\omega}_0}}{e^{j\hat{\omega}_0}-\alpha}\tag{1}$$

For given values of $\omega_0$ and $\hat{\omega}_0$ you can compute $\alpha$ from $(1)$:

$$\begin{align}\alpha&=\frac{1-e^{-j(\omega_0-\hat{\omega}_0)}}{e^{j\hat{\omega}_0}-e^{-j\omega_0}}\\&=\frac{e^{-j(\omega_0-\hat{\omega}_0)/2}}{e^{-j(\omega_0-\hat{\omega}_0)/2}}\cdot \frac{e^{j(\omega_0-\hat{\omega}_0)/2}-e^{-j(\omega_0-\hat{\omega}_0)/2}}{e^{j(\omega_0+\hat{\omega}_0)/2}-e^{-j(\omega_0+\hat{\omega}_0)/2}}\\&=\frac{\sin\left(\frac{\omega_0-\hat{\omega}_0}{2}\right)}{\sin\left(\frac{\omega_0+\hat{\omega}_0}{2}\right)}\tag{2}\end{align}$$

With $\omega_0=2\pi\cdot 60/400$ and $\hat{\omega}_0=2\pi\cdot 120/400$ we get from $(2)$ the value $\alpha=-0.45965$.


Viterbi Decoder without traceback How can I use Viterbi algorithm for a real time scenario where I cannot afford to have delay in the tracing back the best path to the initial state?

A register exchange method implementation of Viterbi decoder avoids trace back and a nominal initial delay (one time delay or boot up delay) is unavoidable. But this method requires more memory resources.


Gaussian signal generation Edit: Could the following be the answer?

Generate WGN-like-signal which is centered around a set dBm value.
Treat that signal like it was a frequency domain representation of an unknown X time-domain signal.
Do some math to get he unit back into Volts/Hz from dBm/Hz (using reference impedance).
Run inverse FFT.

Original Question:

I am looking for some advise on how to generate an array of samples pulled from normal (or Gaussian) distribution.

But, and here's the deal.
I want to control its PSD or ASD "directly", assuming I am using a refernce impedence of 50 Ohms.

I would like to know Ahead-Of-Time (AOT), before signal gets generated, what STD Dev should be applied to achieve a set PSD or at least average spectral level in a given bandwidth/channel.

This value should be set by a user.

Code below should prove I know (more or less) what PSD or ASD is, I am just looking for an algorithm advise, maybe someone did something similar?

Here's how I calculate the PSD or ASD for a generated white-noise-like signal in a standard way using numpy.

import numpy as np
import matplotlib.pyplot as plt
import logging

if __name__ == "__main__":
    logger = logging.getLogger("Main")
    logging.basicConfig(format='%(levelname)s:%(message)s', level=logging.DEBUG)
    DRAW_PLOTS = False
    mu, sigma = 0, 0.001  # mean and std deviation
    N = 1000
    timestamp = 1  # 10 Hz sampling rate

    # This is in Volts and no timestamp
    s = np.random.normal(mu, sigma, N)
    # Because signal is in Volts
    # our result of fft is in Volts/Hz or would be if we say our timestamp for
    # signal `s` is eqault to 1 second
    # This noise_spectrum can be called Amplitude Spectral Density (ASD)
    # If we want PSD we have to go from ASD to PSD:
    # Square the voltage and also maybe apply some reference impedance
    # ... or don't if you want to end up with Voltz (^2) / Hz
    noise_spectrum = np.fft.fft(s)
    # The following is in dBV^2/Hz
    noise_spectrum_db = np.multiply(
        20,
        np.log10(noise_spectrum)
    )
    # The following is in dBW/Hz (divided by 50 Ohm)
    noise_spectrum_dbw = np.multiply(
        20,
        np.log10(noise_spectrum/50)
    )
    freq = np.fft.fftfreq(s.size, timestamp)
    logger.info(
        f"Avg channel power in dBW: {np.real(np.mean(noise_spectrum_dbw))}"
    )
    logger.info(
        f"Avg channel power in dBm: {np.real(np.mean(noise_spectrum_dbm))}"
    )
```

Below is a function which I wrote long back, when I needed to generate AWGN time-domain samples given Noise PSD in dBm/Hz.

AWGN_NOISE() : Generates Additive White Gaussian Noise of PSD power in dBm/Hz

AWGN has Gaussian PDF with 0 mean and $\sigma^{2} = N_{o}/2$

$NoisePSD_{dBm/Hz} = 10.log_{10}(\frac{N_o}{2.BW})$, Why?

Because, Output Noise Power(in dBm) : $$10.log_{10}(\frac{N_o}{2}) = NoisePSD_{dBm/Hz} + 10.log10(BW)$$

Hence,
$$\sigma = \sqrt{\frac{N_o}{2}} = \sqrt{BW.10^{\frac{NoisePSD_{dBm/Hz}}{10}}}$$

So, the python function would be as follows:

def awgn_noise(length, power, Bandwidth):
    sigma = np.sqrt(Bandwidth * 10**(power/10))
    noise = np.random.normal(0, sigma, length)
    return noise

Hope this answers the question.


What Is the Cause for Poor Results with Adaptive Thresholding? I try to run text cleaning process based on the camera captured image.

By using OpenCV adaptive threshold method with parameter (255.0D, AdaptiveThresholdType.MeanC, ThresholdType.Binary, 15, 10.0D), almost all gives satisfactory output but below example gives poor cleaning result (background is not removed well).

(image from http://www.fmwconcepts.com/imagemagick/textcleaner/index.php)

Question 1. I'd like to know which factor cause poor result in this example.

My Question 2 is, which approach can descriminate my attempt(OpenCV) is failed based on the above result? Or, which approach can help to decide whether AdaptiveThreshold will fail or not before running it?

IMO, the pass/fail factor can be a number of background particles but I don't know how to measure it.

I tried Erode method in OpenCV to remove particles, but it gave over whitened result and made too skinny text.

Some guidelines for using Thresholding:

Stretch the image to use the whole Dynamic Range (DR).
Apply some Denoising (Very very gentle).
Median with small radius would be a good idea.
Unless you hand tweak the Threshold, Otsu's Method generally yields good results for this kind of tasks (Text on background).
If one use Adaptive Local Methods (Mean / Gaussian based) the Radius is another hyper parameter to tweak.