Notations to Use in Formulating of Maximum Likelihood Estimation The received
noisy signal $y_n \in \mathbb{R}$ is expressed as:
\begin{align}
y_n = \mathbf{h}^\mathsf{T}\mathbf{u}_n + w_n.
\tag{1}
\end{align}
$\mathbf{h} = [h_0,h_1,\ldots,h_{p-1}]^\mathsf{T} \in  \mathbb{R}^{1 \times p}$ of length $p$ which represents the impulse response of length $p$, and $\mathbf{u}_{n} = [u_{n}, u_{n-1}, u_{n-2},\ldots,u_{n-p+1}]^\mathsf{T}$  having variance $\sigma^2$. $w \sim N(0,\sigma^2_w)$ is the Additive White Gaussian  noise. $\mathbf{y} = [y_0,y_1,\ldots,y_N]$ and $\theta = [\mathbf{h},\sigma^2_w]$
I will be estimating the following terms of the model below using $y$. The estimation of the model is formulated as follows:
 \begin{align}
  z_n = {\mathbf{h}}^\mathsf{T}{\mathbf{u}_n}.
  \tag{2}
 \end{align}

I want to use a variable $z$ to separate what is known and unknown. Therefore, I have used $z$. But, if I use this new variable, the mathematical formulation should have the term $z$.
The condition pdf $\mathbf{y}$ given $\theta$ is given by
$P(\mathbf{y}|\mathbf{\theta}) =\prod_{n=1}^{N}\frac{1}{\sqrt{2\pi \sigma^2_w}}\exp(-\frac{(y_n -z_n)(y_n -z_n)^T}{2\sigma^2_w})$

My confusions are:

1) What notations to use for probability density function is it the one below:

$\mathsf{P}_y(y_n|{\mathbf{u}_n})$ or $\mathsf{P}_z(z_n|{\mathbf{u}_n})$
what goes in the subscript if I want to use $z$?

2) Is $P(\mathbf{y}|\mathbf{\theta})$ correct or $P_y(\mathbf{y}|\mathbf{\theta})$ or $P_z(\mathbf{z}|\mathbf{\theta})$?

3) If I want to use expectation maximization, then would $\xi$ which denotes the complete data set consisting of the input signal and the received noisy signal be written as:
 \begin{align}
 \xi &= \{u_0,\ldots,u_{N-1},y_0,\ldots,y_{N-1}  \} \tag{3}
 \end{align}
or
\begin{align}
 \xi &= \{\hat{u_0},\ldots,\hat{u_{N-1}},z_0,\ldots,z_{N-1}  \} \tag{4}
 \end{align}

4) Would the posterior be written as
\begin{align}
p_{\theta}(y_{1:N}|u_{1:N})p_{\theta}(u_{1:N})  = \prod_{n=1}^{N} \mathsf{P}_u(u_n|{u}_{n+1}) \prod_{n=1}^{N} \mathsf{P}_y(y_n|\mathbf{u}_n).
\tag{5}
\end{align}
 or
\begin{align}
p_{\theta}(z_{1:N}|u_{1:N})p_{\theta}(u_{1:N})  = \prod_{n=1}^{N} \mathsf{P}_u(u_n|{u}_{n+1}) \prod_{n=1}^{N} \mathsf{P}_z(z_n|\mathbf{u}_n).
\tag{6}
\end{align}
Please help what the correct notations are used to represent the pdf, likelihood, log-likelihood if I use the variable $z$ model.

Please point out any other mistakes as well. Thank you.

1) What notations to use for probability density function is it the one below: $\mathsf{P}_y(y_n|{\mathbf{u}_n})$ $\mathsf{P}_z(z_n|{\mathbf{u}_n})$ what goes in the subscript if I want to use $z$?

One way to write the density function is to subscripted upper case letters for the random variables. So $p_{Y_n}$ denotes the density of random variable $Y_n$ and $p_{Z_n}$ is the density of the random variable $Z_n$. Similarly, $p_\mathbf{Y}:=p_{Y_1,\ldots,Y_n}$ denotes the joint density of all the $Y_i$'s.

  2) Is $P(\mathbf{y}|\mathbf{\theta})$ correct or
  $P_y(\mathbf{y}|\mathbf{\theta})$ or
  $P_z(\mathbf{z}|\mathbf{\theta})$?

See my answer to (1) above. In some cases it is ok to drop the subscript if it is understood which random variable the density $p$ belongs to. However, if you have lots of different random variables it is good to explicitly include the subscripts.

  If I want to use expectation maximization, then would $\xi$ which
  denotes the complete data set consisting of the input signal and the
  received noisy signal be written as:
  $\xi = \{u_0,\ldots,u_{N-1},y_0,\ldots,y_{N-1} \}$ or
  $\xi = \{\hat u_0,\ldots,\hat u_{N-1},y_0,\ldots,y_{N-1} \}$

It is not clear what the hats mean.

  4) Would the posterior be written as

You need to clarify the posterior of what conditioned on what. I assume you are interested in estimating the $u$'s so the posterior you are referring to is

\begin{eqnarray}
p_{\mathbf{U}|\mathbf{Y}}(\mathbf u|\mathbf y) &\propto& p_{\mathbf{Y}|\mathbf{U}}(\mathbf y|\mathbf u) p_\mathbf{U}(\mathbf u) \\
&=& p_{Y_1,\ldots,Y_N|U_1,\ldots,U_{N}}(y_1,\ldots, y_N|u_1,\ldots u_N) \; p_{U_1,\ldots,U_{N}}(u_1,\ldots u_N)\\
&=& \left(\prod_{i=1}^N p_{Y_i | U_{i},\ldots U_{1 \wedge (i-p+1)}} (y_i|u_{i},\ldots u_{1\wedge (i-p+1)}) \right) \; p_{U_1,\ldots,U_{N}}(u_1,\ldots u_N).
\end{eqnarray}

That $\wedge$ symbol denotes a min. operation and makes sure that the subscripts don't go below 1. If you have a probabilistic model for the $U_i$'s you can further factor out $p_{U_1,\ldots,U_{N}}$ into a product of conditional densities. For example, if it's a 1st order Markov model with a stationary initial distribution $p_{U_1}(\cdot)$ you can write,
$$
p_{U_1,\ldots,U_{N}}(u_1,\ldots,u_N) = p_{U_1}(u_1)\;\prod_{i=2}^N p_{U_i|U_{i-1}}(u_i|u_{i-1}).
$$

  pdf, likelihood, log likelihood

Finally, I think, you may be interested in writing the "complete data" log-likelihood function i.e. the density of $\xi$ parameterized by $\theta$. This can be written as follows:

\begin{eqnarray}
p_{\xi|\theta}(\xi|\theta) &=& p_{U_1,\ldots ,U_N,Y_1,\ldots ,Y_N|\theta}(u_1,\ldots u_N, y_1,\ldots y_N | \theta) \\
&=& p_{Y_1,\ldots ,Y_N | U_1,\ldots ,U_N,\theta}(y_1,\ldots y_N | u_1,\ldots u_N,  \theta) \; p_{U_1,\ldots ,U_N|\theta}(u_1,\ldots u_N| \theta)
\end{eqnarray}

and simplify this further using your observation model for $Y_i$'s and Markov transition model for $U_i$'s.

It seems that your confusion stems from being unfamiliar with how random variable notation works. I would recommend reviewing a textbook on probability and random processes for electrical engineers.


How to calculate the discrete power of a discrete signal? If I have a time series defined by a series like
$$
x = [a_1, a_2, a_3, ..., a_n]
$$
for time from $t = 1$ to $t = n$

How can I get to the power of the signal in this form:
$$
P_x = [p_1, p_2, p_3, ..., p_n]
$$
Can I use this formula and make it so $N_0 = N_1$

$$
P_x=\frac{1}{N_1-N_0+1}\sum_{n=N_0}^{n=N_1}\left|x(n)\right|^2
$$

Then calculate the $P$ value for every instance in $X(t)$?

Thank you.

The instantaneous power is simply $p(k) = x^2(k)$. The signal energy is given by $W_x = \sum_{k=0}^{k=N} x^2(k) $.


Summation of audio tone Suppose you transmit at a single tone $\omega\textrm{ rad/sec}$ and you have two receiving microphones $i=1,2$ which sample at $r$ samples per second. Assume you have a direct path with coefficients $h_{d,i}$ and a delay of $\tau_{d,i}$ and one reflected component with coefficients $h_{r,i}$ and a delay of $\tau_{r,i}$ at each microphone.

How does the received signal look in analog and digital time domain and also how does the FFT look like?
Do we still perceive only one spike in a particular bin in FFT? Or should you sample appropriately to achieve this?
How does the FFT look if we have two tones?

OK, so suppose your transmitted signal is $s(t)$ and the received signal is $r_i(t)$ in continuous (analog) time .

Then the received signal at the $i^\textrm{th}$ microphone is
$$
r_i(t) = h_{d,i} s(t-\tau_{d,i}) + h_{r,i} s(t-\tau_{r,i})
$$

  How does the received signal look in analog and digital time domain and also how does the FFT look like?

Well if $s(t) \stackrel{\textrm{FT}}{\longleftrightarrow} S(\omega)$ then
$$
r(t) \stackrel{\textrm{FT}}{\longleftrightarrow} R(\omega) = \left [ h_{d,i} e^{-j\omega\tau_{d,i}} + h_{r,i} e^{-j\omega\tau_{r,i}} \right ] S(\omega)
$$
which is simple application of the time shift property.

  Do we still perceive only one spike in a particular bin in FFT? Or should you sample appropriately to achieve this?

Once you start talking about 'bins' they we move to discrete-time. In general, unless you are very careful with sampling, it's almost impossible to sample a single sinusoid and achieve a spike in just two bins (positive and negative frequencies).  Achieving a spike in just one bin requires a complex exponential which is not really possible with real-world sampling unless you use a specific modulation scheme (which you don't mention).

  How does the FFT look if we have two tones?

Again FFT implies discrete time. If you send two tones and sample them, then unless you choose your sampling scheme very cleverly, you will end up with energy in many bins. See this discussion for more details.

  he math of this is what I am seeking. Assume baseband audio such as speech.

@Turbo If the signal is baseband audio, then there will never be just a single tone. Speech is generally modelled as an autoregressive (LPC) model driven by either white noise or a periodic pulse train. See figure 2.40 of this reference (screenshot included below for clarity).


Implementing discrete time systems Discrete time systems work on sampled continuous time signals. We need samples only at discrete time to be fed as input to a discrete time system. Does this mean discrete time systems are clocked so that the values only at discrete time are available to the system? How are these discrete time systems implemented physically?

"Does this mean discrete time systems are clocked so that the values only at discrete time are available to the system?"

yes, but the system can do some mathematics to interpolate values in between the discrete times of the samples.

"How are these discrete time systems implemented physically?"

well, it's called "electronics" and the discipline to learn is "electrical engineering".

the most essential part is a computer of some form.  DSP chips are little computers of a form that are optimized to perform the most common DSP algorithms.  you will likely also need something like an A/D converter and D/A converter, but without more specific information about what the discrete-time system is supposed to do, i wouldn't know that for sure.


Why discrete cosine transform may not maintain locality In paper Convolutional neural network for speech recognition, they say

I didn't understand the highlighted sentence.

