How to find inverse Fourier transform of summ of delta functions? I am practicing for my exam that I have this semester and I stumbled upon this one.
How can i find inverse Fourier transform given:
$$
X(j\omega) = \sum_{k=-\infty}^{\infty}\delta(\omega-2k+1)
$$

Using Duality Property, we have $X(j\omega) = \delta(\omega-\omega_{0})$. By using this and rewriting our function using $1-2k = -\omega_{0}$, we get:
$$\begin{align}
x(t) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}\delta(\omega-\omega_{0})e^{j\omega t}d\omega \\
 &= \frac{1}{2\pi}e^{j \omega_0 t}
\end{align}$$
$$
\sum_{k=-\infty}^{\infty}\delta(\omega-\omega_{k})\leftrightarrow \frac{1}{2\pi}\sum_{k=-\infty}^{\infty}e^{j\omega_{k}t}
$$
Hence:
$$
\sum_{k=-\infty}^{\infty}\delta(\omega+1-2k)\leftrightarrow \frac{1}{2\pi}\sum_{k=-\infty}^{\infty}e^{j(2k-1)t}
$$


DTFT of even and odd samples Here to find DTFT of $h(2n)$ they have scaled omega, while in RHS to find DTFT $x(2n+1)$ they didn't, why is that?

Because they made a mistake. The very first equation on the right-hand side is wrong. It should be

$$\textrm{DTFT}\big\{x[2n+1]\big\}=\sum_{n=-\infty}^{\infty}x[2n+1]e^{-jn\omega}=\sum_{n\textrm{ odd}}x[n]e^{-j(n-1)\omega/2}\tag{1}$$

Using the trick they suggested, $(1)$ can be written as

$$\textrm{DTFT}\big\{x[2n+1]\big\}=\frac{e^{j\omega /2}}{2}\left[\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega/2}-\sum_{n=-\infty}^{\infty}(-1)^nx[n]e^{-jn\omega/2}\right]\tag{2}$$


Advantages and disadvantages in term of radar performance of LFM chirp Which are the advantages and disadvantages of a LFM chirp regarding its radar performance?

Assuming this is not homework, you can summarize LFM performance as

Advantanges

LFM is easy to implement in hardware - even more so if you can use stretch (de-chirp) processing. It has been the most popular type of pulse compression.
Good range resolution which is a function of bandwidth.
Autocorrelation (matched filter output) has good side-lobe levels.
The LFM spectrum is "nice" given its sharp band edges. This make it more robust to signal processing chains with lots of filtering.
LFM is Doppler tolerant. This means that the matched-filter output (think SNR) is robust to Doppler shifts from the return signal.

Disadvantages

LFM exhibits range-Doppler coupling. A Doppler shift will shift the peak of the matched filter output to a different delay (range).


FFT and IFFT : filtering frequencies I am sampling at $30 Hz$.

I have a sample of $N=150$.

Here is an example of how my sample looks like, which starts at 54 and ends on 203 on the X Axis

I am interested in observing frequencies between 15 and 20 per minute. ( Cycles per minute ) Hence, I would derive the cycles per second ( Hz ) like this.

$15*0.0166666666667=0.25 \ Hz$

$20*0.0166666666667=0.33 \ Hz$

When I take the FFT, my data looks as follows. ( X axis are the bins to N / 2, not frequency )

To derive the frequency for each bin, I am using this stack answer to construct this logic:

Frequency = ( bin * sampleRate ) / dft.Points

So I might have some output as follows:

Frequency = 0 * 30/300
Frequency = 1 * 30/300
Frequency = 2 * 30/300
Frequency = 3 * 30/300

If I re-plot the same data, where the X axis now represents frequency, this is what I get :

Here is where my question kicks in:
( Note, if something above is wrong in my logic, please point it out, as it could be the reason why I am getting weird results )

I now want to take the IFFT of the dataset, given that the frequency of the bin is less than . Here is some java code to help illustrate how I am doing this with JTransforms:

final double[] ifft = new double[fft.length];
for (int idx = 0; idx < fft.length; idx++) {
    final double fftPoint = fft[idx];
    final Double sampleCount = Double.valueOf(idx * 30);
    final Double binToFrequency = Double.valueOf(sampleCount / fft.length);
    ifft[idx] = binToFrequency <= 0.33333334 ? fftPoint : 0;
}

analyzer.realInverseFull(ifft, false);

Variable fft above is the DFT/FFT array data points. Since I have 150 samples. fft.length should equate to 300. I am then checking the frequency of the bin, and if it is less than 0.33333334 I am using it, otherwise, it is getting zeroed out.

When I plot my graph, the odd part is that I get two waves:

If I zoom in on the top one, here is what I see:

I also notice that the numbers on the x axis are increments of 2. For example, the wave where the amplitude is around 400K on the Y axis starts at 1, and has points at 3,5,7,9, etc. The lower one starts at 0, and has points 2,4,6,8,etc.

Question:

Should I be seeing two waves? If yes, how do I combine it into one? If no, what step might have I destroyed?

From my cloudy knowledge, the only thing I can think of that can explain this is the notion of imaginary and real numbers. But, I do not see why they should be shown in a time domain, after an iFFT occurs.

FloatFFT has a funny array layout. Real and Imaginary components of each datapoint/bin have to be interleaved in the data array you provide.
The FFT plot you show it is not a common FFT plot: is the real and imaginary part of the FFT interleaved: $[ Re\{X[0]\}, Im\{X[0]\}, Re\{X[1]\}, Im\{X[1]\}, ... ]$.