In a frequency-like transformation (Fourier, discrete cosine, Walsh), one generally ends up with a sequence of coefficients that account for each frequency component over the support of the basis vector (a segment in time, a 2D patch in images). If you have a frequency component located on one half of the support, and zero the other half, or the converse, you often get the same absolute-valued coefficients. MFCCs fall in that context with:

  the short-term power spectrum of a sound, based on a linear cosine
  transform.

Information about locality can usually be traced in the phase or the sign of the coefficients, but they are difficult to read in complex situations, and are not preserved with the absolute value or the power.


Lowest frequency within pitch detection in polyponic sound I perform some tests where I run cepstrum pitch detection algorithm over polyphonic sound. Result shows that there is almost always occurrence of lowest frequency of tones played simultaneously.
I think this is because the lowest frequency has the longest period. But I do not know if it is true or why it is so.
Can you somebody explain this or why there is lowest frequency? Or suggest some publications?

The lowest pitch in a chord is likely to have a large number of its possible overtone frequencies match the actual overtones of other notes in a harmonious chord, leading to a good cepstral magnitude, which looks for spectral periodicity (overtone series's).

However a cepstrum might also find more overtone sequence periodicities at pitches lower than the lowest pitch in a chord.


Approximate the magnitude response of an analog filter with a digital filter, starting from its gain expression I have an analog filter with its frequency response curve in dB described by the following expression:
$$
N_{dB}=20log_{10}\omega t_1 \sqrt{\frac{1+(\omega t_2)^2}{1+(\omega t_1)^2}}
$$
This expression is derived from the series connection of two lowpass filters each associated with the following RC circuit:

where, for each circuit, the time constant $t_i=RC$ and of course $\omega = 2\pi f$, where $f$ is the frequency (this is actually the equalization curve for magnetic tape recording/playback, see Annex B, page 14 of this document).

I would like to obtain an approximation of this frequency response using a digital filter. I don't know if there is a method to exploit our knowledge of the analog frequency response or if I should design the filter myself from scratch.
The end goal is to obtain the impulse response and save it as a a .wav file (I know how to do this last part). I just took a basic DSP course at my Uni but we didn't work with analog filters so I am a little bit lost.

Did you try to plot this curve in matlab.
I am not sure whether it is any standard high,low or band pass filter.
So, to calculate the impulse response there might be two ways.
1. take the inversefft of the NdB and then do the convolution between your impulse input and the inversefft signal of filter usinc conv(x,y) command.
2. Or convert the impulse input into frequency domain by fft transform in matlab and then just multiply the gain and fft of input impulse and then convert them in time domain.
Can you do it using above method?
Do you want a matlab code for this method?


Peak to peak amplitude of sum of sinusoidals (harmonic frequencies) I would like to estimate the peak-to-peak amplitude of a periodic signal whose frequency components are known. This is, I have the frequency spectrum (a peak in the fundamental frequency and other peaks in its harmonics) and I would like to compute the peak-to-peak amplitude. Is that possible without reconstructing the signal in the time domain and then detecting peaks?

Suppose that we have a continuous-time periodic signal
$\displaystyle s(t) = a_0 + \sum_{n=1}^N a_n \cos(n\omega_0 t + \theta_n)$.

What does this mean? Do we have a  trace of $s(t)$ on some recorder chart
and the trace looks periodic?  If we did, the question to be solved would be
trivial since we could simply measure the maximum and minimum values of $s(t)$ on
the chart.  So the problem of interest is:

  Given the values of
  $\omega_0, a_0, a_1, a_2, \cdots, a_N, \theta_1, \theta_2, \cdots, \theta_N$,
  find the values of
  $$M_\max = \max s(t)\quad \operatorname{and} \quad M_\min =  \min s(t).$$

The OP wants to find the peak-to-peak amplitude of $s(t)$ and this is, of
course, just $M_\max-M_\min$.  It is worth noting that since $s(t)$ is periodic
with period $T = \frac{2\pi}{\omega_0}$, the maximum value and minimum value
occur (at least once) in each interval of length $T$ on the time axis.  The key
issue here is since we know the $a_i$ and the $\theta_i$ etc., we can write
down the formula
$$s(t) = a_0 + \sum_{n=1}^N a_n \cos(n\omega_0 t + \theta_n)\tag{1}$$
(whether this constitutes "reconstruction" of the time-domain signal
is a question that I will leave to others to answer) but the
standard method for finding the maxima and minima of $s(t)$
requires us to find solutions to the nonlinear equation $(2)$ below
(that is, values of $t$ for which $(2)$ holds):
$$\frac{ds(t)}{dt} = -\sum_{n=1}^N a_n\cdot n\omega_0 \sin(n\omega_0 t + \theta_n) = 0.
\tag{2}$$
Once we have found the values $t_1, t_2, \ldots, t_k, \ldots$ for which
$\displaystyle \sum_{n=1}^N a_n\cdot n\omega_0 \sin(n\omega_0 t_k + \theta_n) = 0$,
we can simply calculate $s(t_1), s(t_2), \ldots, s(t_k), \ldots$ and search
through this list of numbers (it is not necessary to sort the list
as recommended in hotpaw2's answer) to find $M_\max$ and $M_\min$.  But,
ignoring the difficulty of finding the $t_k$'s makes hotpaw's answer not
particularly useful in solving the OP's problem.  The answer by ethereal
is even worse in this regard since all it boils down to is the assertion
that (for $\alpha > 0$)
the maximum and minimum values of $\alpha\cdot s(t)$ are $\alpha\cdot M_\max$
and $\alpha\cdot M_\min$ respectively. Finding the maximum and minimum
of $s(t)$ from knowledge of its Fourier series is a nontrivial task, not at
all as easy as it is made out to be by ethereal or by hotpaw2.

A related problem has received a lot of research attention in the past
thirty years or so.  For a periodic signal such as $s(t)$, the average power
in the signal is readily computed as
$$\bar{P} = \frac{1}{T}\int_0^T |s(t)|^2\,\mathrm dt = a_0^2 + \frac{1}{2}\sum_{i=1}^N |a_i|^2.$$
On the other hand, the peak power is $P_\max = \max\{M_\max^2, M_\min^2\}$,
and is also of interest, especially to system and power amplifier designers,
and the ratio $\frac{P_\max}{\bar{P}}$, aptly named the _peak-to-average-power
ratio (PAPR) has received much attention. It can be calculated exactly for
a given $s(t)$ but only with great computational effort, and so a lot of effort
has gone into finding bounds on the PAPR.  Some of these could help in getting
bounds for the OP's problem.
only with a great deal of effort.


Normalization of a signal with respect to another signal Imagine that I gave one signal (call it master signal) with respect to which I would like to normalize all the others. How should I treat it mathematically? Any ideas?

I have had the same problem. My signals are coming from the MRI machine. In my case, the scaling is one of my options; in a way, I divided the signal A and B to the maximum of signal A, but I could not use z-score approach because the information has been embedded in the amplitude. If your signal is a non-amplitude modulated like a speech signal, you can use a z-score.


How will the signal $\sum_{n=-\infty}^\infty \Delta (t-n)u_{-1}(t-n)$ look like? How will the signal $\sum_{n=-\infty}^\infty \Delta (t-n)u_{-1}(t-n)$ look like?

I understand that $\sum_{n=-\infty}^\infty \Delta (t-n)$ is a periodic signal with period $1$, where the triangular pulse $\Delta(t)$ is shifted by $1$ units for each value of $n$.

where,
$$\Delta(t) =
\begin{cases}
1-2\vert t\rvert,& \lvert t\rvert<\frac{1}{2} \\
0,&\lvert t\rvert>\frac{1}{2}
\end{cases}$$
From this I infer that $\sum_{n=-\infty}^\infty \Delta (t-n)u(t-n)$ would probably be a signal in which the RIGHT HALF of $\Delta (t)$ will be shifted by $1$ unit for every value of $n$. But here I have a doubt that $u(t-n)$ will be added with the previous values of $u(t-n)$ i.e $u(t-(n-1))$ and will have a staircase effect. Thus the signal will no longer be periodic.(Which is a cause of concern as we have been asked to find the Fourier series of the signal, which means the signal should be periodic)
I also don't understand what the subscript $-1$ in $u_{-1}(t-n)$ means.

Please help me visualize this signal(It would be highly helpful if you would  draw the graph as well to help me visualize the signal)

Initially, don't worry about the sum. Just try to figure out how the signal $\Delta(t)u_{-1}(t)$ looks like. As soon as you know that, you just need to shift it by integer values of $t$, and add them all up. If the shifted versions do not overlap, you simply get copies of the basic signal shifted to integer values of $t$. From the solution in the solution manual (as given in a comment), the basic signal is the right half of a triangle, and the period equals $1$. The symbol $u_{-1}(t)$ refers to the unit step function, which is more frequently referred to as $u(t)$ (see also this question and its answer(s)).


Measuring similarity between two images I have some images (generated from some processing). Some of those images have approximately the same view of the scene with illumination changes. But there are other images which may show the view but with a different viewpoint. I am trying to remove the second set of images mentioned above. I've tried using correlation coefficient as a measure of finding similarity. But most of false image pairs seem to have high correlation coefficient. Then I've tried using mean-square error. On setting  the value at 1500(guess), results have improved. Note all the images have similar size. This setting of 1500 changes for different images. So how do I choose this quantity? Or is there a better way to compare such images?
(The images are gray-scale)

Edit:
I am sorry I was not specific. I have found SURF correspondences between the images and used the transformed images (obtained by projective transform) for matching between the transformed image and part of scene image(the ROI was extracted by transforming coordinates by obtained matrix above).
I have used SSIM as one answer suggested but that metric seems to punish illumination changes which is not what I want.

I don't quite know what you mean by different viewpoints, however, there are a lot of different similarity metrics you could use to compare images. You could try:

extracting features and comparing the features Like Scale-invariant feature transform (SIFT) features
you could try dimensionality reduction and comparing that representation, like taking the discrete cosine transform (DCT).
I know you said you tried mean squared error but you could also try comparing the structural similarity (SSIM).

Edit:

Now that I better understand your question, I believe I found a resource that could be of assistance. Here is a very aptly named paper that is titled "A Novel Algorithm for View and Illumination
Invariant Image Matching".


Thesis Topics for DSP using Matlab/Simulink/LabVIEW I am searching for a good thesis topic on DSP with the use Matlab/Simulink/LabVIEW for my master, however, when I contact my department (we only have 1 DSP professor), he asks me to try to propose one but I cannot think of any.

I wonder if someone can help me with some ideas

If I were you, I'd go for something related to Compressed Sensing. Even though my Ph.D. is not on this subject, I've come across it couple of months ago and it looks quite interesting, and it's relatively a new topic too (correct me if i'm wrong). Violating the Nyquist Sampling rate has never been that scientific and fun I suppose. Give it a shot, read a couple of papers (there are really good review papers out there that can give you a good introduction, you might want to visit this website for a good start).

Good luck and try to have fun while doing it, that's a must!


Power Spectrum in dBFs I'm trying to get the spectrum of my signal in dBFS units. I'm not sure what is the refernce that I should devide by. I computed the power spectrum as follows:

1) FFT normalized to the number of samples

p = abs(fft(signal))/(No_of_samps/2);

2) and then I calculated the power spectrum for half of the samples (one-sided)

p = p(1:Nsamps/2).^2;

3) I  tried to represent it in dBFS as follows:

p_dBFs = 10*log(abs(p)/(vref/2)^2);

Where vref is my reference p-p voltage in the ADC.

I'm just not sure whether I should devide by the reference voltage to get values in dBFS, like when I apply a FS signal I would like to see the highest peak in the spectrum at exactly 0 dBFs.

Note: The time-domain signal is also multiplied by a hamming window before computing the FFT.

Any help is appreciated.

Sorry for the long post. Thank you guys very much.

Hello again guys after a lot of thought, I finally broke this down and I would like to share it with you hoping it'll help.

So what I did was:

1- I used matlab's pwelch to get the power spectrum like this

[AutoPowerSpectrum ,F] = pwelch(signal, Nsamps, 0, Nsamps, fs, 'power');

you can read about the different parameters for pwelch in matlab's documentation right HERE.

2- I calculated Vpeak squared, because as far as I know pwelch produces Vrms values.

AutoPowerSpectrum = 2.*AutoPowerSpectrum;

3- got the Vpeak by squaring the above equation

AutoPowerSpectrum = sqrt(AutoPowerSpectrum);

4- converted it to dBFs

p_dBFs = 20*(log10(AutoPowerSpectrum./(vPerLevel*2^(N-1))));

where: vPerLevel is the voltage corresponding to each level in my ADC. It is calculated like this Vref/2^N.

and N is the resolution of the ADC.

Finally used plot() to display it. This solution worked with me and gave me the signal in dBFS, i'm opened to edits and suggestions.

Hope this helped.

Cheers.


Complex Samples (IQ) - Baseband Filtering We are currently analyzing a large set of IQ samples in a desktop application and we are interested in implementing many different bandpass filters dynamically.
We realized, that working with scipy offers no suppport for complex bandpass filtering.
We have already checked following link which suggests a solution to the problem when it is decided to approach the problem with complex filters:
How to implement bandpass filter on complex valued signal?
We wonder, why that is even necessary, since the data could be transformed to a real format. Instead of $[-f_s/2, f_s/2]$ the range is from $[0,fs]$ (mirrored about the $0\,Hz$ point). This way already implemented filtering tools could be used.

A real signal is complex conjugate symmetric, so any filter over $-fs/2$ to $+fs/2$ is really only unique over the range from $0$ to $fs/2$. This would be a primary motivation for working with a fully complex signal as the unique range is then truly extended over the full $-fs/2$ to $+fs/2$ range.
Another more dominant reason to implement a complex filter besides the bandwidth requirement outlined above (where the choice is really sample twice as much or carry two datapaths; similar complexity and certainly you could map either to be similar if that was the only goal) is when the passband itself is not symmetric: specifically as a bandpass filter when the positive frequency passband is not equal and magnitude and conjugate in phase to the negative frequency, or as a baseband filter- the passband shape itself is asymmetric or perhaps completely one sided. Another dominant reason is when real signal conditions would create images very close to the signal where further filtering would be difficult, such as if conditions require operating very close to the Nyquist boundary.  If the OP doesn't have these conditions, I don't see a strong motivation for complex filtering, which requires 4 real multipliers for every complex multiplication operation (assuming complex inputs and outputs).
Channel equalization is another example where a full complex filter is often necessary (as an implementation on a baseband signal) but not needed for a passband signal.


Are there any good open-source (well, free) handwriting OCR programs? The title asks it all. I have before me a data entry task that I'm not too enthusiastic about:  50-100 pages of handwritten sign out/sign in logs.

The format of the logs may help. Pages are partitioned into clearly delineated rows and columns (13r x 6c with an extra typed header row). Further helping me out is that three columns are date/time related (date, time out, time in). Also, the data in two of the columns (resource and name) are more or less enumerated, so that, for example, the name "Smith" might appear time and again in the name column, each time with the same handwriting. The last column, "Notes", is free form, but if I could automate the previous 6 columns, I wouldn't mind entering Notes by hand.

Any suggestions? (Besides 'start typing'.)

PS If there's a better SE site to ask this, let me know, I'll ask there.

tesseract is probbaly the best and most widespread OCR library.

It has been tested with handwriting and isn't too bad - although handwriting isn't an easy thing to read. http://arxiv.org/ftp/arxiv/papers/1003/1003.5893.pdf


Scatter plot: calculate box where 80 % of the points are I've a scatter plot with x- and y-axis.

Now I want to calculate a rectangle which displays where 80% of the points in this plot are.

I'm not very familiar with statistics, so I don't know the method I have to use. Read something about percentiles. Is that correct?

I know I have to calculate where 80% are on the x-axis and also where 80% are on the y-axis and then I got my rectangle corners, but how do I calculate the 80%?

Sample data:

Your question seems to imply that there will always be just a single cluster of points.  If this is the case, I would go with a fairly simple heuristic like the following-

Calculate the mean and variance of the point positions in both the x and y directions.  The mean is the center point of your box, and you can use the ratio of the x/y variances to determine in which directions the box should grow faster.  For instance, if the x variance is "1" and the y variance is "2", I would have the box expand twice as fast in the y direction.

Now simply expand the box until you have at least 80% of the points contained.

If there is ever more than one cluster of points you may need to use a more sophisticated solution that uses clustering algorithms to identify the clusters and a box for each cluster.


How to add complex WGN to complex damped exponential and compute SNR? I am new to the Signal processing domain, So my question might sound basic. I am trying to reproduce the results from the paper Decimative subspace-based parameter estimation techniques by
Geert Morren, Philippe Lemmerling, Sabine Van Huffel.
I would like to first know what is the unit of the model given here, where $a_k$ amplitudes, $\phi_k$ phases, $\alpha_k$ damping factors, $f_k$ - frequencies in Hz and $f_{sample}$ - sampling frequency.
\begin{equation}
    x(t) = \sum_{k=1}^{K} (a_k e^{(j\phi_k)})(e^{\{(j2\pi f_k - \alpha_k)/f_{sample}\}t}) + Noise(t), \quad t = 0,1,...,N-1
\end{equation}
I am adding circularly symmetric WGN to $x(t)$ using the wgn function from MATLAB (shown below in code), which is in dBW unit.  I am not sure whether this is correct, and I am computing SNR also using a snr MATLAB function.  However, I want to use the variance of noise and scaling of noise to determine SNR.  Could someone explain me clearly how to do this?
f       = [0.02, 0.0205];
a       = [1, 1];
alpha   = [0, 0];
phi     = [0, 0];
K       = length(f);
del_t   = 1;
N       = 1000;
t       = (0:del_t:N-1)';
size_t  = length(t);
f_sample= 1/del_t;
Nruns   = 40;
x       = zeros(size_t,1);
SNR     = zeros(1,Nruns);

for iter = 1:Nruns
    noise = wgn(size_t,1,randi(20),'complex');
    for n = drange(1:size_t)
        for k = 1:K
            x(n) = x(n) + a(k)*exp(phi(k)*1j)*(exp((2*pi*f(k)*1j - alpha(k))/f_sample)).^(n-1);
        end
        x(n) = x(n) + noise(n);
    end
    SNR(iter) = abs(snr((x-noise),noise));
end

One suggestion is to leverage the vector processing in Matlab and to not add the noise in a for loop, but given the vector of your signal of length $N$ you can simply create the white Gaussian Noise using:
noise = std * sqrt(2)/2 * randn(N,1) + std * sqrt(2)/2* 1j*randn(N,1)
(or randn(1,N) if you prefer row vectors instead of column vectors)
Where std represents the total standard deviation of the noise.
Then the noise can be added to your signal simply using
y = x + noise
Next be aware that the total power as a white noise and the square of the standard deviation std in whatever units are chosen is distributed evenly over the entire frequency range given by -$f_s/2$ to +$f_s/2$ where $f_s$ is the sampling rate. SNR of interest is typically that portion of the noise that shares the same spectrum as the signal. (So if the signal only occupied 1/4 of the total spectrum as sampled, then the noise would be 1/4 or 6 dB lower than that predicted by the total variance of noise).


What exactly do we do in non-rigid registration? I found many readings on Non rigid registration, also saw some videos but didn't find a clear explanation anywhere.

What I understood is the following

Rigid Registration: It is the alignment of two point clouds without changing the shape (ie only by rotation and translation). ICP algorithm is used. The closest point in the other point cloud is assumed to be the corresponding point.

Correspondence detection is not necessary for all the points.
All the points are applied the same transformation.
Iterative

Non-rigid registration: In this case, the shape of the point cloud might change. But what exactly are we doing? Are we finding the corresponding points using Spin Images, and then applying separate transformations (translations in this case) for every corresponding point?

Correspondence detection has to be done for all the points.
All the points are applied their separate transformations (might be same).
Not iterative: Just apply a translation

Piecewise rigid registration: We divide the non-rigid point cloud into many rigid parts and then apply separate rigid transformations per part. This lies between rigid and non-rigid registration.

Correspondence detection is done for only some of the points.
Same transformation is applied to those points in a particular rigid part, but different parts might have different transformations.
Iterative

Please correct me if I'm wrong somewhere.

Also, is there any role of SIFT in correspondence detection. I've seen some links for SIFT in 3D.

I have found this survey paper very useful:
http://ralph.cs.cf.ac.uk/papers/Geometry/Registration.pdf

It gives a categorized overview of the field. I think section 3 and 4 are the most relevant ones for you.


Is there a name for $f(x) = |x|^n \cdot \operatorname{sgn}(x)$? I'm looking for a name for the following function

$$f(x) = |x|^n  \cdot \operatorname{sgn}(x)$$

which is the exponent function, but with the sign stripped from the input and reattached to the output.

It seems useful enough for DSP that it ought to have a name, but Google's not much help with this sort of thing. Has anyone seen this function in any DSP software? What was it named?

There are many instances (eg Handbook of Computational and Numerical Methods in Finance, p. 42, Eq. 2.25) where it is simply called the "signed power"  "signed power-$p$ function" or "signed power law".

Some instances are more common. The term "signed square-root" seems to be used, see for instance in this Simulink documentation:

  Square root of the absolute value of the input, multiplied by the sign of the input.

with the signed square as an inverse. It is used in the realm  of statistical distributions, contrast enhancement, or on data features (like SURF) as a compander, close in spirit to the A-law or the mu-law.

For parabolas (and even powers in general), the rewriting into $x\to x.|x|$ could be called a split parabola (see Selected solutions, HW 1, due January 22, 2009, an instance of split functions).


Removing noisy lines from image - opencv - python I want to process different type of  image to extract actual text from noisy image .I am using open cv to do this job it is working fine but the problem is I need to set different lower threshold to convert noisy image to readable text.Is it possible to set single threshold which removes noise from all images by doing some pre processing ?

Main goal : I am trying to remove that wavy horizontal line and make the character clear to read

Code used:

import cv2

# Load an color image in grayscale
img = cv2.imread('it_captcha3.jpg',0)
ret, thresh_img = cv2.threshold(img, 180, 255, cv2.THRESH_BINARY_INV)
cv2.imshow('grey image',thresh_img)
cv2.imwrite("result11.jpg", thresh_img)
cv2.waitKey(0)
cv2.destroyAllWindows()

Captcha1 :

Processed image : threshold used -> low: 180 high ->255,

CAPTCHA : 2

PROCESSED IMAGE: threshold used -> low: 200 high ->255

captcha : 3

processed : low -> 165 high : 255

first of all remember that there is no single solution for all kind of noise and all kind of images. that being said i can think of two solution.
first is using Otsu thresholding:

ret,thresh_img = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY_INV|cv2.THRESH_OTSU)

this will try to guess a good threshold for the image being used.

the other solution would be using "close" morphology transform after thresholding. it will first dilate and then erode the image and using a good kernel( it is also called structuring element ) you can remove the line(although you may remove some of the useful pixels too!!)

kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (4,8))
morph_img = cv2.morphologyEx(thresh_img, cv2.MORPH_CLOSE, kernel)

above is the result of using both otsu threshold and morphology close


Applying a Liftering function to an audio signal Is it necessary to perform the liftering function before applying MFCC to an audio signal?

No, liftering is never applied before computing MFCC because liftering is defined as a windowing operation in the cepstral domain. So you need cepstral coefficients in order to be able to apply liftering.

The effect of liftering is to smoothen the corresponding log magnitude spectrum, which can result in more robust recognition / classification results.

Also take a look at this question, and the references in its answer.


how multi track recording is performed? I have a doubt regarding multi track/multi layer recording of audio signals. in recording studio each instrument or chorus,playback will use a independent layer/track..so after recording a all tracks ..how do they mix all of those into a single track ..any mixing scheme?..or it is the simple time division multiplexing??

I am asking this because while I read any mp3 or wav file in Matlab..I can see only a one channel or 2 channel signal...

any help would be appreciated...

Most of recordings nowadays are being done in stereo - two channels. So when you do a recording in studio, then most of tracks are mono (i.e. microphones, guitars). The main job of sound engineer is to place them in panorama (Left - Centre - Right) to create a nice spatial feeling. For this purpose you have pan knob to set this mono track to the left or right speaker (make it more louder or quieter in either of them). Many tracks are useful when you want to apply different filtering to different instruments or add some reverb, change the volume, etc. In the end all waveforms are being summed and voila!

The main track is obviously a stereo track based on your adjustments of track's pan - that's why you have two channels loaded into your MATLAB.


Getting Frequencies Corresponding to Peaks In FFT Plot - MATLAB So I'm currently working on a project where I need to find the frequencies for specific peaks in my FFT plot of an audio file. However, I'm not sure how to actually get the frequencies that correspond to the peaks I'm interested in.

This is my code: https://pastebin.com/zZW8B5Th

And this is a picture of my FFT plot: https://imgur.com/0ajtful

Essentially, I'm trying to construct a Fourier Series from the signal. Although I know how to find the coefficients of the sine/cosine terms, I'm just a little confused on how to find the frequency.

Thanks everyone!

MATLAB's indexing begins at 1, instead of 0, as most programming languages use.  This causes a little wrinkle in interpreting the fft results.  Index 1 corresponds to bin 0, index 2 to bin 1, etc.  So to get the bin number, subtract 1 from the index.

The bin number corresponds to the frequency in units of cycles per frame, where frame is the section of the signal your are analysing.  Usually, the variable $N$ is used to denote the number of samples in the frame.  You also need to know the sampling rate, usually denoted by $f_s$ which usually has units of samples per second.  CD quality audio has a sampling frequency of 44,100 samples per second.

To calculate the frequency in hertz (Hz) of a bin number, usually denoted with a $k$, do the following calculation:

$$ f_{HZ} = k \cdot \frac{f_s}{N} $$

The units for the equation goes like this:

$$ \frac{cycles}{second} =  \frac{cycles}{frame}   \cdot  \frac{ \frac{samples}{second} }{ \frac{samples}{frame} } $$

Hope this helps,

Ced


Solution to Cross-correlation problem of 2 audio signals There are 2 audio signals, X & Y.
X is of 15 second duration. The duration of Y is not given.
Y is a part of signal X (ie a small part of audio X is extracted & made into audio Y)
I have to find the starting & ending location, in seconds, in the signal X from where signal Y is derived.
I was told to NOT use built-in MATLAB functions like conv, xcorr etc.

So I wrote a MATLAB code which compares sample by sample of X & Y starting from the 1st sample.

clear;
[A1,FS1,NBITS1]=wavread('x.wav');
[A2,FS2,NBITS2]=wavread('y.wav');

for i=1:(length(A1)-length(A2)),
    if(A1(i:i+length(A2)-1) == A2)
        found = 1;
        time1 = i/FS1;
        time2 = (i+length(A2))/FS1;
        break;
    else
        found=0;
    end
end

A1 has 120000 samples, FS1 = 8000, NBITS1 = 16
A2 has 20000 samples, FS2 = 8000, NBITS2 = 16
Therefore,
120000/8000 = 15 sec(length of audio X)
20000/8000 = 2.5 sec(length of audio Y)

time1(starting time) = 6.5001 sec
time2(ending time) = 9.0001 sec

Just to make sure the code was correct, I plotted the impulse response of X & Y using fvtool & found that the impulse response of Y matches that of X between 6.5 to 9 seconds.

Now I want to verify the result using built-in MATLAB functions, namely xcorr
So I used

A3 = xcorr(A1,A2);
plot(A3);

A3 is 239999 samples long. Then using

[max idx] = max(A3)

Using the above index location I found the location of maximum correlation value to be at sample 172000. which gives the midpoint at 172000/8000=21.5sec(clearly wrong). However, the location of the maximum correlation value must come to 7.7501 samples. How do I obtain that from the above value?

Having looked at how xcorr is defined I think I know what is happening.

There are two issues:

The range of xcorr is $2N-1$ (N = length of signal X). The initial values are for negative shifts and the second half is positive shifts. Therefore subtract 120000 to get the correct shift.
xcorr pads Y with zeros to the same length as X. The first 20000 points will be sample Y and the rest of the data zeros. Therefore the value of the maximum correlation will be the start point, not the mid point of the signal (as I originally thought). To prove this think what would happen if Y was the beginning of signal X.

If we now do this on your data we get a start point of (172000-12000)/8000=6.5seconds as expected. The 0.0001 sec comes from how you've defined the first point in time.

However, this is a pretty bad way to do it as any noise in the signal Y will cause your algorithm to fail. I suspect the aim of this task is to manually right code to calculate the correlation function in a similar way to corr does and then find the maximum.


How to determine ALPHA, smoothing constant, of a LPF? I am implementing a simple LPF in Java.

But I have problems in choosing the value of ALPHA.

public class LowPassFilter {

    /*
     * Time smoothing constant for low-pass filter 0 ≤ α ≤ 1 ; a smaller value
     * basically means more smoothing See:
     * http://en.wikipedia.org/wiki/Low-pass_filter#Discrete-time_realization
     */
    private static final float ALPHA = 0.2f;

    private LowPassFilter() {
    }

    /**
     * Filter the given input against the previous values and return a low-pass
     * filtered result.
     *
     * @param input
     *            float array to smooth.
     * @param prev
     *            float array representing the previous values.
     * @return float array smoothed with a low-pass filter.
     */
    public static float[] filter(float[] input, float[] prev) {

        if (input == null || prev == null)
            throw new NullPointerException("input and prev float arrays must be non-NULL");
        if (input.length != prev.length)
            throw new IllegalArgumentException("input and prev must be the same length");

        for (int i = 0; i < input.length; i++) {
            prev[i] = prev[i] + ALPHA * (input[i] - prev[i]);
        }

        return prev;
    }
}

ALPHA is defined to be dT/(dT+RC), where dT is the event delivering rate and RC is the time-constant of LPF.

To find the value of ALPHA, I need to know the values of dT and RC.

My LPF should have a cut-off frequency of 4Hz and the sampling rate is 350Hz.
Anybody can help determine the proper ALPHA value?

Thanks in advance!

All formulas I'm using come from the wikipedia page you cite in your code.

First, calculate RC. We have,

fc = 1/(2*pi*RC)

so

RC = 1/(2*pi*fc)

Now, ALPHA is simply:

ALPHA = dt / ( dt + RC )

where dt is 1/samplerate.

For your values (samplerate=350, fc=4) I got

ALPHA = 0.0670

Caveat: I tried this on a few other values and the shape of the filter was what I expected, but the exact 3dB cutoff frequency was never quite where it should have been. Not sure if this indicates something wrong, or that's just the nature of the beast.


Chunks to Symbols in Gnuradio 3.8 I am trying to follow this tutorial to build a flow graph with GnuRadio 3.8. However, I am unsure what settings should I use for Chunks to Symbols and get some errors when I execute the flow graph:

Traceback (most recent call last):
  File "/home/user/top_block.py", line 186, in <module>
    main()
  File "/home/user/top_block.py", line 164, in main
    tb = top_block_cls()
  File "/home/user/top_block.py", line 129, in __init__
    self.digital_chunks_to_symbols_xx_0 = digital.chunks_to_symbols_ic(1+1j, 1-1j, -1+1j, -1-1j, 1)
TypeError: make() takes from 1 to 2 positional arguments but 5 were given

>>> Done (return code 1)

My flow graph is:

My settings for Chunks to Symbols are as follows:

Input type: int

Output type: complex

Symbol table: 1+1j, 1-1j, -1+1j, -1-1j

Dimension: 1

Num ports: 1

Am I doing something wrong here? Any help is appreciated.

Symbol table: $1+1j, 1-1j, -1+1j, -1-1j$

The Symbol table needs to be a python list/array - [1+1j, 1-1j, -1+1j, -1-1j]


Why should one choose the maximum as reference for power-to-decibel conversion? I have seen that one converts the amplitude to decibels through:
$\text{dB(S)}=10 \text{log}_{10} \big(\text{S/ref}\big)$, where $\text{S}$ is the output of a STFT and ref the reference value for the logarithmic scale, i.e. the value that should be mapped to decibel, right?

However, I've seen implementations with $\text{ref = max(S)}$
and I do not understand why. Can someone clarify on that?

dB is a power ratio, so when we see units in dB we are seeing units on a relative scale. The reason that a 0 dB reference is so common is because this is simply normalizing the number scale to 1. $10Log(1) = 0$ dB

This is similar to using 100% (1) instead of absolute numbers.


Remove Noise from Physiological Signals I have collected a set of physiological signals through emg sensors.
I will link at the end of the post an example (10 sec.) of a signal collected on zygomatic muscle.

When I analyzed these data, I have seen that they cointain a great ammount of noise (as you can see by the plot of the FFT.

Plot of the emg signal

FFT of the emg signal

As you can see, there are several peaks into the FFT. My idea is to delete them using a notch filter on the frequencies of these peaks. Anyway, I am no able to figure out how select automatically these peaks. Do you have suggestions?

Furthermore, are there any better methods to clean the signal?

NB: Each signals have the noises in different frequencies. The data was collected usign a battery and avoiding the current interference.

Thanks :)

File .mat with the data example: https://ufile.io/zk8rf
File .fig with the plot: https://ufile.io/g8owd
File .fig with the fft: https://ufile.io/02m6g

This is a tougher problem than it appears at the first glance. This may be confusing at first, but read carefully, and I hope that it will give a good idea what to do. Here is what is my experience and how I resolved a similar problem:

apply a very good band-pass filter to get the frequency range of interest.
apply a very tight DFT to get a good resolution in frequency domain. You may wish to pad the original signal with zeros if it is not long enough.
determine the bandwidth of the notch filter that you want to use. Be very careful with notch filters, since they have some inconvenient characteristics.
make the running window sum of the DFT results using the width of your notch filter.
set the threshold to decide what is a peak in the results.
the running windows sum should give you values where the peaks are.
apply the notch filter designed for the frequency in mind. If the frequencies are changing, and the filtering has to be done in a real time, adaptive filter has to be applied, but that is yet another question.