Handle two consecutive points everywhere accordingly. For example: $bin[0] = abs( sqrt( a[0]*a[0] + a[1]*a[1] )$ and so on... Also, follow that layout when you perform the inverse.

Documentation is not the best and you will have to read it carefully.

Bonus: the reason you have to provide a $2 \cdot n$ array to some functions is not that the FFT has to double the data set length, it is because each bin is stored in two consecutive array positions.


How to get maxima of rho in hough table i applyed hough transform in my image and this is the result of hough table,  i must get the maxima of rho in two basic angles of theta , the one is around -65 : -60 and the other is round 30:35

How can i get them?
I tried the below code about hough but the houghpeaks function returns all the maxima peaks in the hough table. I want to detect the two basic angles as you can see

[hough,theta,rho] = hough(negfillHoles,'Theta',-90:1:89,'RhoResolution',2);
lenHough=length(hough(:));
maxpeaks=lenHough; %10000;
peaks = houghpeaks(hough,maxpeaks,'threshold',20); %50
x=theta(peaks(:,2));
y=rho(peaks(:,1));
plot(x,y,'o','color','r');

I will really appreciate your help, it is a matter of life or death to find this answer for my home work

Why don't you try to mask your data first. E.g. take the variance for every $\theta$ value, so that you get a $\rho$ variance information for each $\theta$ value. Variance should be high around your peask, since all the other parts are basically black or gray - only in the peak region there is a lot of changes. You could also play around with a auto-correlation of the data along that axis. Then, you could threshold the resulting vector so that you have a better estimate where the peaks are. Then apply the peak detection algorithm only in those areas of $\theta$.


Question about Hilbert transform Hilbert transform of a function $g(t)$ which is defined in time domain, would result in another function in time domain. Is there any other transformation like Hilbert that the results be in time domain?

This question is hard to answer because it depends on your definition of "transform". The Hilbert transform can be written as a convolution with the kernel

$$h(t)=\frac{1}{\pi t}\tag{1}$$

and, consequently, it can be represented by a linear time-invariant (LTI) system with an impulse response given by $(1)$. So the application of any LTI system to a signal can be seen as a transform or mapping which doesn't change the domain of the signal. E.g., the application of an ideal band pass filter to an input signal maps that signal to the space of band-limited functions, and in this sense, it is a transform(ation).


Smoothing "font characters" outlines in bitmap image Adobe Acrobat has this as a feature called ClearScan. While doing OCR it allows image enhancement, by smoothing irregularities on font outlines.
For example, this source:

turns to this:

Are there other tools with such feature, or public algorithms (Matlab/Python) that can smooth contours of arbitrary font-like objects?

I think the open-source command-line program potrace, might do what you want.  It converts bitmaps to bezier curves and has a bunch of options allowing you to trade off smoothness and accuracy.

The open-source Inkscape vector (svg) editor has potrace built in (under the Path->Trace Bitmap menu option.)  The result of applying your example in Inkscape, (I used the mode Inkscape calls "Brightness cutoff" with the default threshold of .45, since your image was already black-and-white rather than grayscale) scaling the size of the resulting bezier curves up a bit, and then having Inkscape output a bitmap is:

One of the cute side-effects of converting to bezier and then back to grayscale is that you get nice anti-aliasing around the curves.

With grayscale input images I find that the "Edge detection" mode works a little better.  There are also options for choosing less or more smoothing, and less or more speckle elimination.


Z-transform of a downsampler In this paper or multirate filtering, the author establishes the following mathematical relationship. Let $y_D$ be the output of a downsampler such that

$$y_D[n] = x[Mn]$$

where $M$ is the downsampling factor. In other words, we keep every $M$-th sample of the original signal. The author then goes on to state the following:

  ...the z-transform of $y_D[n]$ is given by

  $$Y_D[z]=\frac{1}{M}\sum_{k=0}^{M-1}X[z^{1/M}W^k]$$

  where $W^k$ is the $M$-point Discrete Fourier Transform kernel, namely
  $e^{(-j2\pi k)/M}$.

How can we go from the former expression to the latter? What is the relationship between DFT and the Z-transform that allows for such transition?

This derivation is a tricky one. The approach suggested before has a flaw. Let me demonstrate this first; then I will give the correct solution.
We wish to relate the $\mathcal{Z}$-transform of the downsampled signal, $Y_D(z) = \mathcal{Z}\{x[Mn]\}$, to the $\mathcal{Z}$-transform of the original signal $X(z) = \mathcal{Z}\{x[n]\}$.
The wrong way
One could think of simply plugging the expression for the downsampled signal into the expression of the $\mathcal{Z}$-transform:
$Y_D(z) = \sum\limits_{n=-\infty}^{+\infty}{x[Mn] z^{-n}}$
A change of variable $n'=Mn$ seems obvious:
$Y_D(z) = \sum\limits_{n' \in M\mathbb{Z}}{x[n'] z^{-n'/M}}$
However, it is important to realize that even though the new summation index $n'$ still runs from $-\infty$ to $\infty$, the sum is now over 1 out of M integer numbers. In other words,
$n' \in M\mathbb{Z}=\{...,-2M,-M,0,M,2M,...\}$,
while the definition of the $\mathcal{Z}$-transform requires
$n \in \{...,-2,-1,0,1,2,...\}$.
Since this is no longer a $\mathcal{Z}$-transform, we cannot write:
$Y_D(z) = X(z^{1/M})$
The right way
Let us first define a 'helper' impulse train signal $t_M[n]$ as:
$\begin{align*}
t_M[n]
& = \sum\limits_{k=-\infty}^{+\infty}{\delta[n-kM]} \\
& =
\left\{
     \begin{array}{lr}
       1 & : n \in M\mathbb{Z}\\
       0 & : n \notin M\mathbb{Z}
     \end{array}
   \right.
\end{align*}
$
This function is $1$ at one out of every $M$ samples, and zero everywhere else.
Equivalently, the pulse train function can be written as:
$ t_M[n] = \frac{1}{M} \sum\limits_{k=0}^{M-1}{e^{j2\pi k n / M}} $
Proof: We need to consider separately the cases $n \in M\mathbb{Z}$ and $n \notin M\mathbb{Z}$:
$$\begin{align*}
 t_M[n]
& = \frac{1}{M} \sum\limits_{k=0}^{M-1}{e^{j2\pi k n / M}} \\
& =
\left\{
     \begin{array}{lr}
       \frac{1}{M} \sum\limits_{k=0}^{M-1}{1} && : n \in M\mathbb{Z}\\
       \frac{1}{M} \frac{1-e^{j2\pi n}}{1-e^{j2\pi n / M}} && : n \notin M\mathbb{Z}
     \end{array}
   \right. \\
& =
\left\{
     \begin{array}{lr}
       \frac{1}{M} M && : n \in M\mathbb{Z}\\
       \frac{1}{M} \frac{1-1}{1-e^{j2\pi n / M}} && : n \notin M\mathbb{Z}
     \end{array}
   \right. \\
& =
\left\{
     \begin{array}{lr}
       1 && : n \in M\mathbb{Z}\\
       0 && : n \notin M\mathbb{Z}
     \end{array}
   \right.
\end{align*}$$
In the case $n \notin M\mathbb{Z}$, we used the expression for the finite sum of a geometric series.
Now let's go back to our original problem of finding the $\mathcal{Z}$-transform of a downsampler:
$Y_D(z) = \sum\limits_{n=-\infty}^{+\infty}{x[Mn] z^{-n}}$
We apply the substitution $n'=Mn$, keeping in mind that this makes the summation run only over integer multiples of M:
$Y_D(z) = \sum\limits_{n' \in M\mathbb{Z}}{x[n'] z^{-n'/M}}$
We can now use the above impulse train function to safely rewrite this as a summation over all $n \in \mathbb{Z}$:
$Y_D(z) = \sum\limits_{n=-\infty}^{+\infty}{t_M[n] x[n] z^{-n/M}}$
Using the above formulation for the impulse train function as a finite sum of exponentials, we get:
$\begin{align*}
Y_D(z)
&= \sum\limits_{n=-\infty}^{+\infty}{(\frac{1}{M} \sum\limits_{k=0}^{M-1}{e^{j2\pi k n / M}}) x[n] z^{-n/M}} \\
&= \frac{1}{M} \sum\limits_{k=0}^{M-1}{\sum\limits_{n=-\infty}^{+\infty}{e^{j2\pi k n / M} x[n] z^{-n/M}}} \\
&= \frac{1}{M} \sum\limits_{k=0}^{M-1}{\sum\limits_{n=-\infty}^{+\infty}{x[n] (e^{-j2\pi k / M} z^{1/M})^{-n}}}
\end{align*}$
The summation on the right is a summation over all integers, and is therefore a valid $\mathcal{Z}$-transform in terms of $z'=e^{-j2\pi k / M} z^{1/M}$. Therefore, we can write:
$
Y_D(z) = \frac{1}{M} \sum\limits_{k=0}^{M-1}{X(e^{-j2\pi k / M} z^{1/M})}
$
This is the formula for the $\mathcal{Z}$-transform of a downsampler.


Stability of the Low pass filter after Spectral Transformation Using a spectral transformation a filter $G(z^{cap})$ is
designed to be a low pass filter from an existing lowpass filter $G_{lp}(z)$ via $G(z^{cap}) = G(z)z=H(z^{cap})$
.  Which property should H(z) have in order to preserve
stability of the designed filter?

The filter is stable if the all the poles of the transformed filter, i.e., filter obtained after applying spectral transformation, lie inside the unit circle. For example, if $H(z)$ is the original filter and $F(z)$ is the filter response of the transformation filter, then the spectrally transformed filter response is given as
\begin{align}
G(z) = H(F(z)) = H(z)|_{z^{-1}\leftarrow 1/F(z)}
\end{align}
where
\begin{align}
\frac{1}{F(z)} = \frac{1-\beta z}{z - \beta}.
\end{align}

Check if the poles of $G(z)$ lie inside the unit circle. If yes, then $G(z)$ is stable. Note that $F(z)$ is an all-pass filter. Therefore, the regions of stability and instability of $H(z)$ are preserved in $G(z)$.


Given f(x) and its FFT F(u), I need to prove that df/dx = F(u)*2iu(pi)/n Given $f(x)$ and its FFT $F(u)$, I need to prove that $df/dx = F(u) * 2iu\pi /n.$

$df/dx = f'(x)$ and $n$ is the number of pixels of the one dimensional image $f$.

I tried to use the convolution theorem using conv($f[1 -1]$) which is the same as $df/dx$ and in the frequency domain it's $F(u) * FFT([1 -1])$ however there sizes don't match.

This should have been a comment on the question but is too long to fit within
the 500-character limit.

The $f$ that is the subject of this question is a sequence or vector of $n$ numbers
$$[f_0, f_1, \ldots, f_{n-1}] ~~ \text{or} ~~ \left[f[0], f[1], \ldots, f[n-1]\right]
~~ \text{or} ~~ \left[f(0), f(1), \ldots, f(n-1)\right]$$
depending on choice of notation, and its Discrete Fourier Transform (DFT) $F$ (called the
FFT by the OP) is another sequence or vector of $n$ numbers. Thus,
it is not clear what the derivative $\frac{\mathrm d}{\mathrm dx}f(x) = f^\prime(x)$ means
since the argument of $f(\cdot)$ is an integer and not a continuous variable.  One
could attempt to make sense of what is intended by looking at the OP's approach to the
problem via (cyclic or periodic) convolution with the sequence $[1, -1]$ (actually should have been $[1,-1,0,0,\ldots, 0]$ which results in replacing $f_i$ by $f_{i+1}-f_i$.
Now, we are getting closer to the notion of differentiation since the derivative of
a continuous function is obtained by looking at the limit of $\frac{f(x+h)-f(x)}{h}$ as
$h \to 0$. Here we have the closest we can get since we are
using $\frac{f(i+1)-f(i)}{1} = f_{i+1}-f_i$, and of course $1$ cannot "approach" $0$.
Unfortunately, this does not give the desired result.  The $k$-th term of the
DFT of $[1,-1,0,0,\ldots, 0]$
is readily computed from the definition (no FFTs needed)
as $1 - e^{-j2\pi k/n}$, and so the DFT of the "derivative" sequence
$f_{i+1}-f_i$ has $k$-th term $F_k(1 - e^{-j2\pi k/n})$ which is not
what the OP wants to prove.

So, what does the OP want to prove? It is
$\frac{\mathrm df}{\mathrm dx} \mathbf{=} F(u)*j(2\pi u/n)$.
This needs clarification on two grounds. First, the left side
of this equality is a function of $x$, the right side a function
of $u$. So what does it mean to say the two are equal?
Second, what the right side means needs clarification. Does
the $*$ mean multiplication of the two complex numbers $F(u)$ and
$j(2\pi u/n)$, or is the $*$ denoting complex conjugation of $F(u)$
before the multiplication, or does $*$ denote convolution so
that $F(u)*j(2\pi u/n)$ means the $u$-th term of the (cyclic or periodic)
convolution of the sequences $[F_0,F_1,\ldots,F_{n-1}]$ and
$[1, j(2\pi/n),j(2\pi 2/n), \ldots, j(2\pi (n-1)/n)]$?
The closest match occurs if we assume that $=$ is a typographical
error and that what is meant is not the equality that the OP
claims is needed. Instead, suppose that what was really asked
for but mistranslated by the OP is that the $k$-th term of
the DFT of $\frac{\mathrm df}{\mathrm dx}$
by which is meant the sequence with $i$-th term $f_{i+1}-f_i$
is $F_ke^{j2\pi k/n}$ (another typo, the missing exponential!). As
noted above, this is not quite right: it should be
$F_k(1 - e^{-j2\pi k/n})$. Note, incidentally, that the inverse
DFT of the sequence with $k$-th term $F_ke^{j2\pi k/n}$ is a
cyclic shift of the sequence $f$; see, for example,
this
Wikipedia entry

But maybe as a non-expert in image processing I am missing something completely
and so I await clarification from the OP, or perhaps a more detailed
explanation from @JasonR of the method outlined in his comment on the OP's
question.


VOFDM (Vector OFDM) advantages and disadvantages What are the advantages and disadvantages of using VOFDM over OFDM?
I know that one advantage is frequency diversity, what are others?

edit : By VOFDM I mean Xia's method as shown here: arxiv.org/pdf/1507.06833.pdf

Referring to VOFDM as described in this reference, some thoughts on that (I wrote this summary some years ago, due to an issue in my own research). We can understand VOFDM as quite similar to interleaved SC-FDMA, (this means that a single carrier occupies several frequency bins, which are not grouped together but spread around the full bandwidth of the system). So, assuming the system has $N=KM$ OFDM subcarriers (FFT window length N, N frequency bins), the frequency indices occupied by an interleaved SC-FDMA carrier would be $[l, l+K, l+2K, l+(M-1)K]$, where $M$ is the number of frequency bins per SCFDMA carrier and $K$ is the number of different carriers possible in the full bandwidth. $l=0,\dots,K-1$ is the index of the SC-FDMA carrier.  Then, after FFT at the receiver, the received signal for the $l$ carrier is given by

$$Y_l=H_ld_l+n_l$$ where $Y_l$ is the received signal at the frequency bins for the $l$ carrier, $d_l$ is the transmit data for the $l$ carrier and $n_l$ is AWGN. $H_l$ is the equivalent channel matrix, which can be considered a circulant matrix.

For the considerations below, VOFDM is actually doing exactly the same thing (compare eq (4) of this reference). So, how can we compare OFDM and VOFDM?

The big advantage of OFDM is low complexity at the receiver. A single DFT plus one-tap equalization is sufficient for the detection. In contrast, VOFDM requires additional M-point DFTs on each carrier to perform one-tap equalization.
In OFDM, each frequency bin can be treated independent from the others (it's orthogonal). This allows great implementation speedups in terms of parallel processing. VOFDM can be only parallelized to $K$ "threads".
The complexity aspect comes even more into play, when we consider MIMO transmission (e.g. spatial multiplexing). With today's hardware we actually have the chance to use a sphere-decoder in real-time on each carrier to get the optimal receiver performance. For VOFDM, the dimension of the sphere-decoding would be M-times larger, which quickly becomes infeasible for today's hardware.
VOFDM can exploit frequency diversity, since one data symbol is spread around multiple, far apart, frequency bins. However, in order to harvest this diversity, you'd at least need an LMMSE equalizer. However, even better performance can be obtained by employing a sphere-decoder on the equivalent channel $Y_l=H_ld_l+n_l$. But, this requires very high implementation complexity.  Also, it needs to be investigated, if additional (smart) soft-output decoding of the (certainly) encoded input data stream can gain the same/similar diversity in the OFDM case.
In principle VOFDM can achieve better FER performance when considering iterative receivers due to its frequency diversity and concatenation of a channel code and the constellation constraint of the equivalent channel. However, again, to achieve this performance high calculation complexity needs to be accepted.
regarding robustness against synchronization errors and channel estimation errors, I'd assume both systems are similar. However, Channel estimation itself can be more challenging for VOFDM, since no orthogonal frequency bins exist which can be nicely used for dedicated pilot frequency bins.

So, to summarize: VOFDM is a nice idea and has potential to perform better in terms of FER than OFDM. But, it needs a lot more complexity.


select local maxima in energy signal I have the following plot which represents the RMS energy of a signal over time.

The energy goes up and down with varying maximum amplitude. I want to keep the portions of the signal where the energy is quasi-steady (i.e. the top of each "hill", like the brown line). I cannot use thresholding because the local maxima are varying. How can I do it?

I am using Matlab.

I would suggest to apply a median filter like

y=filter(1,0.05*ones(1,20),x)

And then apply Matlab function "findpeaks"

As an alternative to second step you perform a segmentation based on different from zero areas like in

ne0 = find(A~=0);                                   % Nonzero Elements
ix0 = unique([ne0(1) ne0(diff([0 ne0])>1)]);        % Segment Start Indices
ix1 = ne0([find(diff([0 ne0])>1)-1 length(ne0)]);   % Segment End Indices
for k1 = 1:length(ix0)
    section{k1} = A(ix0(k1):ix1(k1));
end
celldisp(section)                                   % Display Results

And then perform calculation on each segment. For example sum all samples and divide by number of samples.


How does time shift correspond to phase change in a discrete signal? I was watching this video where the presenter remarks:

  For a discrete signal, time shift corresponds to phase change in a discrete signal but not vice versa.

I was trying to figure out how this could be proved. For example, if our discrete signal is

$f[n] = A \cos(\Omega_o n + n_o)$ where $ n,n_o \in Z$

If the signal undergoes a time shift $n'$, we have

$ A \cos(\Omega_o (n + n') + n_o) = A \cos(\Omega_o n + \Omega_o n' + n_o)$

This would correspond to a phase shift only if $\Omega_o n'$ is an integer. But that means $\Omega_o$ needs to be an integer, which is not always the case.

What is the correct explanation here ?

If you have a signal

$$f[n]=\cos(\Omega_0n)$$

and you apply a time shift of $n_0$ you get

$$f[n+n_0]=\cos(\Omega_0(n+n_0))=\cos(\Omega_0n+\Omega_0n_0)=\cos(\Omega_0n+\phi)$$

where $\phi=\Omega_0n_0$ is the phase shift.

The other way around, if you have a phase shift of $\phi$, this is not always equivalent to a time shift of the original signal:

$$g[n]=\cos(\Omega_0n+\phi)=\cos(\Omega_0(n+\phi/\Omega_0))$$

which only corresponds to an integer time shift if $\phi/\Omega_0$ is integer, i.e. if $\phi/\Omega_0=n_0$ which results in

$$g[n]=f[n+n_0]$$

EDIT:
Re-reading your question, I think that the misunderstanding lies in the fact that you believe that the phase of a discrete-time signal must be integer. This is not the case. Imagine two continuous-time signals
$$x(t)=\cos(\omega_0t)\quad\textrm{and}\quad y(t)=\cos(\omega_0t+\phi),\quad\phi\in\mathbb{R}$$

Note that the following holds for any value of $\phi$:

$$y(t)=x(t+\phi/\omega_0)$$

So we can always express $y(t)$ as a shifted version of $x(t)$. Now imagine that we construct two discrete-time signals by sampling $x(t)$ and $y(t)$ at times $t_n=nT$ with some real-valued $T>0$:

$$f[n]=x(nT)=\cos(\omega_0nT)=\cos(\Omega_0n)\\
g[n]=y(nT)=\cos(\omega_0nT+\phi)=\cos(\Omega_0n+\phi)
$$
with $\Omega_0=\omega_0T$. Note that the phase $\phi$ is of course the same as before (and for this reason not necessarily integer). The difference between the continuous-time and the discrete-time case is now that $g[n]$ cannot in general be obtained from $f[n]$ by time-shifting because we can only shift by integers. The condition under which $g[n]$ can be obtained by shifting $f[n]$ is if $\phi/\Omega_0$ is integer because then we can write

$$g[n]=\cos(\Omega_0(n+\phi/\Omega_0))=f[n+\phi/\Omega_0]=f[n+n_0],\quad n_0\in\mathbb{Z}$$


The Standard Deviation of The Derivative of a Signal Given a signal with zero mean and a standard deviation of 0.1 sampled at 5000 Hz.
What would be the Standard Deviation of its 1st, 2nd and 'n' derivative?

For instance, let's say we measure the distance to a constant range target.
The STD of the sensor measurements of the range is 0.1 [Meter] (Assume the error is white) and its sampling rate is 5000 Hz.

What would be the standard deviation of the acceleration process?

Thanks.

As pointed out in pichenettes' answer, you need to know more about the original signal $x(t)$. If the signal's power spectrum $S_x(\omega)$ is known then the power spectrum of its $n^{th}$ derivative is

$$S_y(\omega)=\omega^{2n}S_x(\omega)$$

The power of the $n^{th}$ derivative can then be computed as

$$E\{|y(t)|^2\}=\frac{1}{2\pi}\int_{-\infty}^{\infty}\omega^{2n}S_x(\omega)d\omega\tag{1}$$

(assuming convergence of the integral, e.g. by assuming that $x(t)$ is band-limited). Since the mean of the derivative is zero (even if the mean of the original signal isn't), the variance of $y(t)$ equals its power, and so the standard deviation is the square root of the expression in (1).


Slicing and Recognition of book title and author from bookshelf image I am trying, for my own learning purposes, to develop an implementation of an algorithm that would list books, given an image of a bookshelf like the following:

The first step is to slice the image into individual books.

My algorithm, in Mathematica, is:

    img = ColorConvert[Import["http://i.stack.imgur.com/IaLQk.jpg"], "GrayScale"]

do a basic edge detection &
remove the text and try to keep the long lines

edge = DeleteSmallComponents[EdgeDetect[img, 3],Last[ImageDimensions[img]]/5]

then remove unwanted horizontal lines

lines = Sort[ImageLines[img5] /. {{0., _}, {_, _}} -> Sequence[]]
Show[img, Graphics[{Thick, Orange, Line /@ lines}]]

The results, though, are less than great:

My questions are:

How can I improve this to yield better results?
Is there a smarter way to do this?
how further should I process the images to increase accuracy at the (later) OCR phase?
How to use the color information to improve the segmentation?

Here is the link to a research paper that tries to do the same thing as you wanted. It might help you.using image features
Also a cool video on the youtube


Help or suggestions with Line detection in microscopy images I'm doing a personal project in which I want to select some very specific lines in a microscopy image.
The first picture is the raw image, the second picture is the image for the lines (in yellow) I would like to detect. Not straight lines yes, but lines indeed.

Now, I made a very basic code for line detection that kind of "works" detecting lines, but it takes a lot of undesired detected lines in the process. See an example with another image (and the code below it, in Python, using OpenCV).

import cv2 as cv
import numpy as np

#Load and normalize .tif image:
img = 'imagenes.tif'
img = cv.imread(img, 0)
img_scaled = cv.normalize(img,  np.zeros((800, 800)), 0, 255, cv.NORM_MINMAX)

#Resizing images:
pct = 30
w = int(img.shape[1] * pct / 100)
h = int(img.shape[0] * pct / 100)
dim = (w, h)
resized = cv.resize(img_scaled, dim, interpolation = cv.INTER_AREA)

# First stage: Gaussian blur.
gaussian = cv.GaussianBlur(resized,(5, 5),0)
# Second Stage: edge detection with Canny.
low_threshold = 90
high_threshold = 150
edges = cv.Canny(gaussian, low_threshold, high_threshold)
# Third stage: Hough Line transform.
rho = 1
theta = np.pi / 180
threshold = 15
minLineLenght = 15
maxLineGap = 15
myBlankLines = resized.copy()
lines = cv.HoughLinesP(edges, rho, theta, threshold, np.array([]),
                    minLineLenght, maxLineGap)
for line in lines:
    for x1,y1,x2,y2 in line:
        cv.line(myBlankLines,(x1,y1),(x2,y2),(255,0,0),5)

cv.imshow('image_resized', resized)
cv.imshow('Stage 1: Gaussian', gaussian)
cv.imshow('Stage 2: Edges', edges)
cv.imshow('Stage 3: Lines', myBlankLines)
cv.waitKey(0)
cv.destroyAllWindows()

Now. I don't mind building a filter from zero and avoid using OpenCV or any other image processing/computer vision library. But I feel a bit lost on what direction I should go. I'm not sure if:

Doing more pre-filtering.
Adjust parameters for Hough Lines Transform.
Something else related with Machine Learning perhaps.

What would you suggest? Thanks for any heads up you might have on this interesting task.

I used Median Blur instead of Gaussian Blur, and used Sobel Filter to detect the edges, the results are not perfect, but you can improve them if you tuned the parameters well(and they are many!)

#==================
# Import Libraries
#==================
import numpy as np
import cv2
import matplotlib.pyplot as plt

# Read Image
img = cv2.imread('1.png', 0)

# Median Blur
img = cv2.medianBlur(img, 3)

# Sobel filter for Edge detection
sobelx = cv2.Sobel(img,cv2.CV_64F,1,0,ksize=5)
sobely = cv2.Sobel(img,cv2.CV_64F,0,1,ksize=5)
sobel = np.sqrt(sobelx**2+sobely**2)

# Threshold Sobel to locate the lines
maxi = 0.30*sobel.max()
sobel[sobel<maxi] = 0

# Normalize
sobel = (255*sobel) / (sobel.max()-sobel.min())
sobel = sobel.astype('uint8')

# Hough Line transform.
rho = 1
theta = np.pi / 180
threshold = 15
minLineLenght = 18
maxLineGap = 13
myBlankLines = img.copy()
lines = cv2.HoughLinesP(sobel, rho, theta, threshold, np.array([]),
                    minLineLenght, maxLineGap)
for line in lines:
    for x1,y1,x2,y2 in line:
        cv2.line(myBlankLines,(x1,y1),(x2,y2),(255,0,0),1)

# Visualize Results
plt.imshow(myBlankLines, cmap='gray')
plt.axis('off')
plt.show()


In-phase and Quadrature demodulation I'm new with signal processing and have started to learn it for physics research I am involved in. I have a signal which I want to process to find the amplitude of one frequency of. To do this, I'm multiplying the signal numerically by sin(2πft) and cos(2πft), so Y = signal * sin(2πft), X = signal * cos(2πft).  then taking the square root of the sum of these squared for the amplitude and arctan of Y/X for the phase over time. What does the plot of arctan(Y/X) actually show me? For an example, I have a code in Python I'm using for testing this.

from matplotlib.pyplot import *
from numpy import *
from scipy.signal import butter, lfilter

def butter_lowpass_filter(data, cutoff, fs, order=5):
    '''Creates a lowpass filter for the data'''
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False) #filters data
    y = lfilter(b, a, data) #filters data
    return y

def sin_stat(time, f0, method = 'sin'):
    '''Creates sine and cosine data to be used for manipulating data'''
    if method == 'sin':
        return [2.0 * np.sin(2*pi*f0 * t) for t in time] #2 in front to cancel out 1/2 later
    if method == 'cos':
        return [2.0 * np.cos(2*pi*f0 * t) for t in time]

def stats(time, signal,f0,sample_freq, cutoff = 1): #cutoff is in hz
    ''' first creates list of sine and cosine data to numerically multiply all
    data by. Then Makes X, Y, which are data times cosine data and sine data
    respectively. '''
    C = sin_stat(time, f0, method = 'cos')
    S = sin_stat(time, f0, method = 'sin')

    X = [d * c for d,c in zip(signal,C)]
    Y = [d * s for d,s in zip(signal,S)]

    X_low = butter_lowpass_filter(X, cutoff, sample_freq) #lowpass each
    Y_low = butter_lowpass_filter(Y, cutoff, sample_freq)

    amp = [np.sqrt(x**2 + y**2) for x,y in zip(X_low, Y_low)]
    phase = [np.arctan(y/x) for y,x in zip(Y_low,X_low)]

    phase = [2 * (x + np.pi/2) for x in phase]
    phase = np.unwrap(phase)
    phase = [x/2.0 for x in phase]
    return (amp,phase)

sample_freq = 100 #hz
end = 1000
t = linspace(0,end,end*sample_freq)
s = 5*sin(2*pi*10*t + 2) + 2*sin(2*pi*7*t + 1) + sin(2*pi*4*t +.5) + .2*sin(2*pi*2.3*t+.7)
plot(t,s)
show()
close()
amp, phase = stats(t,s,2.3,sample_freq,.5)
plot(t,amp)
show()
close()
plot(t,phase)
show()
close()

The amplitude plot shows the amplitude converge to the amplitude I define in the signal, but the phase plot does not. If my guess for the frequency in my In-phase and Quadrature signals is the same as in the simulated signal with a low-pass cutoff of .5 hz, I receive this:

However, if I change the guess frequency to something not exactly the frequency, say 2.4 Hz, I receive this:

Could someone explain what I am actually seeing? What does the slope of this represent? In my real data, I always have a slope of some sort due to not knowing the exact frequency of what I am looking for, only an approximation. I would like to calculate frequency vs time to perform an allan deviation plot, but am unsure of how to calculate frequency from this. I feel I am on the right track, but am a little confused. I've been using this for understanding a little of what I'm doing: http://www.thinksrs.com/downloads/PDFs/ApplicationNotes/AboutLIAs.pdf

Changes in frequency are equivalent to a linear phase sweep, so you're seeing exactly what you'd expect. Consider a reference sin wave with frequency $\omega_1$, and a signal with frequency $\omega_2=\omega_1+\Delta\omega$. Well,
\begin{align}
\sin\omega_2 t &= \sin((\omega_1 + \Delta\omega)t) \\
 &= \sin(\omega_1t+\Delta\omega t) \\
 &= \sin(\omega_1 t+\phi(t))
\end{align}
Where $\phi(t) = \Delta\omega t$. This is the line you're seeing in the unwrapped phase.

So if you do lock in at $\omega_1$, your're going to pick up that frequency component, but the frequency offset will be seen in the phase varying linearly with time.


How do I implement a nonorthogonal quadratic spline wavelet into discrete wavelet transform? I am trying to use the discrete wavelet transform for signal processing (time series data from a plant). I would like to use a mother wavelet that is not in Matlab. Mallat and Zhong 1992 described a nonorthogonal quadratic spline wavelet I would like to use. They name some coefficients and impulse responses, but I'm not sure how to use these (reference below). To detect and treat abnormal sudden changes, the wavelet must satisfy a very particular property. We call a smoothing function any function θ(x) whose integral is equal to 1 and that converges to 0 at infinity. This function must be twice differentiable and define ψ_1 the first order derivative and ψ_2 the second order derivative.

ψ_1=(dθ(x))/dx  and ψ_2=(d^2 θ(x))/(dx^2 )

The wavelet transform is computed by convolving the signal function f with a wavelet function (dilated at a scale s).

W_s^1 f(x)=fψ_(1,s) (x) and W_s^2 f(x)=fψ_(2,s) (x)

The extrema points of the first wavelet transform (WT) gives us the zero crossings points of the second WT. The extrema of W_s^1 f(x) and the zero-crossings of W_s^2 f(x) are the inflection points or the sharp variations points of f(x) smoothed by θ(x).Those edges in the signal can only be found if the wavelet is the second derivative of a smooth function. Furthermore, the first- and second-order WT are used in the detection of steady-state periods.

So I would like to code a quadratic spline wavelet function and a cubic spline scaling function like those in Mallat and Zhong (1992) article. Could anyone help me ?

Mallat, S., & Zhong, S. (1992). Characterization of Signals From Multiscale Edges. IEEE Transactions on Pattern Analysis and Machine Intelligence, 14(7), 710-732. doi:10.1109/34.142909

I would think it is available in the Matlab codes from the Wavelab toolbox. It hosts codes for papers and books. Especially, there is a directory \Wavelab\Books\WaveTour\ for Stéphane Mallat's book: A Wavelet Tour of Signal Processing. You can check images chapter by chapter from
A Wavelet Tour/Images. Those images can be regenerared. For instance, from Chapter 5,  Wavelab codes in \Matlab_Wavelab\Books\WaveTour\WTCh05 say from instance in wt05fig04.m:

  disp('The dyadic wavelet is a quadratic spline with one vanishing
  moment.') disp('These graphs display the modulus square of its Fourier
  transform') disp('dilated by 2^j for 0

which may suit your needs.


Duality in the discrete-time Fourier series Suppose $g[n]$ is periodic with fundamental period $N$ and $f[k]$ being its Fourier coefficients i.e. $$ f[k] = \frac{1}{N}\sum_{n=<N>}g[n]e^{-j\frac{2\pi}{N}nk}$$ In more convenient notation $$g[n]\leftrightarrow f[k]$$  It's easy to prove $$f[n]\leftrightarrow \frac{1}{N}g[-k]$$ This is called "Duality in the discrete-time Fourier series". My problem is applying duality in the case of $$x[n-n_0] \leftrightarrow a_ke^{-jk(\frac{2\pi}{N})n_0}$$ And $$\sum_{r=<N>}x[r]y[n-r]\leftrightarrow Na_kb_k$$
I think for the first case duality should gives$$x[n]e^{jk(\frac{2\pi}{N})m} \leftrightarrow \frac{1}{N}a_{m-k} \tag{1}\label{*}$$ and the second gives $$x[n]y[n] \leftrightarrow \frac{1}{N^2} \sum_{l = <N>}a_lb_{-k-l} \tag{2}\label{**}$$
Edit:
Let $g[n] = \sum_{r=<N>}x[r]y[n-r]$ and $f[k] = Na_kb_k$. So according to the mentioned property $$Nx[n]y[n] \leftrightarrow \frac{1}{N}\sum_{r=<N>}a_rb_{-k-r}$$
Dividing by $N$ leads to $\eqref{**}$.
Similarly if $g[n] = x[n-n_0]$ and $f[k] = a_ke^{-jk(\frac{2\pi}{N})n_0}$ I think then $$x[n]e^{-jn(\frac{2\pi}{N})n_0}\leftrightarrow a_{-k - n_0}$$
Let $m = -n_0$ and we have \eqref{*}.
Obviously these results are wrong but I don't know what's my mistake.
This is taken from Signals and systems by Alan V. Oppenheim:

Based on the references and the definitions given, here is a proof.
Given:
$$x[n] \leftrightarrow a[k], \ and$$
$$y[n] \leftrightarrow b[k]$$
Fourier Series expansion of $x[n]e^{j\frac{2\pi}{N}nm}$ would be derived as follows:
$$\frac{1}{N}\sum^{N-1}_{n=0}x[n]e^{j\frac{2\pi}{N}nm}e^{-j\frac{2\pi}{N}nk} = \frac{1}{N}\sum^{N-1}_{n=0}x[n]e^{-j\frac{2\pi}{N}(k-m)n} = a[k-m]$$
The factor $\frac{1}{N}$ gets consumed in the definition of Fourier Series coefficients. Therefore, the following relationship can be stated:
$$x[n]e^{j\frac{2\pi}{N}nm} \leftrightarrow a[k-m]$$
Similarly, for the second case of multiplication in time-domain, you need to figure out that the inverse Fourier Series relation or the synthesis expression will be the following:
$$x[n] = \sum^{N-1}_{k=0}a[k]e^{j\frac{2\pi}{N}kn}$$
I am not adding the proof of the above which can be easily derived by plugging the expression for $a[k] = \frac{1}{N}\sum^{N-1}_{n=0}x[n]e^{-j\frac{2\pi}{N}kn}$ on the RHS.
Then all we are required to do is proving that the inverse Fourier Transform of convolution of $a[k]$ and $b[k]$ would be product $x[n]y[n]$. Which is as follows:
$$\mathcal{FS}^{-1}(a*b)[k] = \sum^{N-1}_{k=0}(a*b)[k]e^{j\frac{2\pi}{N}nk} = \sum^{N-1}_{k=0}( \sum^{N-1}_{l=0}a[l]b[k-l]) e^{j\frac{2\pi}{N}nk}$$
$$=>\mathcal{FS}^{-1}(a*b)[k] =  (\sum^{N-1}_{l=0}a[l]e^{j\frac{2\pi}{N}ln})(\sum^{N-1}_{k=0}b[k-l]e^{j\frac{2\pi}{N}(k-l)n})$$
$$=> \mathcal{FS}^{-1}(a*b)[k]  = \mathcal{FS}^{-1}a[k] \mathcal{FS}^{-1}b[k] = x[n]y[n]$$
I have replaced $m = (k-l)$ in the proof above to get $y[n]$.
Hence both the equations (5.69) and (5.71) are true and correctly given in the book.
I think you are trying to use equation (5.67) to prove equation (5.69) and (5.71). That is not correct, because (5.67) is stating Duality of Fourier Series expression between Continuous Time and Discrete Time. Equation (5.69) and equation (5.71) are duality of Shift-Property and Convolution Property. The second reason you cannot use Duality Expression (5.67) is because (5.67) states that if $\mathcal{FS}g[n] = f[k]$, then $\mathcal{FS}f[n] = \frac{1}{N}g[−k]$. That is it takes Fourier Series of the Fourier Series itself assuming $f[n]$ in time-domain. But in (5.69) and (5.71), the time-domain sequences remain $x[n]$ and $y[n]$ only, so (5.67) does not apply here. There are 3 Duality Properties stated in equation (5.67), (5.69) and (5.71). They are not the consequence of one another. They are a consequence of similarity in definition of Continuous time Fourier Series and Discrete Time Fourier Series.


What is the default carrier frequency of matlab's pskmod? I'm using the matlab function pskmod for modulating a bit string. What is the carrier frequency used by default by the function?

The function pskmod.m returns a vector of complex baseband symbols. Each element of that vector has the form $y_k=e^{j(2\pi m(k)/M+\phi) }$, where $M$ is the alphabet size, $m(k)$ is an integer between $0$ and $M-1$ which depends on the message symbol at time $k$, and $\phi$ is some initial phase (which defaults to $0$).

There is no carrier involved, and for simulation purposes it's usually also not necessary to use one because a digital communication system can be simulated by using the equivalent baseband system.


DSP - Sinewave: what is $2\pi f t $? I am a beginner in DSP and is familiar with $\sin(\theta)$ declaration (as studied in trigonometry) - where $\theta$ is the angle it makes with X - axis. Please tell me what does $\sin(2\pi f t  + \phi)$ equation means and how $\theta = 2\pi ft + \phi$.

Thanks in advance

A 1D sinusoidal wave, time dependent, is described by the next formula:

$$
y(t) = A\sin(2\pi f t + \phi)
$$

The trick here is just the concept. When you use the formula $A\sin(\theta) $ you aren't specifying if $\theta$ is time-dependent and you are not thinking in the formula as a wave, but as a mere function. But when you use the argument $2\pi f t$ you are considering the angular frequency $2\pi f$ ($rad/s$), where $ f $ is the normal frequency, and the wave's phase $\phi$ (which represents a shift in the wave).

If you are working with a simple sine wave, you can assume a zero-shift, $\phi = 0$, and select a working frequency. Sample the wave and work normally with your DSP.

I hope it helps!

You can check sine wave description in WIKI


understanding image error metrics, mean squared error i am trying to understand the error metrics for images. i had no statistics class or anything and i'm having a bit of problem with even the simplest things.

i have started reading this http://unitstep.net/downloads/school/image-qa-presentation/

just for starters, it says that mean squared error is (page 3)

sum_over_i_j( |ImageA(i,j)-ImageB(i,j)|^2 ) / amount_of_pixels

i understand that this is trying to be something like euclidean distance between two points, but why is there no square root then?

then when you go to the next page of the presentation it tells you that this is a special case of minkowski distance and there IS the square root included in case of mean squared error!

and then, shouldn't we calculate distance of two points by summing the powers of the differences of their dimension components? (`[ (x1-x2)^2 + (y1-y2)^2 + ...] ) and only then square root the sum? then do it for all the pixels and average it.

while on the page 4 it says that we power, sum and square root the pixels first and then do this for all the channels (dimensions) separately?

i don't understand :(

First of all, the mean squared error should not have a square root. As the name implies, it is the mean of the squared difference. What you are referring to is the root mean squared error. Second, the this is a pixel-wise comparison between to images. Therefore the number of pixels is the dimension, and the sum is performed on it. There is even a comment on that page about color images.


How to go about implementing a program to measure the "intensity"/"saturation" of a color in an image? I have a task of developing an application to measure the intensity of a color in an image. I am researching about how to go about it.

  Intensity refers to the purity of a hue. Intensity is also known as
  Chroma or Saturation. The highest intensity or purity of a hue is the
  hue as it appears in the spectrum or on the color wheel.
  Source

The image will be RGB but I can convert it HSV or any other color space.

I am not asking for code, but just high level steps on how to do it.

I am thinking of using OpenCV library for Java.

You can find mapping from RGB space to HSV here:
https://math.stackexchange.com/questions/556341/rgb-to-hsv-color-conversion-algorithm

For other spaces, you can refer to vision modelling textbooks, i.e: Digital video quality, Stefan Winkler


Soft question: seeking online course reference for CS Is there a good online course (with video) for Image/Signal Processing or Matrix Computations that includes a section on the principles of Compressed Sensing?

EDIT: There is no shortage of links to go through:
http://www.quora.com/Compressed-Sensing/What-are-some-good-resources-for-learning-about-compressed-sensing

If anyone can suggest a good one from experience I would appreciate it.


How does carrier frequency offset effect the constellation diagram in OFDM? In the case of OFDM modulation, if the carrier frequency is slightly off  then how does that effect the constellation diagram ?

Is there some general rule how it will change (such as rotate) ?

Or is it more or less arbitrary ?

so let's act as if we had only one active subcarrier "at a time", just to ignore interference between the different subcarriers.

The frequency offset in time domain means a time offset in frequency domain, and probably a fractional one at that. Since that offset is constant for each carrier, the resulting sample stream (i.e. the output of your DFT) is simply shifted in time.

Now comes two things:

A frequency offset usually means that your sampling rates don't match perfectly, either. In any case, this means you're losing orthogonality of the OFDM subcarriers – which is why OFDM synchronizers need to be more accurate than (non-sinc-shaped) filterbank multicarrier systems, simply because each DFT bin "sees" energy of its adjacent, slightly shifted, neighbor.

That's why OFDM systems typically make use of a proper preamble – if you enjoy reading classics, go for Schmidl and Cox' Robust frequency and timing synchronization for OFDM (Google tells me you can get the PDF from here) – in fact, that paper goes from understanding how to timing-synchronize an OFDM system to fixing the carrier offset (because they are so closely linked together).


Seperation of wideband and narrowband - Adaptive Filter I have the following diagram for the adaptive seperation of a narrowband and wideband signal using LMS algorithim,

The way it was explained to us was that "The narrowband signal is correlatied over significant time whereas the broadband signal is not correlated over significant time" and in the delayed line, the ADF can only predict the narrowband.

That doesn't make much sense to me. Why is the narrowband correlated over significant time and the broadband isn't and also what's the point of the delay in this?

A narrowband signal seems like (almost) periodic as indicated by

$$ x[n] = m[n] \sin( w_0 n) $$ where the message $m[n]$ has such a low bandwidth that the peak amplitude (the envelope) of the carrier sine wave changes very slowly compared to how fast the sine wave oscillates between those +/- envelope limits. This makes its autocorrelation sequence also to be long tailed and almost periodic. And it means that the narrowband signal is highly correlated for long lags (long delays).

The opposite is true for the wideband signal, its correlation sequence ceases very rapidly even for a small amount of lag.

The Wiener filter works by estimating the supplied desired response signal from its filter input signal, whose success depends on the amount of correlation between the desired signal and the filter input.

If there is large correlation between the desired signal and the filter input, then the filter will successfully estimate the desired signal and the error (which is desired signal minus the filter ouyput) signal will be small. On the other hand if there is not enough correlation between the two, then the estimation will be poor.

Based on this operational principle of the Wiener filter, and based on the fact that the narowband and the wideband signals have different correlation characteristics, from the given adaptive filter scheme, you can explain its operation as follows.

First, you see that the filter input signal and the desired response signals are the same (except the delay) and are a mixture of an auto-correlated narrowband component plus an uncorrelated wideband component.

The delay is essentially used to decorrelate the desired response and the filter input. However, by carefuly chosing the amount of delay (the lag), you can insure that the wideband component is fully decorraleted whilst the narrowband component is still corelated enough. If you chose the delay to too long, both components will lose their cross-correlation to the desired response. If you chose the delay too short, you would not be able to suppress the wideband component enough. Some intermediate delay will be sufficiently decorrelating the wideband component while still maintaining the cross-correlation between the narrowband components at the desired signal and the filter input.

When this is the case, the filter will successfully predict the correlated narrowband component in the desired response from the delayed version of it at its input, and it will therefore be subtracted from the supplied desired response signal which yields as the narrowband component as the error signal.


z transform property question I have a few questions about this assignment I am solving.
It is q3.33(a) from oppenheim's DSP book 2nd edition.

basically, the question gives you a pole zero plot, from which I was able to determine the correct x[n]. (Couple classmates confirmed my answer on that).

and it also gives you the y[n] of the form $(0.5)^n \cdot x[n]$.

And wants you to work out a pole-zero of y[n] in the z plane.

my question is this: given that I already know X[z], can i just do z transform of $0.5^n$ as if it was $a^n \cdot u[n]$?

Then I will have $[\frac{1}{1-0.5z^{-1}}] \cdot X[z]$

Thanks everyone.

HINT:

Use the property

$$z_0^nx[n]\Longleftrightarrow X(z/z_0)$$

with $z_0=1/2$. You can find this property in Table 3.2 on page 126. It's quite straightforward to show this property by using the definition of the $\mathcal{Z}$-transform.


Role of Power Metric and LPF in AGC(Automatic Gain Controller) I am trying to understand & implement simple toy AGC system . Few questions while reading about it,

1) How POWER metric came into picture, why not other metric like absolute value to control amplitude of input signal?

2) What's the role of LPF(Low Pass Filter) in AGC? LPF is an attenuator or blocker to high frequency components. How such filtering effect can control input amplitudes?

Any reading references will be helpful.

Either power or rms value is typically used for AGC as a peak value would be much more variable given it is an instantaneous measurement. The low pass filter typically refers to the Loop Filter (and if so, would be better described as such as it more typically done as an integrator) which is a core part of the control loop design. A band pass filter prior to the power detector is also important since any signals within the bandwidth of the filter/detector would control the AGC level.

Below are example block diagrams depicting typical AGC implementations in the analog and digital domain.  The log and alog interfaces for the digital implementation provide consistent loop bandwidth over the operational range of the AGC due to using multipliers for the gain adjust.

Further AGC considerations for digital implementation are covered in this post: How to adjust receiver gains to avoid saturation and quantization noise to optimise post digital processing?


Are SOS structures better for linear phase FIR filters? I've seen many discussions suggesting that SOS's are often/usually the preferred structure for fixed point hardware filter implementations, since they preserve conjugate pairs, etc., but I've never seen any specific arguments explaining why they are better for FIR filters.

Do these advantages still apply for direct form, linear phase FIR filters which allow the number of multipliers to be reduced by a factor of 2?

Finally, are there any good resources discussing these trade-offs?  I'm trying to implement fairly long FIR filters with small multipliers, i.e. 40+ taps with 8-10b coefficients, and am looking for resources that might shorten this effort.

well, i see little value in factoring an FIR into quadratics and implementing them in cascade.  it won't be computationally cheaper and it's not better from a quantization noise POV if your FIR has access to an accumulator that is double width.

perhaps it will help with limiting coefficient range so that it's less likely that you have non-zero coefficients (at the tails of the FIR) with values of, say, $10^{-8}$ the size of other coefficients (in the middle).  there are other, more efficient ways of dealing with that problem (like computing the FIR from the outside in, doing the tails first and the middle last).


Gain at given frequency from $z$-plane zero-pole plot. Two methods gives different results I have two zeros at $z=-1$ and two complex conjugate poles at $z=A\cos\theta\pm jA\sin\theta$

This gives me the next transfer function

$$H(z)=\frac{1+2z^{-1}+z^{-2}}{1-2A\cos\theta z^{-1}+A^2z^{-2}}$$

To get filter's DC gain I'm trying to substitute $z$ with $1$:

$$G=\frac{4}{1-2A\cos\theta+A^2}$$

Following graphical method described here (gain as a product of lengths of vectors from every zero to given point divided by product of lengths of vectors from every pole to given point) I've got next formula:

$$G=\frac{2}{\sqrt{1-2A\cos\theta+A^2}}$$

Where is my mistake? Why I've lost square root in the first case?

I think you lost the square root somewhere in your graphical method: Considering the nominator (product of lengths of vectors from every zero given point), I end up with

$$N=2*2$$

i.e. you have two vectors of length 2 (two zeros at z=-1).

Then, for the denominator you have

$$
D = \sqrt{(A\sin\theta)^2+(A\cos\theta-1)^2}^2=A^2-2A\cos\theta+1
$$

And finally $G=N/D$ yields the same result as your analytical solution.


Detection of Troughs and Notches in a PPG Signal I'm working a project which tries to determine Blood Pressure from PPG signals.I'm trying to extract the features as shown below...

I'm having problem in finding the troughs and dicrotic notch for each cycle and also the cycle start and end duration. Are there any manual method methods which can be applied ( I'm a newbie to signal processing )thereby which we can extract the trough and notch location in each cycle ?

Much thanks in advance

Have you tried finding where the derivative of your signal is 0 (or crosses between positive and negative values)?

Simple, approximate, derivative filters are the first difference:

$y[n] = x[n] - x[n-1]$

and the central difference:

$y[n] = x[n] - x[n-2]$

You'll have to compensate for filter delay when mapping the zero crossings of the derivative back to the local minima and maxima of the original signal.
For the first difference, the output is delayed by about half a sample period, and for the central difference, the output is delayed about a sample period.


Calculating Cross Power Spectral density between two complex signals I have a small misunderstanding of cross power spectral density between two complex signals. I know that it is the Fourier transform of cross correlation between two signals. Let's say the complex signals are of the form (a+ib) and (c+id). Then the cross-correlation should result in something like

The real part result => ac - bd i.e. real1*rea12 - imag1*imag2;

complex part result => ad + bc i.e. real1*imag2 + real2 * imag1;

where real1, real2 represents the real part of both the signals and
      imag1, imag2 represents the imaginary part of both the signals.

But my professor said that the cross-correlation should result in

The real part result => ac + bd i.e. real1*rea12 + imag1*imag2;

complex part result => ad - bc i.e. real1*imag2 - real2 * imag1;

Where I am going wrong? Is my understanding of the equation form of cross-correlation correct? And also, I just want to know how to calculate cross PSD from this result. Can I just compute it like

Cross PSD result = |real_part result + imaginary_part result|.^2 ?

or Cross PSD result = |real_part result|.^2 + |imaginary_part result|.^2 ?

Cross correlation is defined as $R_{xy}(\tau) = E\{x(t)y^*(t+\tau)\}$. You missed to take the conjugate of $x_2$. So $x_1x_2^* = (c+id)(a-ib) = ac+bd +i(ad-bc)$.

Cross PSD is the fourier transform of $R_{xy}(\tau)$, not the two equations you mentioned.


Difference between Uniform Linear Array (ULA) 3 dB beamwidth and bearing resolution Is there  a difference between linear array 3dB beam width and linear array bearing resolution ?

The answer will be NO, but to be precise especially in case of conventional beam forming.

Theoretical Derivation

The antenna electric field pattern of array antenna can be given by

where N is the number of antenna elements

d is the spacing between antenna elements

In order to find the beamwidth (3 dB), the above equation should be equated to
  and solve for

The solution will  come to be as

where D is the total aperture distance and can be approximated as

For spacing of   the equation simplifies and can be approximated as

Thus beam width of planar antenna can be represented as

Angular Resolution Using FFT Over Antenna Elements

In order the prove the lemma that the angular resolution obtained by performing the FFT across antenna dimension equals to the beamwidth of the antenna array, we have to obtain the angular resolution using FFT.

The below diagram shows the frequency obtained due to path difference between the antenna elements which occurred due to the angle of arrival other than the broadside angle

The frequency resolution using FFT can be represented as  where N is the number of samples but in our case it is equal to number of antenna elements.

The angular resolution   can be found by equating the difference in frequency of different angles
 and

The frequency resolution becomes:

Solving we get

and if we put  and
we get the same resolution as the beamwidth i.e

Summary

Thus for summarizing the above discussion using FFT we can only achieve the max angular resolution equals to the beamwidth of the antenna array which is equal to angular resolution = antenna array beamwidth = 2/N

Kindly note that this equality is true for the case of conventional beam forming only, there exist also some advance beam former where your angular resolution is better than the 3db beam width

This is the same answer which I provided on
ULA beamwidth and angular resolution, might this will be helpful.


Smooth signal with Hanning I have an acceleration record and i want to smooth it's spectrum by using a hanning window with Bandwidth=0.5 Hz.

How could i do that?

Example in python/matlab will be appreciated!

You can generate Hann window of size $N$ and use it to window your signal $x$ using the MATLAB command

w=hann(N);
x_windowed = x.*w;

Since you specified BW of 0.5Hz, you need normalize it with your sampling rate $F_s$. That is $w_s = \pi 0.5/F_s$. You need to choose $N$ such that your main lobe width of the Hann window falls at $w_s$. It is not an easy problem and you may have to do trial and error initially. You can use freqz(w) command in MATLAB to see the frequency plot of Hann window. You can hover an see at what normalized frequency your main lobe width is.


Chopping up and averaging a signal I have a 1D signal that is roughly periodical (but not exactly) - think of some kind of motion that someone performs repeatedly, and which is measured by some sensor.

I am looking for a way to

automatically chop up the signal into the individual cycles, and
somehow calculate an "average" that preserves the characteristics of the measured motion

How could this be done?

What you need to achieve is to estimate the typical periodicity of a cycle. What is required is to do harmonic analysis on the same.

There are few things - i am making assumption which is critical to your solution.
1. Signal you have is not a fixed analytical expression but some arbitrary format.

We don't quite know a prior what could be the range of periodicity.
It is periodic but not perfectly - so we should rather treat it is a random process rather than a pure analytical function. There is probably noise at play as well - but we will worry about it a bit later.
We need to implement such a system - so we are more interested in discrete time representation of signal and a device to detect such periodicity rather than analytic solution of the given arbitrary function in continuous time.

To describe the solution - we shall start the other way and include each assumption.

a. Case for perfect periodic signal.
Suppose if you have a pure sine wave of a frequency say Omega - W. Since it is expressed as an infinite time continuous sequence, the frequency domain is a finite spectrum - which is pure impulse. So If we can observe the signal in a frequency domain - the position of the said impulse is nothing but the periodicity of the wave form. That's quite easy.

When the signal is sampled (as it would be in real world) as opposed to an analytical form - the resultant frequency response is Discrete Time Fourier Transform - where one see continuous lobes rather than impulses in frequency domain. Here the dominant and subordinates lobes are equivalent to the frequency impulses earlier.

b. what if the signal is not analytic form but arbitrary
Now consider a case of the square wave or any such similar function. For a square wave what you will see is that now you will have many more harmonics. But you will notice one important thing. The primary harmonic you will observe will be also the one which has highest energy - and THAT is also the periodicity. This is usually known as Dominant harmonic. For any arbitrary signal we need to consider DOMINANT harmonic which is the best approximation of the frequency.

c. when signal is not perfectly periodic
When the signal is not perfectly periodic, the observation in frequency spectrum is time varying. In such cases, the primary dominant frequency component will fluctuate around the most expected periodicity. If the periodicity is firm - only shape varies - the location of the dominant component will be static but magnitude of all critical composition will vary. If the periodicity also fluctuates - the same fluctuation will be visible in the position of the frequency component and consequently one needs to "track" such a change over time which is the information we are seeking essentially.

So looking at a, b, and c - the conclusion is simpler - find the dominant principle component (and track it).

d. when the signal is finite time
When the signal is finite in time - the spectrum is infinite - in such cases one usually cannot hope to plot the infinitely long frequency spectrum. However, typically the approaches with finding Segment wise is to have apply the STFT domain

Now here is the catch. What should be the Finite time window for which Finite frequency spectrum should be calculated? The simple answer is the larger the periodicity - even larger needs to be the size of the window. This is equivalent of Nyquist criteria but not exactly same. In order to capture the real periodicity - multiple dominant cycles must be visible in the given Time window span - to ensure that frequency of measured accurately. However, this means that period for which periodicity is measured is larger than a few cycles at least; but this means resolution of change estimate is lesser. Higher the period higher the accuracy - lesser the resolution of change! This is known as uncertainty principle.

Of course, a large number of approaches to do pitch detection from a speech signal does exactly this - and most tools use STFT.

Over time- of course, more modern way of doing things are based on wavelet based approaches. Wavelets doesn't eliminate the uncertainty principle but just hope to get you closest combination of resolution and accuracy.

Dipan.


Face recognition using independent component analysis (ICA) Reading about independent component analysis (ICA), I learned that one of its applications is face recognition. I think in this problem we have a database of images and a test image to be recognized. However, I can't figure out what the independent components are and what the mixture (multivariate signal) is.

For your data set of images, first vectorize the images by raster scanning them, and making them vectors. Thus, say you have $M$ images, each of size 64*64 pixels. Then the total number of pixels per image is $N=64^2$, which means $N=4096$. Now, you have an image matrix of size $M\text{x}N$.

For this image matrix, what you want to do is find the $M$ independent components, of again, length $N$. In other words, you want to decompose this set of images into another set of images, except those images make up the independent components of the original faces. (Similar to eigen-faces, but different).

To answer your question about the mixture, the set of faces that you start with is the mixture. In other words, the set of faces that you have, is assumed to be a mixture of independent faces that you are trying to find.

ICA will find those independent features for you, and for natural images, independent components turn out to be features with extremely high kurtosis, namely, the edges.


How to Combine / Cascade two `3 x 3` Filters into One `5 x 5` Filter If you have a uniform 3x3 box filter T which is:

1 1 1

1 1 1

1 1 1

And an 3x3 laplacian filter W which is:

1 -2 1

-2 4 -2

1 -2 1

Can these 2 filters be combined into one 5x5 filter, the order being first T and then W ? (The combination should be possible) And what is the math behind this combination. I can make a 5x5 filter of all one's and after that stride the W filter over that 5x5 and add everything up, is it that simple?

Let's say you have a convolution kernel $ f $ and another convolution kernel $ g $.
We also have an input signal (1D or 2D) $ x $ and we are after the result of the cascaded convolution:

$$ y = g \ast \left( f \ast x \right) $$

The nice thing about convolution is its associativity property.
Which means:

$$ y = g \ast \left( f \ast x \right) = \left( g \ast f \right) \ast x = h \ast x $$

So the equivalent is $ h = g \ast f $ which with the commutative property of convolution means:

$$ h = g \ast f = f \ast g $$

Pay attention that for the discrete case only the Full Convolution (In MATLAB formulation) is commutative.

So in the above:

mF = [1, 1, 1; 1, 1, 1; 1, 1, 1];
mG = [1, -2, 1; -2, 4, -2; 1, -2, 1];

mH = conv2(mF, mG, 'full')
mH = conv2(mG, mF, 'full')

mH =

     1    -1     0    -1     1
    -1     1     0     1    -1
     0     0     0     0     0
    -1     1     0     1    -1
     1    -1     0    -1     1

mH =

     1    -1     0    -1     1
    -1     1     0     1    -1
     0     0     0     0     0
    -1     1     0     1    -1
     1    -1     0    -1     1


Resampling FIR coefficients I have received from a client a 32-coefficient FIR filter running at 1 kHz. I would like to adapt the filter to a 24-coefficient FIR filter running at 750 Hz while preserving the 0-375 Hz frequency response as much as possible, especially the 0-30 Hz band. Is there a way to resample my FIR coefficients?
Edit :
The FIR has linear phase, is a low-pass filter, and also has zeroes at specific frequencies (multiples of 60 Hz)

As mentioned in a comment, I would probably try to re-design the filter at the new sampling rate. Take equidistant samples of the magnitude of the existing filter between DC and the new Nyquist frequency as a desired response. Since you want to make sure that the new filter has zeros at integer multiples of $60$ Hertz, split your new filter response into two parts:
$$H(z)=P(z)G(z)\tag{1}$$
where $P(z)$ is a polynomial with zeros at integer multiples of $60$ Hz. Your desired magnitude is then defined by equidistant samples of
$$M_D(e^{j\omega})=\left|\frac{H_D(e^{j\omega})}{P(e^{j\omega})}\right|\tag{2}$$
where $H_D(z)$ is the transfer function of the original filter. In your frequency grid, avoid frequencies at which the zeros of $P(z)$ occur. Of course, those zeros are cancelled by the zeros of $H_D(z)$, but you might nevertheless have some numerical problems otherwise.
Now you find a linear phase transfer function $G(z)$ approximating $M_D(e^{j\omega})$ on the unit circle, and your final filter transfer function is given by Eq. $(1)$.
I would use a weighted least squares approximation, which just requires the solution of a system of linear equations. If the range between DC and $30$ Hz is especially important, you can assign those frequencies a higher weight, such that the approximation is better in that range (at the cost of the approximation outside that range).


How to improve the FIR filter performance degradation when the amplitude of the input signal decreases? I'm tring to design a filter to restore the original signal from strong noise (low SNR). I have tried low pass FIR and it could work well as the amplitude of orignial signal is much higher than the background noise. My matlab code of a demo are as follows:
clear
clc
close all

time = 1:10000;
fre1 = 50;
fs = 8000;
amplitude = 100;      % amplitude of the noise

sin_signal1 = sin(2*pi*fre1*time/fs)*50;
noise = rand(1,length(sin_signal1));
noise = amplitude * sqrt(10^(-3/10))*noise;  %noise

sin_signal = sin_signal1+ noise;
figure()
plot(sin_signal);

fir_low = fir1(128,55/8000);
sin_after_filter = filter(fir_low,1,sin_signal);
figure()
plot(sin_after_filter)

When the amplitude of noise is low (eg. amplitude = 0.5), namely a high SNR, the low pass FIR works very well. But when the amplitude of noise is high (eg. amplitude = 150), namely a low SNR, the low pass FIR works very bad.
For a well-designed filter, only changing the amplitude of noise, but not the frequency, why will the performance of the filter degrade? Can't even filter the signal? How to improve the problem? Thanks!

A linear filter just attenuate some frequencies more and others less. If you alter the SNR, the output of the FIR filter should change accordingly.
Many applications use more elaborate noise reduction methods in stead of (or in addition to) linear filtering. If you know that your signal is going to be a sine of known frequency, you may estimate its phase and amplitude embedded in noise using a matched filter. If you know that your noise is shot-noise, you could try to exploit that.
You noise seems to be uniform (non-zero mean) noise. Did you intend to call the normal distributed randn()?
-k


Least Angle Regression (LARS) without Matrix Inversion Sorry if this is too damned long. I did what I could to abbreviate it.

The question is about Least Angle Regression (LARS).

I'm new to numerical work with matrices.

I believe I have a way to compute Least Angle without explicit matrix inversion. I'm hoping someone with more experience can check my work to see if I'm doing anything dangerous, illegal, or otherwise ignorant.

This paper: Least Angle Regression, says:

  "...we compute $m−k$ inner products $c_{jk}$ of the non-active $x_j$
  with the current residuals to identify the next active variable, and
  then invert the $k \times  k$ matrix $\mathcal{G}_k = X'_kX_k$ to
  find the next LARS direction."

(Emphasis mine.)

They update the Cholesky decomposition in each iteration, so it seems they're doing something like this:

At $k$th iteration:
$X_k$ is a matrix formed from $k$ predictors.
$X'_k X_k = \mathcal{G}_{k}$
We already have $\mathcal{G}_{k-1}$, and $\mathcal{G}_{k} = \mathcal{G}_{k-1}$ except for a new row and column of dot products, so let's say $$\mathcal{G}_{k} = \left[ \begin{array}{cc}
\mathcal{G}_{k-1} & \mathbf{p} \\
\mathbf{p}' & x \end{array} \right]$$
We can find the Cholesky decomposition $L_k$ from $L_{k-1}$:
$$L_{k} = \left[ \begin{array}{cc}
L_{k-1} & \mathbf{r}=L_{k-1}^{-1}\mathbf{p} \\
\mathbf{r}' & y=(x - \mathbf{r}\cdot\mathbf{r})^{-\frac12} \end{array} \right]$$
Find $L_k^{-1}$. Note this explicit inversion.
$\mathcal{G}^{-1}_k = {L'}_k^{-1} L_k^{-1}$
$A_k = (1'_k \mathcal{G}^{-1}_k 1_k)^{-\frac12}$, where $1_k$ is a vector of $k$ 1s.
${\mathbf w}_k = A_k {\mathcal G}_k^{-1} 1_k$
Use ${\mathbf w}_k$ to select the $(k+1)$th predictor.

Now here's my approach:

At $k$th iteration:
Find $L_k$, $\mathbf{r}$, and $y$ as above
Define $L_k {\mathbf q}_k = 1_k$
Find $\mathbf{q}_k$:
$${\mathbf q}_{k} = \left[ \begin{array}{c}
{\mathbf q}_{k-1} \\
\text{$\scriptsize{(1 - {\mathbf q}_{k-1} \cdot \mathbf{r})/y}$} \end{array} \right]$$
$\begin{align}A_k &= (1'_k \mathcal{G}^{-1}_k 1_k)^{-\frac12} \\
                     &= (1'_k L'{_k}^{-1} L_k^{-1} 1_k)^{-\frac12} \\
                     &= \left[(L_k^{-1} 1_k)' (L_k^{-1} 1_k)\right]^{-\frac12} \\
                     &= ({\mathbf q}'_k {\mathbf q}_k)^{-\frac12} \end{align}$
$\begin{align}{\mathbf w}_k &= A_k {\mathcal G}_k^{-1} 1_k \\
                               &= A_k L'{_k}^{-1} {\mathbf q}_k \end{align}$
So solve for ${\mathbf w}_k$: $\frac{1}{A_k} L'_k {\mathbf w}_k = {\mathbf q}_k$

Comments? Does it even matter if I'm right?

If you want to solve for single value of $ \lambda $ in the model:

$$ \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} $$

Then you can use Coordinate Descent method which is the fastest and simplest and doesn't require any matrix inversion.

I have a MATLAB code for in my $ {L}_{1} $ Regularized Least Squares Analysis GitHub Repository:


How to find the convolution kernel in frequency domain? I have a two vectors of spatial data (each about 2000 elements in length). One is a convolved version of the other. I am trying to determine the kernel that would produce such a convolution. I know that I can do this by finding the inverse Fourier transform of the ratio of the Fourier transforms of the output and input vectors. Indeed, when I do this I get more or less the shape I was expecting. However, my kernel vector has the same dimensionality as the two input vectors when in reality the convolution was only using about one fifth (~300-400) of the points. The fact that I am getting the right shape but the wrong number of points makes me think that I am not using the ifft and fft functions quite correctly. It seems like if I were really doing the right thing this should happen naturally. At the moment I am simply doing;

FTInput = fft(in);
FtOutput = fft(out);
kernel = ifft(FtOutput./FTInput).

Is this correct and it's up to me to interpret the output vector correctly or have I oversimplified the task? I'm sure it's the latter, I'm just not sure where.

If you have noise present in your signal, the straightforward Fourier domain division will cause plenty of errors in your result. Some ways to avoid it are by using the so-called dual channel FFT (Part 1 and Part 2). I can also suggest deconvolution via adaptive filters, LMS or NLMS ([Normalized] Least Mean Squares) filters in particular are easy to understand are not hugely expensive in terms of CPU cycles in case your signals are long. LMS adaptive filters are very robust to noise.


How to normalize output of $\tt scipy.signal.correlate$ I have 2 different signals and I'm trying to cross-correlate then using Python 2.7 and scipy.signal.correlate. How do I normalize my results (such that the max amplitude is 1.0? I tried the following:

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(0,2*np.pi,num=1000)
x1 = 10*np.sin(2*np.pi*t)
x2 = np.sin(2*np.pi*t+np.pi/2)
x12 = scipy.signal.correlate(x1,x2,'full')
plt.plot(x12)

results in the following plot

When you say normalized cross-correlation I guess you mean the Pearson correlation.
Anyways you just divide the cross correlation by the multiplication of the std(standard deviation) of both signal, or more conveniently:
$ \rho_{xy}  =\frac{<x,y>}{\sigma_x\sigma_y}$

and in code:

x1 = x1/x1.std()
x2 = x2/x2.std() and then as you did it


Is it possible to replace an integrator system with an equivalent differentiator? I have a system whose input-output relation is as follows

$$y(t)=x(t)+\int_{-\infty }^{t} x(\tau) \,\mathrm d \tau$$

Can I create an equivalent system by using differentiators rather than integrators?

I think something like taking derivative of both sides of the equation but improper integral makes it hard. Is there any nice way to convert this system? Thanks in advance.

Typically you can devise equivalent block diagrams in terms of just integrators or just differentiators.

A simple example is on page 2-15 of

https://web.stanford.edu/~boyd/ee102/systems.pdf

which should look remarkably familiar. You should be cautioned that mathematical equivalence doesn't mean that there aren't reasons where one form would be prefered over the other.  OP Amp integrators have better noise characteristics. You also have to band limit an analog differentiator to make it physically realizable.

Another issue is that they will not in general have the same complexity.

When modeling ordinary differential equations systems,
$$
\mathbf{M}(t) \frac{d \mathbf{y}(t)}{dt}= \mathbf{g}(t,\mathbf{y})
$$
where $\mathbf{M}(t)$ (mass matrix) is invertible is often more sparse for physical systems than
$$
\frac{d \mathbf{y}(t)}{dt}= \mathbf{h}(t,\mathbf{y})
$$
which is a more general way of writing the common state space representation.

One can write :
$$
\frac{d \mathbf{y}(t)}{dt}= \mathbf{M}^{-1}(t)\mathbf{g}(t,\mathbf{y})
$$
but if $\mathbf{M}(t)$ is sparse, $\mathbf{M}^{-1}(t)$ probably isn't ,which makes a big difference when the dimension of $\mathbf{y}$ is large.

The Matlab Control systems toolbox also includes state variable equations of the form

$$
\mathbf{E} \frac{d \mathbf{y}}{dt}= \mathbf{A}\mathbf{y}+ \mathbf{B}\mathbf{u}$$


inverse fourier transform of magnitude and phase I stuck this question. Frequency response is written as magnitude and phase and I don't find inverse fourier given signal which as magnitude and phase.How can I solve it ? Can you explain the solution way in order? magnitude of signal is : y=−ω while −1≤ω<0,
y=ω while 0≤ω<1 phase of signal is: y=−3ω

What is the Fourier transform? I mean I have phase and magnitude information of a signal that has been transformed from fourier. I want to find inverse fourier transform. And phase and magnitude information are above

You can express a complex signal $y$ as $|y|.e^{i \phi_y}$ where $|y|$ is the magnitude and $\phi_y$ is the phase. You can use Inverse the Fourier transform formula:
$$\displaystyle y(t) = \frac{1}{2\pi}\int_{-\infty}^\infty Y(\omega) e^{j\omega t} d\omega$$
You can express $Y(\omega)$ in $|y|.e^{i \phi_y}$ form and then perform the IFFT. But I am not sure if a minus comes in magnitude unless its in dB.


What does the amplitude in a frequency response of a signal signify? Well according to my understanding it tells the concentration of signals at a particular frequency but what does the number really mean? Suppose in a frequency response graph i have such a relation ::
"a1"(amplitude) at "f1" (frequency)
a2 at f2 , a3 at f3. What do i infer from it?

Short answer, it depends what exactly you're looking at.

If you have a time series and compute the DFT, this will give you the concentration of energy at each discrete frequency point you have chosen. When you compute the PSD of a signal in MATLAB (using pwelch or cpsd), you'll get units of power/frequency. This describes the amount of energy or power in that specific frequency range (bin). If instead your goal is to relate the amplitudes of sinusoids in a time series directly to the values that you're reading off the peaks in the PSD plot, you need to ensure that your FFT is appropriately scaled. In the MATLAB functions cpsd() or pwelch(), the FFT is scaled such that the $\textit{total energy}$ in the signal is preserved (done by dividing by the FFT by the sampling period, $t_s$). If you want to recover the original amplitudes of a signal (e.g., $x=A_1 sin(\omega_1*pi*t)+A_2 sin(\omega_2*pi*t))$, you would want to ensure that the FFT is scaled by the duration of the signal, $T$.

Basically, to make any accurate inference from the frequency response plot, you need to know exactly how it was computed (units).

Some refs. that I found helpful:
What is the difference between the PSD and the Power Spectrum?


Apply windowing before FFT or after? I am computing a spectrogram and I've found a code example online which goes like this:

for jj = 1:size(signal_framewise,2)
    current_frame = signal_framewise(:,jj);
    dtf = fft(current_frame).*gausswin(window_length_s);
    %nfft is half of the fft results (since the fft is symetric)
    out_buffer(:,jj) = dtf(1:nfft);
end

But I intuitively did it like this:

for jj = 1:size(signal_framewise,2)
    current_frame = signal_framewise(:,jj).*gausswin(window_length_s);
    dtf = fft(current_frame);
    out_buffer(:,jj) = dtf(1:nfft);
end

The difference is, that the guy online applies the window to the finished FFT and I apply it to the signal frame. Question is: What's the right way to do and does it make a difference? I've attached the outputs of the spectrogram and the full code.

[![%% my implementation][1]][1]

clc; clear all; close all;
Fs = 44100;
t_max = 3;
T = 1/Fs;

time = 0:T:(t_max-T);
input = chirp(time,1500,1,8000);

window_length_t = 0.01; %10ms window length
window_length_s = round(0.01 * Fs); %window length in samples
if mod(window_length_s,2) == 0
    window_length_s = window_length_s + 1; %make sure we have odd window size
end

%activate the following for a no overlap implementation
%signal_framewise = buffer(input , window_length_s);
signal_framewise = buffer(input , window_length_s , floor(window_length_s/2));
nfft =((window_length_s-1)/2)+1;

out_buffer = zeros(nfft,size(signal_framewise,2));
for jj = 1:size(signal_framewise,2)
    current_frame = signal_framewise(:,jj).*gausswin(window_length_s);
    dtf = fft(current_frame);
    out_buffer(:,jj) = dtf(1:nfft);
end

F = linspace(0,Fs/2,nfft);
T = linspace(0,3,size(signal_framewise,2));
surf(T,F,20*log10(abs(out_buffer)), 'EdgeColor', 'none');
axis xy;
axis tight; colormap(jet); view(0, 90);
xlabel('Time');
colorbar;
ylabel('Frequency(Hz)');
set(gca, 'YTickLabel', num2cell(get(gca, 'YTick')));

Element-wise multiplication in the spatial/temporal domain with a windowing function serves to reduce the effect of the potentially large jump you get at the signal edge, when making the sampled signal periodic. This jump would otherwise introduce lots of frequencies that might not be present in the signal.

Element-wise multiplication in the frequency domain with a windowing function applies a low-pass filter. Multiplying by a Gaussian window is equivalent to applying a Gaussian filter in the spatial/temporal domain.

If the code you found online were to compute the inverse transform after applying the window, then whoever wrote it would be low-pass filtering. With the code as you posted, I tend to think that person doesn't know what they're doing (or maybe it's just a typo).


Calculating incoherency in compressive sensing CS is based on a choice of a sensing basis $\Phi$ relative to a representation basis $\Psi$. There are many well known pair matrices for $\Phi$ and $\Psi$ like random Gaussian and FFT, and also they should have low coherency which can be calculated with:

$$\mu(\Phi, \Psi) = \sqrt n\cdot \max\limits_{1\leq k, j \leq n} \left|\langle \varphi_k, \psi_j \rangle\right|$$

Some of these Coherency values which are calculated are as: between noiselets and Haar wavelets is $\sqrt 2$. Coherence between noiselets and Daubechies $\rm D4$ and $\rm D8$ are $2.2$ and $2.9$, respectively; random matrices are largely incoherent with any fixed basis $\Psi$.
My question is how can I calculate these values in matlab? can anyone provide a code example for me?

As a starting-point, you might try something like this:

Step 1: create random (complex-valued) matrices

rows = 32;   %number of output samples from compression matrix
cols = 2048;  %number of input samples supplied to compression matrix
matrixA = randn(rows, cols) + 1i*randn(rows, cols);
matrixB = randn(rows, cols)+ 1i*randn(rows, cols);

Step 2: Compute Maximum Dot-Product Magnitude

%compute all possible cross-correlations
xc = abs(matrixA * matrixB');
%find the maximum
result = max(max(abs(matrixA * matrixB')));
%todo:  scale by sqrt(n)
%...

Step 3: Beware of Amplitude Scaling Impacts

Since the coherence equation that you provided depends upon dot-products, the amplitude scaling of your matrix will effect the result.  Most likely the equation you provided assumes the matrices are scaled in amplitude.  For example, each column in matrixA & matrixB might need to be scaled by it's "Euclidean Norm" or "L2-Norm".  You might find it helpful to read about correlation coefficients which are scaled accordingly, or revisit your equation reference just to be safe.


Lowpass Butterworth filter design (own poles' calculations) I'm trying to calculate coefficients for my Butterworth filter of my own. I've found this article that looks really reliable. But I'm trying to implement it in MATLAB and I'm getting bandstop filter.
Please find attached matlab script below:

clc, clear all, close all;
%DIGITAL BUTTERWORTH FILTER DESING
%Butterworth LP with -3dB at fc and N order and G gain
fs = 10; %frequency sampling
fc = 1;   %cut off frequency
N = 2;      %filter order
G = 1;      %filter gain
wc = 2*pi*fc;

%analog filter design
%poles and zeros
k = [1:1:N];
%poles searching
%sk = wc * exp(((j*pi)/(2*N))*(2.*k+N-1))
sk = wc * exp(j*(pi/2)) * exp(j*(2.*k+1)*pi/2*N);
z = 0;
G = 1;
[NUM,DEN] = zp2tf(z,sk,G);
%no zeros for LP
Hs = tf(1,DEN)

%filter analysis
figure(1);
subplot(3,1,1)
pzmap(Hs);
subplot(3,1,2)
impulseplot(Hs);
subplot(3,1,3)
stepplot(Hs);
fvtool(1,DEN);

I've calculated poles according to equation $10.3.47$ in the article and created transfer function with this poles as a parameters. But fvtool shows bandstop filter instead of lowpass (as described in article). In pzmap I'm getting also pole on right side of $s$-plane so it seems to be unstable. I dont really know how to interpretet this results.

Maks

Ok, I've solved this problem myself. I've used equation from this article mentioned in first post. I've just implement it in matlab not enough precise.

clc, clear all, close all;

%ANALOG BUTTERWORTH FILTER DESING
%https://engineering.purdue.edu/~ee538/AnalogFilterDesigns.pdf
%Butterworth LP with -3dB at fc and N order and G gain

wc = 1000;   %cut off frequency
N = 2;      %filter order
G = 1;      %filter gain

%analog filter design
%poles and zeros

%poles searching
for k = 0:N-1
   sk(k+1) = wc*exp(j*(pi/2))*exp(j*(2*k+1)*(pi/(2*N)));
end

[NUM,DEN] = zp2tf([],sk,G);
b0 = DEN(end)*G      %gain normalization

Hs = tf(b0,DEN)

figure();
freqs(b0,DEN)
figure();
plot(sk./10^2, 'x')
xlim([-10,10]);
title('Poles placement in s-plane')

Many thanks for your help,
Maks


What exactly does ordfilt2 do? Can you please explain in simple terms what do the input parameters indicate in the ordfilt2 function in matlab?

B=ordfilt2(A,Order,Domain)

I have seen people use this function as J = ordfilt2(I, 9, true(5)) or J=ordfilt2(I,25,ones(3,3))

But I do not understand what each input does to the image I to give J as output..

Thanks a lot in advance

I'm pretty new to this myself, so please correct me if I get this wrong.

Using your example, J = ordfilt2(I, 9, true(5)).
ordfilt2 will move over the 2d array I in blocks of the same size as true(5). For each of these 5x5 blocks, sort all the elements from smallest to largest. Now fill in the corresponding block in J with a bunch of copies of the 9th smallest element.

To use a smaller example, so I have room to type it:

I = [ 1 2 4 5 ;
      5 3 5 1 ;
      0 3 5 2 ;
      2 1 7 7 ];
J = ordfilt2(I,3, ones(2,2));

Now, let's go through a few blocks one at a time. The first ones(2,2) block is [1 2; 5 3] in the top left corner. If we sort these elements, we get [1 2 3 5], and since we're looking for the 3rd smallest, we receive a value of 3 for the (1,1) position of J.

J = [ 3 ? ? ? ;
      ? ? ? ? ;
      ? ? ? ? ;
      ? ? ? ? ];

next up is the I(1:2,2:3) block. Ordering those elements gives [2 3 4 5], so the third smallest is 4. Now we replace J(1,2) with a 4.

J = [ 3 4 ? ? ;
      ? ? ? ? ;
      ? ? ? ? ;
      ? ? ? ? ];

Go ahead and run this command and compare I, J to understand what's happening. One other thing to know is that the input matrix is padded by default with zeros at the lower and right sides. So the I(4:5,4:5) block is [ 7 0; 0 0].


A few questions about the output value of dsp.SpectrumAnalyzer I am testing the output values of dsp.SpectrumAnalyzer for a sine wave mixed with white noise. The $\sigma^2$ of the white noise is 1e-4, so I use 0.01*randn(1024,1) to generate it. The sine wave is generated by dsp.SineWave with default amplitude (=1). The sampling rate is 44100 Hz. If the 'SpectrumType' property of dsp.SpectrumAnalyzer is set by default ('Power'), the scope output is as follows:

According to the analysis in "Convert the Power in Watts to dBW and dBm" section of Estimate the Power Spectral Density in MATLAB, The peak is at 26.9897 dBm, which is what the above screenshot shows us. But other values in the figure confused me. I will ask them one by one.

(1) The power spectrum of the white noise is $\sigma^2$=1e-4 (right?), which is 10*log10(1e-4/1e-3)=-10dBm according to the definition in the above link. But the figure shows that the power spectrum of the white noise is at -40dBm. How is this value obtained? PS: if we calculate the power spectrum in dBW, 10*log10(1e-4)=-40 which is close to what the figure shows, but here I am using dBm unit!

(2) If I change the 'SpectrumType' property of dsp.SpectrumAnalyzer to 'Power density', the output of the System object is as follows:

According to the explanation of this property in documentation dsp.SpectrumAnalyzer System object, now the peak of the sine wave should be "the magnitude squared of the spectrum normalized to a bandwidth of 1 hertz.". The spectrum is 0.5, so its square is 0.25, what I expected to get in dBm should be 10*log10(0.25/1e-3)=23.9794dBm, but what I can see in the figure is 12.6948 dBm. How is this value obtained?

(3) Let's look at the power spectrum density of the white noise again. According to the definition mentioned in the above question, the value should be 10*log10(1e-4^2/1e-3)=-50dBm. I'm not sure if the calculation is correct, but it is close to the figure, although there seems to be some distance. The question is, the section "Convert the Power in Watts to dBW and dBm" of Estimate the Power Spectral Density in MATLAB says the $P_{whitenoise}$=1e-4*(22050/21.53)=0.1024, but if we calculate the dBm, we get 10*log10(0.1024^2/1e-3)=10.2060. To my surprise, the webpage then says "the total power of the white noise is -39.87 dBm." How is this -39.87 dBm obtained?

I'm totally confused about these values. I will greatly appreciate it if you can help me clarify them. Thanks a lot!

The following are some reference code for your experiment. I am using the latest Matlab R2016a.

Fs = 44100;
Sineobject1 = dsp.SineWave('SamplesPerFrame',1024, 'PhaseOffset',10,...
    'SampleRate',Fs,'Frequency',5000);
SA = dsp.SpectrumAnalyzer('SampleRate',Fs,'SpectrumType','Power',...
    'PlotAsTwoSidedSpectrum',false,'Window','Hann',...
    'ChannelNames',{'PSD of the input'},'ShowLegend',true);
for Iter = 1:5
    Sinewave1 = step(Sineobject1);
    Input = Sinewave1;
    NoisyInput = Input + 0.01*randn(1024,1);
    step(SA,NoisyInput);
end

I don't really know what Matlab's spectrumanalyzer does internally, but I suspect it mimics the behaviour of an actual spectrum analyzer – which does not at all do a DFT and then square the output.

Instead, a SA shifts a filter across the spectrum. This filter's equivalent noise bandwidth is exactly how much white noise power gets let through and gets measured to derive a single point in that display. So, for white noise, there's always the factor "filter bandwidth" between your power spectral density and the power displayed in your plot.

Also, filters don't have to be unitary, i.e. there might be yet another factor, the filter gain, between your sine's power and the displayed power.

Also, in a SA you typically need to configure the "frequency advance step" separately from the filter bandwidth, hereby defining how much the windows "overlap". A quick skim of the matlab doc says the default overlap is 0 – so you might really just be missing your sine wave's "peak".

All in all, it's quite questionable why you'd want to use this visualization if you could also directly estimate the spectrum all at once, without losing any info, using any of the PSD estimators that Matlab can do – from the simple (windowed) FFT estimator, to the more complex parametric estimators.

This tool definitely only has the purpose of illustrating how measurements on a Scope can be modeled and how they differ from the mathematical idea of the PSD. Unless you need to model a SA, I really don't see any practical value in this tool.


How many times does a certain sound appear in an audio clip? I have an audio clip which can be found here  or here. If you hear the audio clip you will notice that there are six times a kind of hit sound in the clip.

My questions are:

How can I count these hit sounds (I should get six in this clip)?
What is the best tool/software/app for this problem?
Where can I find a good source that explains the theory behind this?

How can I count these hit sounds ( I should get six in this clip).

In this case you can simply rectify the signal (i.e. run it through an absolute function, then through a simple running average low pass filter and at the output of this add a threshold. Essentially, when the power of the input signal exceeds the threshold you will be receiving a pulse whose width is approximately equal to the duration of the "hit". Alternatively, you can do away with rectification and filtering and simply pass the whole thing via a threshold operator. That is, when the power of the input signal exceeds the threshold, you get a "pulse". In either way, count the "pulses" and you can count the number of "hits" that you have in that wave file.

  What is the best tool/software/app for this problem.

This is not an incredibly challenging Digital Signal Processing task, so you could simply write a Python script that does it in 3-5 lines of code (I am not implying that Python is not capable of serious processing here). Alternatively, you could install something like scilab and again do it in 3-5 lines of code, using one of its toolboxes (or, this one).

  Where can I find a good source that explains the theory behind this.

This "rectify and threshold" technique is a sort of a "standard" in counting pulses, alternatively, you might want to have a look at the moving average filter as an example of a sliding window technique where you process a long signal in small "frames" of samples. Other than this, this website has an extensive list of book references if you search for "books".

Hope this helps.


What is the algorithm to do frequency modulation? I'm familiar with audio synthesis, oscillators etc. I'm now trying to implement my own software algorithm to do FM as described by https://en.wikipedia.org/wiki/Frequency_modulation.

What I have implemented outputs the following (spectrogram produced with sox):

This was produced by my algorithm set with a property of:
A 1200hz sine wave modulated by a 12hz sine wave.

However, I don't think this is actually 'correct' in what FM should do. I'm assuming I would see a 'signal' increasing and decreasing in frequency.

Is this correct? Or am I doing something wrong?

The way I am calculating the amplitude for each time point (code is written in Haskell):

fmSine :: Frequency -> Signal' -> Signal'
fmSine (Frequency hz) (Signal sng) = Signal (\t -> sin ( (sng t) + (hz* 2 * pi * t)))

Essentially I'm doing sin ((amplitudeOfModulationFrequency t) + hz * 2pi * t)
Where hz is the frequency and amplitudeOfModulationFrequency t is the amplitude at time t of the modulating frequency.

What your plot shows is noise from varying size discontinuities in the waveform.

When you change a modulating frequency, you have to make sure the resulting FM waveform does not have discontinuities.  You have to do this by incrementally changing the phase of the argument to the sin() function, not by naively multiplying the modulating frequency by time, which will jump the phase by some (psuedo)random angle when the frequency changes.