Sharp sample transition in OFDM signal I am implementing OFDM PHY for Visible Light Communication in FPGA. (I am quite new to communication system and signal processing, so please understand my lack of background knowledge.)

The FFT size is 64 and it use Hermitian symmetry. The preamble symbol (64 sample) is added in front of the data symbol for synchronization, so for one OFDM symbol consists of 128 time domain samples. This sample is sent to DAC every 200ns (DAC clock is 5 MHz). The DAC maximum update rate is 125 MSPS and the output is connected to the LPF that has 50 MHz bandwidth (-3dB). The DAC output is shown in the following figure.

As shown in the figure, the OFDM signal is ringing on every transition. I think this ringing is caused by sharp transition that is filtered by LPF, so the high frequency harmonics is removed. In the receiver, the ADC sampling clock is 125 MHz, so I have to undersampling it by a factor of 25. The undersampling process is done by pick up one sample in every 25 samples. The OFDM signal after undersampling is shown in the following figure.

As shown in the figure, I think because of the ringing, and the ringing signal is picked up in the undersampling process, then the signal amplitude is decreased or increased.

Question 1: Is the DAC implementation correct, or I need to add oversampling and pulse shaping before DAC? so the output is not ringing?

Question 2: Is there any solution how to remove this ringing in the receiver?

You're running a DAC at 5 MS/s and have a reconstruction filter of 50 MHz bandwidth.

There's your error: You're not doing proper reconstruction filtering! Reduce the filter bandwidth to reflect your sample rate.

You can of course also interpolate your signal 10× and set your DAC clock to 50 MS/s. Interpolating an OFDM signal is potentially very easy: Just add enough unused carriers to the side (note: you might need to adjust your CP length and synchronization algorithms).


Determining the average probability of error I'm trying to solve the following problem:

In a binary PAM system, the input to the detector is
$$y_m = a_m+n_m+i_m$$
where $a_m = \pm1$ is the desired signal, $n_m$ is a zero-mean Gaussian random variable
with variance $\sigma_n^2$ and $i_m$ represents the ISI due to channel distortion. The ISI term
is a random variable which takes the values, $\frac{-1}{2}, 0, \frac{1}{2}$ with probabilities $\frac{1}{4}, \frac{1}{2},\frac{1}{4}$ respectively. Determine the average probability of error as a function of $\sigma_n^2$.

My try: I think the probability of error is $$P_e = \mathbb{P}(y_m<0 |a_m = 1) + \mathbb{P}(y_m\gt 0|a_m = -1)$$According to the definition of conditional probability, we have $$\mathbb{P}(y_m<0 |a_m = 1) = \frac{\mathbb{P}(y_m<0 \cap a_m = 1)}{\mathbb{P}(a_m = 1)} = \frac{\mathbb{P}(n_m+i_m<-1)}{\mathbb{P}(a_m = 1)}$$ By law of total probability $$\mathbb{P}(n_m+i_m<-1) = \mathbb{P}(n_m+i_m<-1\cap i_m = \frac{-1}{2}) + \mathbb{P}(n_m+i_m<-1\cap i_m = 0) + \mathbb{P}(n_m+i_m<-1\cap i_m = \frac{1}{2}) = \mathbb{P}(n_m<\frac{-1}{2}) + \mathbb{P}(n_m<-1) + \mathbb{P}(n_m<\frac{-3}{2})$$
Each term can be written as a function of $\sigma_n^2$ easily. So it seems we need $\mathbb{P}(a_m = \pm1)$ instead of $\mathbb{P}(i_m = \pm\frac{1}{2},0)$ which can't be true because we should certainly use values of $\mathbb{P}(i_m = \pm\frac{1}{2},0)$. What's my mistake here? Also is there any difference between "average probability of error" and "probability of error"?

Note that
$$
\mathbb{P}(n_m+i_m<-1\cap i_m = \frac{-1}{2})
$$
is equal to
$$
\mathbb{P}(n_m < -0.5)\mathbb{P}(i_m = -0.5) = 0.25\mathbb{P}(n_m < -0.5)
$$
since $n_m$ and $i_m$ are independent.
Another way to obtain the solution is to add the the probability of these three disjoint events:

$n_m < -1.5$,
$i_m = 0$ and $-1.5 < n_m < -1$,
$i_m = -0.5$ and $-1 < n_m < -0.5$

In this case, we can add the three probabilities because the events are disjoint. To keep them disjoint, it is important to keep the noise ranges disjoint too; otherwise, you would be adding them multiple times.


Phase Correlation and Negative Shifts I am implementing phase correlation algorithm to determine shift between two images.

It generally works, but I am not sure how to interpret the resulting shift.

Pseudocode:

ft1 = FT(input1)
ft2 = FT(input2)

for each pixel p
  product = ft1[p] * conj(ft2[p])
  cross_power[p] = product / abs(product)

cross_correlation = IFT(cross_power)

where FT and IFT are Discrete Fourier Transform and its inverse, respectively.

Here are input images and resulting phase correlation taken from Wikipedia:

The peak is located at (20, 21) which corresponds to shifting of the second image.

I have obtained similar result, but it comes out flipped:

When I flip the input images the peak appears at the top left corner of the phase correlation image but how to interpret this as a negative shift (-20, -21)? Based on image quadrants?

I am using Fourier transform routine from Math.NET Numerics package. It provides only 1D FT, so I am doing transforms of rows and then columns to obtain 2D FT.

How to extract shift from the peak position allowing for negative shifts as well?

Why is the peak located differently in my result?

the matter is due to this one:

product = ft1[p] * conj(ft2[p])

You can change to:

product = ft2[p] * conj(ft1[p])

I guess input1 is the reference image.


How to Actually Use FFT to Implement Low Pass Filter have searched high and low for a practical guide to implement a simple 50Hz low pass filter on a signal.I understand the theory but its the practical but that's missing. If I was using MATLAB it would be trivial because there are so many examples but I am coding this by code.

This is where I am at:

1: I have a series of 8192 bytes sampled at 1000Hz in an array (0-256)

2: I have generated the frequency domain representation of a simple step filter in a symmetric array (+/-F, Ire, Iim) :

0,1,0
0.1220703125,1,0
0.244140625,1,0
0.3662109375,1,0
0.48828125,1,0
0.6103515625,1,0
.....
-0.6103515625,1,0
-0.48828125,1,0
-0.3662109375,1,0
-0.244140625,1,0
-0.1220703125,1,0

3: I have used Inverse FFT of the filter to obtain the time representations of the filter as (t, Ire, Iim) [as expected the im parts are ~0 and the filter gives the expected decaying ringing]:

0.000,819,0
0.001,805.600848087199,4.11892742135933E-12
0.002,766.191083986647,7.80009390410896E-12
0.003,703.07869137741,1.07704956064936E-11
0.004,619.930866652693,1.26768595620774E-11
0.005,521.519198753447,1.32877042702262E-11
0.006,413.390428478234,1.26576527037514E-11
0.007,301.48763797423,1.08400510789863E-11
8,0.008,191.750273447113,7.79765141345479E-12
....
8.186,413.390428478235,7.79385628257856E-09
8.187,521.519198753448,9.84392145575441E-09
8.188,619.930866652693,1.17116023545805E-08
8.189,703.07869137741,1.32913150485692E-08
8.19,766.191083986645,1.44921696865197E-08
8.191,805.600848087196,1.52439926792702E-08

OK so what do I do with this data???

From what I have read in the literature the thoery is that in the freq domain you multiply the spectra but because this is the DFT the usual course is to convolve in the time domain but the question remains:

1: Do I simply multiply the two time series together (doesn't feel right)?

2: If not how is the weighting applied?

3: How do I normalise any output back to 0-256 values?

I realise this is a pretty basic question but there doesn't seem to be a really basic answer or example out there :-(

Thanks

iced

this is not really a question but is an answer to a question I asked here and the answer is a little long winded so I thought I would post the answer here as a guide to others. I have done this because I found it demonically difficult to find a real "how to" for do basic filtering. I found a few guides to this helpful here and here So here goes, I hope it helps some others:

I had a data set of biometric data that was full of hum and noise, most of which was >30Hz and I had tried various "averaging" techniques to get rid of this with varying success so I turned to FFT to do the job. While I had studied transforms in general at uni and understood the basic theory and maths I needed to code this for a small app for a tri club (no its not commercial!).

So here is the psuedo code with comments:

SF: Sample Frequency
Samples: is array with samples 0..255
N: is number of samples (8192)
NP: is root of 2^NP = N (13)
CF: cutoff freq of filter (30)
W: number of coefficients in filter (40)
dare: real data input
daim: imaginary data input

for i = 0 To N – 1
{
// time axis for graphing results
time(i) = i / SF

// freq axis for graphing, note symetric FFT so need -ve freq
if i <= N / 2 then // Positive frequencies
freq(i) = i * SF / N
else // Negative frequencies
freq(i) = i * SF / N - SF

// copy input data to local array
// and normalise to 127 = 0, centering the input/transform

dare(i) = Samples(i) – 127
daim(i) = 0
}

// fft the input data
FFT(NP, N, dare, daim)

// Create the impulse response of the filter in the frequency domain
// this could be any filter but I use a simple box here: 1 in the passband and 0 elsewhere. Could use any filter here:

// Note that I used symmetric complex FFT: it is necessary
// fill in the negative frequencies as well.

for i = 0 To N - 1
{
re(i) = 0
if i <= (CF * N / SF) or i >= (N - CF * N / SF)
re(i) = 1
im(i) = 0
}

// Inverse filter into the time domain
// Note: imaginary parts should be ~ 0 after iFFT
IFFT(NP, N, re, im)

// Hann window to truncate the filter series at say +/-20
// can use other window formulas (e.g. bartlett) if desired
for i = 0 To N – 1
{
window(i) = 0
if i <= Won2 then
window(i) = (0.5 + 0.5 * Cos(i * 2pi / (W + 1)))
else if i >= N - Won2 then
window(i) = (0.5 + 0.5 * Cos((N - i) * 2pi / (W + 1)))
}

// Apply the Hann window to the filter
for i = 0 To N – 1
{
//careful here if complex multiply required
re(i) = re(i) * window(i)
im(i) = 0
}

// can IFFT and check for freq response and back again

// Shift the filter to make it i.e. +ve times only
//shift lower coefs up
for i = Won2 To 0 Step -1
{
re(i + Won2) = re(i) //
}
//wrap -ve coeffs around
for i = 0 To Won2 – 1
{
re(i) = re(N - Won2 + i) //wrap N -> 0, N - 1 -> 1 etc
re(N - Won2 + i) = 0
}

// Transform filter back into the frequency domain
FFT(NP, N, re, im)

//multiply input data by filter COMPLEX MULT REQUIRED!!
// may be able to shorten this step if desired
for i = 0 To N – 1
{
tmpre = dare(i) * re(i) - daim(i) * im(i)
tmpim = dare(i) * im(i) + daim(i) * re(i)
dare(i) = tmpre
daim(i) = tmpim
}

//inverse back to time domain
IFFT(NP, N, dare, daim)

// copy back to input normalised and re centered
// may need to normalise/scale at this stage
for i = 0 To N – 1
{
Samples(i) = dare(i) + 127
}


Transmission capacity in "duobinary coding" In "duo-binary signalling" how exactly is the transmission capacity increased?

Because , when we encode the binary bits, the bit duration remains the same, but the amplitude will change. Then how can the transmission capacity be better?

If you want to transmit a symbol sequence $A_k$ using baseband pulse amplitude modulation (PAM), the transmitted signal is

$$s(t)=\sum_kA_kp(t-kT)$$

where $p(t)$ is the transmit pulse, and $T$ is the symbol interval. If you want to avoid intersymbol interference (ISI), the pulse function $p(t)$ must satisfy the Nyquist criterion, which says that its value must be zero at multiples of the symbol interval:

$$p(0)=1,\quad p(kT)=0,\quad k\neq 0\tag{1}$$

The minimum bandwidth required to satisfy (1) is $1/2T$, but the corresponding pulse $p(t)$ is an ideal low-pass filter, i.e. a sinc function, which is of course impractical. So in practice, Nyquist pulses with a higher bandwidth (excess bandwidth) are used. This means that for a given bandwidth, the possible symbol rate is reduced as compared to the theoretically maximum possible value.

With duobinary signalling practical systems can be built with zero excess bandwidth, because the data symbols are filtered such that a spectral zero is introduced at half the baud rate. Note that this is a very simple filter unlike the ideal low-pass filter required for the minimum bandwidth Nyquist pulse. The price paid is the deliberate introduction of ISI and the introduction of a ternary (instead of binary) signal, which requires a higher SNR to achieve a given error probability.


Calculating SFDR from FFT of ADC output A few months ago, as part of my summer internship, I had performed interfacing of an ADC with FPGA (Zynq-7000). The goal was to take an analog signal from a signal generator, feed it to the ADC and present the digitized output to the FPGA. The FPGA was supposed to perform FFT on the data through an algorithm which we had written on Vivado Suite. The magnitude spectrum was successfully obtained on the waveform viewer.

We then decided to calculate the SNR and SFDR values from the data. The way I proceeded for computing SFDR was to record the value of the amplitude of fundamental frequency (which kept fluctuating, but not much) and value of the next highest peak (the highest spur) in the first Nyquist zone. The ratio of these two values would fetch me the SFDR.

The input to the ADC was a 200 kHz sinewave.
To check for consistency, we performed three different sized FFTs, 2k-, 8k- and 32k- and in each trial the strength of the analog signal was varied from -60 dBm to 0 dBm (in steps of 10dBm).  The value of SFDR always hovered around 23 dBc.

The value of SFDR given in the ADC datasheet is around 80 to 90 dBc, but we ended up with a value much less than that (23dBc).

I understand that the conditions under which the datasheet values were calculated were quite different from ours in terms of strength and frequency of input signal.

My major concern which I wish to clarify from the experts sitting on this platform is whether or not my procedure was proper enough, or atleast worth appreciating? Was there any major flaw(s) which I may have failed to notice in my approach? I am sure there certainly were. This concern arises because whenever I think of mentioning this project in my resume, I wonder if the examiner will ridicule me of missing anything fundamental. There wasn't a subject expert at that time who could point out if I was going on the right path.

I would tend to believe that there is an indeed a flaw in the approach as the datasheet SFDR is achievable given the test conditions which are usually specified in the datasheet as well.

23 dB is SIGNIFICANT and I suspect your sine wave is not really a sine wave and you are seeing a harmonic (or an aliased harmonic).

Here are some possible culprits, some of which you may have ruled out already:

You did not specify what frequency the next highest peak was, but harmonics are usually excluded from SFDR.  (400KHz, 600KHz, 800KHz, etc)

Did you sufficiently filter the analog signal prior to sampling? If not any higher aliases will fold into band. (Changing the sampling rate and retesting, not the FFT size, would give clues if this is occuring with such a strong spur).

Did you observe your test signal in the analog domain to see if it is cleaner than the SFDR you are trying to measure? If you did not look at this (with a spectrum analyzer for example) is possible your test signal itself has such a spur.

Similarly, did you confirm that your clock signal is cleaner than the SFDR you are trying to measure? Spurs on the clock will translate to spurs in the output signal as well.

Did you set the input signal level appropriately and is your ADC interface per the vendor's recommendations? Typically the best SFDR is when the input signal is one to two dB below full scale where clipping starts to occur, otherwise non-linearities the front-end of the ADC with start to reduce SFDR. If your input signal was too low then you will clearly not get the SFDR of the datasheet (typically the spurs stay at the same level, so you maximize SFDR with an input level as stated above).

As far as your resume, I don't think anyone will ridicule you given what you accomplished in this project, as long as you present what you did, what you think you understand from your efforts and admit openly that you still have a lot to learn and show that you are actively learning and interested (which you are by coming here!). That is probably one of the best skills someone new to the field coming out of an internship can have. I think it was a great project you got to do and would be impressed hearing more about it. At your level, if I was a hiring manager interviewing you, I would look for your interest and enthusiasm and complete honesty in what you understand and willingness / ability to learn. If you thought or acted like you knew something that you didn't really understand, that would typically come across as a big negative and a reason to move on to the next candidate for a technical position. Best to "under-promise and over-deliver" by immediately saying you are not an expert, you still have a lot to learn etc and then while you have the floor dive into what you did and what you learned and how excited you were about doing it. (My 2 cents being on both sides of that fence).


How to scale a signal to get desired variance I am new to dsp, and I'd really appricate it if someone please help me with this problem. I guess it's a basic simple question, but I cannot get it write. Here is the question:

The variance of a random signal generated using MATLAB command randn is approximately 1. We need a random signal with variance 0.01. With which factor you need to scale the signal?

With amplitude? Like we need to multiple the amplitude of the signal by 0.1?

standardizing a Normal random variable $x$
$$
\frac{x -m}{\sigma} \sim \mathcal{N}(0,1)
$$
so if $y$ is a random variable that is Normal with mean zero, with standard deviation $1$,
$$
\sqrt{\sigma^2}y \sim \mathcal{N}(0, \sigma^2)
$$
This standardization actually works for more than just Normal Distributions.


Checking if u[n] and v[n] are orthogonal Let $x[n]$ and $y[n]$ be two orthogonal baseband signals, with bandwidth B. Check if the following signals are orthogonal:

$u[n] = x[n]*\cos(2*\pi*f_c*n)$

$v[n] = y[n]*\cos(2*\pi*f_c*n)$

I'm not sure what the problem asks to be honest.
1) What exactly does he mean by bandwidth?? I only know the definition of a bandwidth of a sinc function, but how can I think of it generally?

2) I assumed the case for a rectangular function and reduced the problem to:

$\sum_{n=-\frac{1}{B}}^{n=\frac{1}{B}}\frac{1}{2}\cdot x[n]\cdot y[n]\cdot \cos(4\cdot \pi\cdot f_c\cdot n)$

But I don't know where to go from there. It's also weird for me to evaluate since $\frac{1}{B}$ is likely not an integer.

EDIT: also, $f_c$ >> B

It's essentially asking you to show that two signals that are orthogonal in the baseband (their original form) remain orthogonal after amplitude modulation process. In other words, amplitude modulation preserves orthogonality.

Think about the question this way. You have two signals that are orthogonal, meaning that

$$ \sum_{m=-\infty}^{\infty}x[m]y[m] = 0 $$

This is a bit ugly, since we have to deal with infinities. However, your problem also says that they both have bandwidth $B$, which means that in frequency domain they will have finite width.

Let's define $\mathscr{B}$ as

$$ \mathscr{B} = \frac{B}{F_s}$$

which means that in your discrete Fourier transform signals $X[k]$ and $Y[k]$, the only non-zero samples exist for $-\frac{\mathscr{B}}{2} \le k \le \frac{\mathscr{B}}{2}$. Since DFT is an orthogonal transform, signals orthogonal in time domain will remain orthogonal in frequency domain, i.e.

$$ \sum_{l=-\infty}^{\infty}X[l]Y[l] =\sum_{l=-\mathscr{B}/2}^{\mathscr{B}/2}X[l]Y[l] = \sum_{m=-\infty}^{\infty}x[m]y[m] = 0 $$

We know that multiplying a signal by a cosine shifts is away from DC in frequency domain, and the fact that the modulating frequency $f_c$ is larger than $B$ means that this won't cause any problems with overlapping modulated bands (think about what would happen if this weren't the case).

I think this gives you enough hits to take it from here.


How do scaling of the coefficients for the poles of a biquad filter affect the gain? Given a digital filter consisting of a chain of digital biquad filters.
Given a digital biquad filter with a normalized Direct form 1 difference equation: $ y_n = b_0x_n + b_1x_{n-1}+b_2x_{n-2} - a_1y_{n-1} - a_2y_{n-2}$

My goal is to implement the biquad filter with few multiplications. For this end I play with the coefficients. A first question is:

Can I change coefficients $a_1$ and $a_2$ proportionately, and still retain the frequency response of the biquad filter?"

The rest of this post assumes an affirmative answer.
Many of my filters will have $b_0 = b_2$. As a start I change $b_0, b_1, b_2$ so that $b_0$ (and $b_2$) are powers of 2. This changes the gain of the biquad filter, this will be compensated later in the signal chain.
I want to change the coefficients $a_1$ and $a_2$ in a similar way, to eliminate yet another multiplier. I reason that this should be possible - the feedback is dependent on the proportion of $a_1$ and $a_2$. I am, however, unable to fathom how this affects the gain of the biquad filter. So my second question is:

How does proportional changes to $a_1$ and $a_2$ affect the gain of the biquad filter?

A concrete example

Assume the following coefficients:

a0 =  1.000000
a1 =  0.305805
a2 = -0.192467
b0 =  0.236979
b1 =  0.412704
b2 =  0.236979

If used "as is", four multipliers are needed (we would add $x_n$ and $x_{n-2}$ before multiplying with $b_2$).
To remove one of the four multiplications is easy. I modify the non-recursive coefficients b0, b1 and b2 in such a way that b0 is a power of 2.

a0 =  1.000000
a1 =  0.305805
a2 = -0.192467
b0 =  0.236979/0.236979*0.25 = 0.25
b1 =  0.412704/0.236979*0.25 = 0.435380
b2 =  0.236979/0.236979*0.25 = 0.25

I now have a gain of $b_2/0.25$ I now modify $a_1$ and $a_2$:

a0 =  1.000000
a1 =  0.305805/0.305805*0.25 = 0.25
a2 = -0.192467/0.305805*0.25 = -0.157344
b0 =  0.250000
b1 =  0.435380
b2 =  0.250000

Experimentally, it seems the gain now changes with a factor 1.25, but how to derive this?

Sorry, that doesn't work. You have to scale all numerator (the "b") or all denominator coefficients (the "a") at the same time. That doesn't work for the denominator since you all practical implementations require $a_0=1$. Scaling $a_1$ and $a_2$ without scaling $a_1$ will move the poles and change the overall shape of the transfer function and not just the gain.


Tracking using color I am working on tracking an object in a video by using its color. what i did is that i have converted the color of a selected object into LAB color space since it is a uniform color space. and then in each frame of video i am searching for same color by finding color difference in LAB color space. Now main problem is arising when more than one object of same color is appearing in the video scene. one method i adopted for avoiding detection of wrong object is to search for same color in limited area  around last detected position of the object. but this is causing failure to detect fast moving objects and leaving my methods suitable for slow moving targets only. so i want to know is there any method i can apply for avoiding wrong detection.

You can check the area of the same-color region. It should be close to the last one and it should be around some value (you can filter out very small and very big regions).
You can check rectangularity, circularity, convexity, angle, width, height depending on your object to be tracked.
If your false-detected regions are very small, you can try to a bit smooth the image first and then check for values.
Besides that if the false detected objects are very similar to the correct ones, your method seems to be the most feasible way to do it. You can consider making it faster by choosing different platform or optimizations.


what is the reason for the Change in DFT analysis after downsampling firstly the original signal's dft analysis is this the signal is sammpled at 1000Hz

the signal is then downsampled under different rates such as 25Hz

this is the dft analysis of the signal downsampled under 25Hz

what is the reason for this phenomenon?

Because you have resampled the signal to a sample rate below twice the highest frequency in the signal, you are seeing the signal plus an aliased image of the signal both summed together.  Two partially overlapped rectangles will sum to produce the "head-and-shoulders" shape you see in the aliased spectrum.


Differences between FFT from acceleration, velocity and displacement data of a dynamic system I am new to vibration mechanics and signal processing. Recently, I’ve performed an experimental test with an accelerometer mounted on a plate (shaker). The plate vibrates with exact frequencies 10 and 20 Hz. The purpose of the experiment is to extract the acceleration vs. time and FFT diagrams. I have a basic and fundamental question.

What are the differences between applying FFT to acceleration, velocity or displacement diagram of a dynamic system? Does FFT give the same peaks for the three mentioned diagrams so that I assume that the peaks are natural frequencies of the vibrating system (10 and 20 Hz in this problem)?

I would be appreciated if you could help me with these questions.

I have the same problem and I made some simulation and calculations to understand this phenomenon.

Facts:

velocity is an integral of acceleration
displacement is an integral of velocity
sine wave having frequency f is sin( 2*pi * f * t ) where t is time in seconds
integral of sin( n*t ) is -1/n * cos( n*t )
double integral of sin( n*t ) is -1/n^2 * sin( n*t )
To make equations cleaner, I'll use k = 2*pi

Calculations:

If you have hypothetical acceleration measurement signal having frequencies 1Hz, 10Hz, 100Hz with equal amplitudes, then:

acceleration: sin(k*t) + sin(k*10*t) + sin(k*100*t)
velocity: -1/(k) * cos(k*t) - 1/(k*10) * cos(k*10*t) - 1/(k*100) * cos(k*100*t)
displacement: -1/(k^2) * sin(k*t) - 1/(k^2*10*10) * sin(k*10*t) - 1/(k^2*100*100) * sin(k*100*t)

So, it seems that:

Both acceleration, velocity and displacement signals have the peaks at the same place: 1Hz, 10Hz, 100Hz
the 10Hz and 100Hz peaks of velocity are lower and the peaks of displacement are MUCH lower that the peaks of acceleration
the 100Hz peak of displacement is 10,000 times lower that 100Hz peak of acceleration

Simulations:

I made short R script which created the acceleration signal above. Then I integrated it (cumulative sum) to get velocity and displacement signals and finally took fft of them. The results verified the calculations


How to segment frequency bands after FFT and calculate power? I have been working with EEG signal. I have raw EEG data for a 2-second (i.e. $t=2$) duration. My objectives are:

Transforming these raw data from time domain to frequency domain.
Segmentation into two frequency bands: $\left[4 - 8\textrm{ Hz}\right]$ and $\left[9 - 13\textrm{ Hz}\right]$.
Calculating the average power of these two bands separately.

I am doing the signal analysis in MATLAB. Sampling frequency is $512\textrm{ Hz}$. I have set the frequency resolution to $0.5\textrm{ Hz}$. For the first objective, I have used the fft function to transform the data.

freq = fft(x);  // x is my eeg data

I'm not sure on how to proceed in the second objective but I did in the following way. Since the frequency resolution is $0.5\textrm{ Hz}$, the bin numbers ($9$ to $17$) represent $\left[4 - 8\textrm{ Hz}\right]$. So, I did

band1 = freq(9:17);

Finally, to calculate the average power of a band, I used bandpower function.

power = bandpower(band1);

Are the proceedings correct? If not, please help me. I am new to signal processing.

Your procedure is mostly correct. You should also add the bins corresponding to the negative frequencies -8Hz to -4Hz, which appear in the final bins of the fft result: the first bin corresponds to frequency 0Hz, the second to frequency $f_s/N$ = 0.5Hz, and so on. The last bin corresponds to frequency $f_s(N-1)/N$, which also corresponds to frequency $-f_s/N$.

Also note that bandpower calculates the average power of a signal, not of its Fourier transform. It should be equivalent (following Parseval theorem), except for your missing negative frequencies.

I'd suggest using bandpower directly in the form p = bandpower(x,fs,freqrange). Here you provide the signal, the sampling frequency, and the frequency range, which is what you want.


Time Frequency Signal Analysis in Mathematica (MMA) I want to know " Whether there is any Tool Box in Mathematica (MMA) for the Time-Frequency (TF) Signal Analysis".

I am well-versed in MMA Programming, so want to do TF Signal Analysis in MMA. I think if there is toolbox of TF Analysis, then it will be of very much great help in long programming.

Thanks in Advance!!!!

If you are using an older version of Mathematica (pre v.8) and are interested in wavelets - yes, you need an add-on to perform wavelet analysis. More about it here. If you are using v.8 or above then everything wavelet-related is built-in. If you are interested in Fourier analysis - it is built in (since v.1). To showcase the capabilities of the time-frequency analysis functionality:

Suppose we have the signal $2 \exp \left(-\frac{t^2}{10}\right) \cos (5 t)$

We can visualise its Gabor transform

Wigner transform (notice the artifacts)

and the (somewhat) corrected Gabor-Wigner transform

A second test signal

and its spectrogram calculated using partitions of length 256, offset 1 and BlackmanHarrisWindow

its periodogram (with the same options as above)

Now we move onto wavelets
Scalogram after a continuous wavelet transform with a Morlet wavelet

and the respective scalogram after performing a discrete wavelet transform

Just because I find it pretty - a scalogram of a simulated noise in 3D


FFT of image data: "mirroring" to avoid boundary effects I load and display an image of some rice in Matlab:

g = imread('rice.png');
imshow(g);

I take the FFT of this image and shift it:

G = fft2(g);
imshow(log(abs(fftshift(G)) + 1), []);

If I place a x and y axis trough the center of the image; I find that the image is symmetric g(-x,-y)=g(x,y). For a 1D signal we have that the FFT of a real signal has a symmetrical real part and an asymmetrical imaginary part. I guess that is what we see here in 2 dimensions?

Since the original image is darker at the bottom than at the top, there is a strong horizontal discontinuity at the periodic boundary causing the vertical line in the FFT.

I want to get rid of this boundary effect. A common approach to this seems to be windowing.

However I want to solve this problem by a technique I found in a paper called "mirroring".
The paper was not very specific so I need your help in figuring out this approach :-).

First I create a symmetric "tile" from the original image:

tile=[flipdim(g,2) g; flipdim(flipdim(g,1),2) flipdim(g,1)];
imshow(tile);

Now I take the FFT of this "tile":

Tile=fft2(tile);
imshow(log(abs(fftshift(Tile)) + 1), [])

The vertical line seems to be (almost) gone: good.
However the mirroring seem to have introduced more symmetry.

What is the correct result: the FFT of the original image or the FFT of the "mirrored" image?

Is there a way I can "mirror" so that I both get rid of boundary effects and get a purely real FFT?

Thanks in advance for any answers!

The FFT of the original image is correct. The artifacts you are seeing are typical for the DFT, because the DFT basis functions have difficulty representing non-periodic signals. Even though the DFT is of finite length, it actually represents a periodic signal that extends to negative and positive infinity. The DFT's basis functions are all sinusoids with periods that are integer divisors of the DFT size, so they all start and end at the same value, and it is hard for it to fit signals that are not periodic in this way. So, imagine tiling your image: there is indeed an abrupt discontinuity at the edges. The edges do not match up well, and that is the reason for the boundary effects.

In audio signal processing, where the FFT is used often, a very common approach is to use a window function, as you noted. This tapers the edges off towards 0, and reduces this effect.

In image processing, the DCT is usually used rather than the DFT, because it imposes different symmetry constraints that give better results. In other words, it's used to prevent exactly the issue you're seeing. In contrast to the DFT's basis functions (remember how they all start and end at the same value), the DCT's basis functions are integer divisors or half of integer divisors, so many of them start and end at different values. As a result, the implied boundary conditions are different: the signal is assumed to be symmetric about its edges, which is actually somewhat similar to the "mirroring" you were experimenting with. Here is another post by me with a little more information on the DCT: https://dsp.stackexchange.com/a/362/392

So, as a short summary, if you are absolutely set on using the FFT, you might want to try windowing. However, the best choice is probably to use the DCT, which implies boundary conditions similar to your "mirroring" idea, and as a result, handles images better.


Finding the convolution $a^n u[n]$ with $b^n u[n]$ Hi I'm trying to solve the problem when studying for an upcoming test. The given solution is $$y[n]= \frac{1}{b-a}(b^{n+1}-a^{n+1}) \quad \text{for } n\ge0 \ .$$ However, I'm not sure how to reach this point. We have to also verify using the z-transform. Any help would be appreciated so much.

Just do it the hard way by writing out the individual formulas for the general convolution sum
\begin{align}
y[n] &= \sum_{k=-\infty}^\infty a[k]b[n-k]\\
&= \sum_{k=0}^\infty a[k]b[n-k] &\scriptstyle{\text{because }a[k]=0~\text{whenever } k < 0,}\\
&= \sum_{k=0}^n a[k]b[n-k] &\scriptstyle{\text{because }b[n-k]=0~\text{whenever } k > n,}
\end{align}
which gives
\begin{array}{rclcl}
y[0] &= &a[0]b[0] &= &1\\
&&&= &\displaystyle\frac{b^1-a^1}{b-a}\\
y[1] &= &a[0]b[1] + a[1]b[0] &= &b+a\\
&&&= &\displaystyle\frac{b^2-a^2}{b-a}\\
y[2] &= &a[0]b[2] + a[1]b[1] + a[2]b[0] &= &b^2 + ba +a^2\
&&&= &\displaystyle\frac{b^3-a^3}{b-a}\\
&&&= &\displaystyle\frac{b^3-a^3}{b-a}\\
\vdots &= &\ddots~~~\scriptstyle{\text{Hey, Ma! I think I see a pattern developing here}}\\
\\
y[n] &= &a[0]b[n] + a[1]b[n-1] + \cdots + a[n]b[0]
&= &b^n + b^{n-1}a + \cdots + ba^{n-1} + a^n\\
&&&= &\displaystyle\frac{b^{n+1}-a^{n+1}}{b-a}
\end{array}
if you can recall the formula for the sum of a geometric series.


SNR of a 16-bit DSP with 12-bit ADC, 40-bit accumulator According to the 6dB-per-bit rule, a 16-bit DSP would provide 96 dB of dynamic range, approximately.

If the ADC and DAC are both 12-bit long, the DSP is still considered to have 96 dB, or the dynamic range is now 72 dB?

Also, I would like to know if the guard bits provided by the accumulator (40-bit) increase in some way the dynamic range.

Thank you very much.

Let me add the following additional important consideration to dynamic range on the ADC side of the equation to the other answer provided.

The SNR as characterized with a full scale sine wave for a typical (uniform quantization noise) ADC is well established to be given by the following relationship:

$$SNR = 6.02 dB/bit + 1.76 dB + 10 Log(N)$$

(For more on how that equation is derived, see What are advantages of having higher sampling rate of a signal?)

Which demonstrates the consideration given for oversampling to dynamic range.  (And just to mention, there are other ADC typologies such as Delta-Sigma with noise shaping offering even more dynamic range improvement due to oversampling; however in any case the sampling and the distribution of the quantization noise is limited to the architectural implementation of the ADC itself).

The example below illustrates this with a 10 MHz BW signal at a digital IF of 25 MHz, sampled at 100 Msps.  This example was done with Analog Devices HMCAD1520 which is a 14 bit ADC with a SINAD (Signal to Noise and Distortion) of 72.5 dB. 72.5 dB translates to an ENOB (Effective Number of Bits) of 11.8 (using the same 6dB/bit+1.76 dB equation). This means that with a full scale sine-wave, the additive noise due to the ADC sampling process will be 72.5 dB lower than the power level of the sine-wave itself. However this noise, well approximated as a white noise source is spread evenly across the sampling spectrum. Thus as shown in the figure below, the portion of the noise that is within the bandwidth of the signal is 1/5 or -7 dB lower.

The rest of the noise, out of the band of interest, can be filtered out in the DSP, down to the limits of the datapath (unless special provisions are made to effectively increase the datapath precision as suggested in the other answer). For this filtering purpose, extended precision accumulators are particularly important, as described in more detail in other posts.  (for example, see what is the suitable design Method to the filter?)

Some of the distortion is in spurious frequencies, but referring to the datasheet for this device, the SFDR (Spurious Free Dynamic Range) is 83 dB, so in this specific case, with this oversampling rate, the dynamic range will not be limited by spurs.

So given a sine wave at full scale (to establish a consistent metric point), we can establish the ADC noise floor.  However for dynamic range, we also need to consider the headroom required for a modulated signal. For this reason I have created the following graph which is applicable to Gaussian distributed waveforms in magnitude (such as OFDM and CDMA both of which are well approximated by Gaussian distributions), and shows the optimum set point in rms level below full scale in order to maximize dynamic range for the case of a modulated waveform based on the number of bits. This is applicable throughout the digital system everywhere in the datapath, not just at the ADC, so is a very important consideration. This shows specifically the effect of clipping noise as the rms signal level approaches full scale, and the effect of quantization noise increasing relative to the signal, as the rms signal level is reduced from full scale.  The term "dBFS DC" is the level in dB below full scale of a DC signal (a sine wave at "0 dBFS" is -3dB below '0 dBFS DC"). The 0 dBFS DC mark is also applicable to full scale for a single tone in a complex datapath (with I and Q).


Power spectrum of a pure sine problems I am writing a C++ program to use to analyse oscillations and create a power spectrum.
I am using FFTW3.

To test the system, I am using a perfect sine function.

But no matter what parameters I try, I run into the same problem. I find the correct frequency, but also many 'overtones'.

First I assumed the problem was my way of chopping up the signal or how I calculate the running mean. But the problem is there when I just have all data in a single big window, bypassing my suspect code, still creating the pattern in the pictures.

I then assumed it was an aliasing problem, but even 1024 datapoints in a period, and they don't go away. They just lose shape and appear randomly. I assumed it was too little periods in my window. But more periods per window and they just become less frequent and stronger.
Then I assumed it may be spectral leakage, so I assumed I may need to use a windowing function like a Hamming, but that creates an odd slant. Is that the trade-off?

Or do I just set everything that is small enough to zero? I was just expecting a single strong point peak and the remainder to be so near zero, the log is -6 or lower everywhere else uniformly.

If I can't figure out how to get the right parameters for the best signal-to-noise ratio, I am not sure how I will achieve that when I start to work with the actual data.

I assumed it would be simple to pick a window that contains 32 periods with 32 discrete points in each period, use Welch's method and be done with it.

I am not an engineer who took many many classes on this kind of stuff. So I just don't have the proper training. I am working on a project of 4 weeks. There is so much info out there, it is confusing, but I don't know where to start as I don't have the time to learn stuff properly.

It can't be caused by the limited numbers of decimals of the input sine values, right?

Since FFT treats the signal as if it is periodic you need either to apply a window function (for example hanning) on your signal or make it coherent. In the image3 you attached you may make it coherent by only using data for a number of cycles; use samples 0..31 or 0..63. If the signal is non-coherent it will be seen by the FFT as a concatenated signal with discontinuites that will cause overtones.

EDIT:
A coherent signal may be concatenated with itself without generating any discontinuites.

The blue line is the input data to FFT

And a non-coherent:


Why autocorrelation can be more efficiently calculated using the fft Can anyone explain why autocorrelation can be more efficiently calculated using the fft ?

The cross-correlation between two functions $f(t)$ and $g(t)$ can be seen as a convolution of $f(t)$ and $g(-t)$. The auto-correlation is of course a special case where $f=g$. This operation in the frequency domain corresponds to a multiplication. Therefore, an efficient way to compute the cross-correlation is to compute the Fourier Transform of both $g(t)$ and $f(t)$, multiply them in the frequency domain and the anti-transform the result, to go back to the time domain.

This approach apparently looks more computationally expensive, however in the discrete case, if we use the fast algorithm of the FFT (Fast Fourier Transform) to compute the DFT (Discrete Fourier Transform), which uses $N\log N$ arithmetic operations, it leads to more efficient algorithm than the actual computation of the cross-correlation function in the time domain, which requires $N^2$ operations, where $N$ is the number of samples. So the number of operations is much lower if $N$ is particularly high, which is often the case ...

You can watch this video if you want to know more about this https://www.youtube.com/watch?v=iTMn0Kt18tg


Finding the Amplitude of a Sinusoid in Noise I'm trying to solve for the amplitude and frequency of a sinusoid embedded in zero-mean gaussian white noise.

I am supplied a sample file of a 40000 element array. I first took the autocorrelation function in matlab and then took the fft of this to find the power spectral density. From there I was able to find the frequency of the sinusoid due to the peaks. I am unsure how I could go about finding the amplitude, as I am unsure how I can go about normalizing the PSD or Autocorrelation function.

ps this is all done in matlab, thanks!

If you know you have a single sinusoid in Gaussian noise, then a least squares type parametric estimator might be one method of estimating the amplitude. See  DSP Related, specifically:


Particle filter, resampling step implementation I am not very good in statistic neither in probability theory, and I do not completely understand the mathematical background of Monte Carlo simulation or SIR sampling, but I'd like computer vision and I'm experimenting with it hoping to get some nice results with already known techniques (I don't want to reinvent the wheel neither to find some new method).

I have done an approximate implementation of particle filter with Bhattacharyya distance as liklihood score function. What I mean with approximate refers to the resampling part.

Now, my rough algorithm does something like that, in python pseudo-code:

for each particles:
    p.position = move with a random walk motion model

    p_histogram = calculate histogram that is under the particle
    p.weight = bhattacharyya_distance( p_histogram, model_histogram )

sort(particles, less_weight_comparator)
# now particles[0] has less weight --> less distance
# and particles[nParticles-1] has more weight --> more distance

# find best k particles
# i.e., the particles that has weight less then a threshold 0.25
k = -1
for each particles:
    if p.weight < 0.25:
         k++;
    else:
         break

# resample all the other particles near those k best ones
if k!= -1:
    for i in range(k, nParticles):
        # choose a random number between 0 and k
        new_index = rand(0,k)
        # resample the particle near the k_th
        # with a gaussian random with 0 mean and 0.6 variance
        particles[i].position = particles[k].position * random_gauss(0, 0.6)
else:
    print "target is lost"

predicted_position = particles[0].position

this work quite well, but not every time, so I'm trying to improving it following some famous implementation, like the one in the Arulampalam_etal_2002 paper.

But.. I'm confused and I don't understand what is happening in the resampling stage.

With my rough and superficial algorithm I'm always aware of what is happening: particle with the weight less then my threshold are too far from the histogram model, and if k == -1 the target is lost. It is clear, in my mind.

In the paper mentioned (and in all the other particle filter papers) I'm reading about normalizing weight as all weights sums up to 1, and calculate a degeneracy value $N_eff$ (Equation 51) and.. I'm lost.  My n_eff has values from 249.444 (at the first frame after the particle filter initialization so that some particles has to be near to the target) to 245.452 (more frame later, when the target is lost) and here, again I cannot find a threshold. If I normalize the weights, after the sorting, I have particles weight that goes from values like 0.00418279 to 0.00240669, and I cannot easily find a threashold value that can help me in understanding which are near the target or not, and it depends on how many particles I have so if I add more particles I'm more lost.

More over, in paper like Nummiaro_2002 the Bhattacharyya is calculated inside a gaussian (Equation 10) and I'm not sure on what this means, and I lost an easy way to understand which particle has an histogram near to the model or not.

I'm totally lost in what including the Bhattacharyya into a gaussian means, and why normalizing and how to find the threshold for understand how many particles has to be resamples and which are out which are in.

I hope someone can explain me a little bit, with simple and easy words.

Your question is about the resampling step, and let's focus on that.

Resampling is used in particle filtering to counteract "sample impoverishment" that is the fact that some particles may have very low weights. This is a waste of resources as you want to describe your probability distribution ate best. You should note that it should not interfere (a priori) with the measurement of prediction steps.

The optimal situation is to have all particles having roughly uniform weights and that's the goal set by most resampling algorithms. I propose you use this function which uses something similar to histogram equalization for a ratio resampleof particles (a nd-array were the last line represents the weights):

def resampling(particles, resample):
    """
    Resample weighted particles.

    The particles are grouped in N_pop batches to represent multiple velocities.
    The probability is given for all particles but corresponds to the proba-
    bility of one batch.

    Input
    -----
    N particles as a particles.shape array

    Output
    ------
    N particles as a particles.shape array

    Parameters
    ----------
    resample : gives the ratio of particles that get resampled at every frame (between 0 and 1)

    """
    particles_out = particles.copy()
    N = particles.shape[1]  # number of particles

    # resample a percentage of the particles (the total number keeps to N)
    if resample > 0:
        # draw resample (in %) random addresses of particles that we will reassign
        N_resample = int(N*resample)
        address_resample = np.random.permutation(np.arange(N))[0:N_resample]

        # draw from this subset some addresses uniformly over their pdf using histogram equalization
        proba_resample = particles_out[-1, address_resample]
        proba_resample /= proba_resample.sum()
        address = np.interp(np.linspace(0, 1, N_resample, endpoint=False)+1/2./N_resample,#np.random.rand(N_resample),
                            np.concatenate(([0.], np.cumsum(proba_resample))),
                            np.arange(N_resample+1))
        address = [int(k) for k in address]

        # reassign these particles and set their weight to a uniform value
        particles_out[:, address_resample] = particles_out[:, address_resample[address]]
        subset_weight = proba_resample[address].sum() # should be \approx N_resample / N, that is, resample
        particles_out[4, address_resample] = subset_weight / N_resample #1/N_resample #particles_out[4, address_resample].sum()

        # ultimately, normalize weights
        particles_out[4, :] /= particles_out[4, :].sum()

    return particles_out