What's the difference between spreading and scrambling I have a question regarding the scrambling code, What's the difference between spreading and scrambling? As I know, following the below figure, the scrambling doesn't affect the signal bandwidth:

As shown above, the chip rate after the scrambling code is still the same,I mean it's not spread. My question, does scrambling is performed only by multiplying the sequence resulted from spreading, element by element, with scrambling code? or how is it processed ?

Thank you

Spreading refers to the process of multiplying your data sequence with a higher rate PRN sequence often but not necessary such that one complete duration of a repeating sequence occupies the time duration of one data symbol. For example one such modulation using BPSK to send a “1” you send the entire sequence before it repeats and to send a “0” you send the inverted sequence. Thus the resulting spectrum is “spread” since the resulting BPSK pattern is at a much higher rate.

Scrambling in contrast is multiplying your data with a rate matched PRN sequence. This does not increase the pattern rate but does “whiten” your data which is very useful in many applications such as avoiding long runs of “0” or “1”, and can also help provide security— although linear codes should be avoided for such purposes as they are not very secure.

So to process something as in your figure, the sequence would need to be descrambled by multiplying it with the scrambling code and then it can be despread by multiplying it with the spreading code. This can be done in one step (although less likely) by multiplying with a scrambled spreading code.

The diagram as shown does add a layer of security as long as the scrambling code is not repeating synchronous with the spreading code and the adversary has no access to plaintext.


Square Wave with Random Edge Shifting I took a square wave signal and individually shifted all the rising and falling edges by a random amount of time, within certain +/- bounds (a random time for each edge, not shifting all edges by the same time). This is what the FFT looked like before the edge shifting:

After the random edge shifting, the FFT looks like this:

I did not expect to still see sharp peaks at the signal frequency and odd harmonics. Rather, I expected that the peaks would become wider, with a width corresponding to the maximum frequency deviation that the time shift bounds allowed. Can someone explain what is going on here?

Edit: How is it different than this?

By separately changing both the leading edge and the trailing edge of each pulse, you have converted the square wave into what can be viewed as an asymmetric double-edge pulse-width-modulated (DEPWM) signal. Analysis of the Fourier spectra of various types of deterministic PWM signals (including the type mentioned in the previous sentence) can be found in the paper Z. Song and D. V. Sarwate, "The Frequency Spectrum of Pulse Width Modulated Signals", Signal Processing, vol. 83, pp. 2227-2258, 2003. (The article is available behind the publisher's paywall but there are some free copies floating around the web).  The results show that the spectra of PWM signals do have the spikes shown in the OP's diagram just as the square wave spectra do, but there is additional clutter between the peaks. That is, the PWM Fourier spectrum is not a Fourier series like the square wave spectrum is; PWM spectra are of a continuous-time signal that contains periodic components (giving those spikes) as well as non-periodic components giving the clutter.  The details (and the nuances) are too messy to transcribe here; interested readers should read the original.


Would the ratio of Autocorrelation value and autocorrelation of noise for different data types be different? in Matlab Considering a model of the form:

$y[n] = x[n] + noise[n]$ and

$y[n] = u[n] + noise[n]$

where the term noise represents Additive White Gaussian Noise of variance $\sigma^2_{noise}$.
I want to understand if there is any difference in the value of the terms: $\sigma^2_x/\sigma^2_{noise}$ and $\sigma^2_u/\sigma^2_{noise}$ by varying the noise level.

$x$ represents a signal generated from Gaussian distribution whose variance is $\sigma^2_x$. In Matlab, I have used the command rand() to generate $x$.

$u$ represents another signal taking values zero or 1 having variance $\sigma^2_u$. In Matlab, I have used rand()>0.5 to generate thie binary valued data.

Based on my understanding the numerator in the ratio is the autocorrelation of the data mathematically expressed as $E[\mathbf{x}[n]\mathbf{x}[n]^T]$ and  $E[\mathbf{u}[n]\mathbf{u}[n]^T]$. Why are these terms giving same result?

In most papers, the variance of the data irrespective of whether it takes numeric or symbols is considered to be 1.

Question1) Is my approach of implementing the ratio correct? I cannot understand what the value of $\sigma^2_x$ and $\sigma^2_u$ should be in practice. Implementing this way would give the same value for the ratio for both the different kinds of data $x$ and $u$. Would the auto correlation for numeric values x and binary data u be the same value?

Question2) If the ratio is the same then what is the need for using binary signal?

N = 128;
var_x =1;
var_u = 1;
signal_x = rand(1,N);
signal_u = rand(1,N)>=0.5; % this creates 0/1 data
noise = 0;
index = 1;
for noise = 0:5:30
noise_var  = 10^(-noise/10);
ratio_x(index)= 10*log10(var_x/noise_var); %in dB
ratio_u(index)= 10*log10(var_u/noise_var);
index = index+1;
end

I am not really sure what you want to do so please comment and I will modify my answer.

The case that $x$ and $u$ seems statistically identical is that $x$ is standard Gaussian random variable and $u$ has $\pm 1$ values.

Use your notations, $x$ is Gaussian random variable with zero mean and unit variance (deducted from the fact that you used the standard rand function, even though the randn() function must be used for Gaussian pdf).

Let's stay with real $x$ for the sake of simplicity :

$$p(x) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$$

with $\mathbb{E}[x] = \mu = 0 \textrm{ and } \mathrm{var}(x) = \sigma^2 = \mathbb{E}[x^2] = 1$.

Thus the probabilites $\mathrm{Pr}(x >= 0) = \mathrm{Pr}(x < 0) = 0.5$.

The processed data $u = 1 \textrm{ if } x \ge 0 \textrm{ and } u=-1 \textrm{ otherwise}$. Thus $\mathbb{E}[u] = (1) \times \mathrm{Pr}(x >= 0) + (-1) \times \mathrm{Pr}(x < 0) = 0.$

Then $\mathrm{var}(u) = \mathbb{E}[(u - \mathbb{E}[u])^2] = \mathbb{E}[u^2] = (1)^2 \times \mathrm{Pr}(x >= 0) + (-1)^2 \times \mathrm{Pr}(x < 0) = 1$.

The means and variances are just coincidentally identical.

MATLAB/Octave code

octave:6> N = 65536;
octave:7> x=randn(1,N);
octave:8> u=2*(randn(1,N)>0) - 1;
octave:9> var(x)
ans =  1.0020
octave:10> var(u)
ans =  1.0000

The case of OP question $x$ is (0,1)-uniform (by rand()) and $u$ is binary:

$\mathrm{Pr}(x >= 0.5) = \mathrm{Pr}(x < 0.5) = 0.5$.

$\mathbb{E}[x] = 0.5$

$\mathrm{var}(x) = \int_0^1 (x-0.5)^2 \mathrm{d}x = \frac{2}{3}\times 0.5^3 \approx 0.8$

The processed data $u = 1 \textrm{ if } x \ge 0.5 \textrm{ and } u=0 \textrm{ otherwise}$. Thus $\mathbb{E}[u] = (1) \times \mathrm{Pr}(x >= 0.5) + 0 \times \mathrm{Pr}(x < 0.5) = 0.5$.

$\mathrm{var}(u) = \mathbb{E}[(u - 0.5)^2] = (1 - 0.5)^2 \times \mathrm{Pr}(x >= 0.5) + (0 - 0.5)^2 \times \mathrm{Pr}(x < 0.5) = 0.25$

MATLAB/Octave code

octave:13> N = 65536;
octave:14> x=rand(1,N);
octave:15> u=rand(1,N)>0.5;
octave:16> var(x)
ans =  0.083424
octave:17> var(u)
ans =  0.25000

The variances are different !!!


Correct way to add AWGN to a signal I have a signal, S whose bandwidth is bw Hertz and is sampled at fs Hertz.

To contaminate it with a noise corresponding to SNR x dB, I used the matlab function

out = awgn(S,x);

But I have came across some codes doing it instead like the following.

out = awgn(S,corrected_x);

Where corrected_x = x + 10*log(bw/fs);

Which is the correct method and why?

Short answer

10*log(bw/fs) to take into account the oversampling operation because the awgn() function specifies the signal-to-noise ratio per sample, in dB.

Longer answer

The discrete time AWGN model is

$$Y = X+N$$

where X is data from continuous time $X(t)$, N is noise sequence from AWGN process $N(t)$ and Y is receive symbols.

If $X(t)$ is characterized by its baseband equivalent limited between $[-W/2, +W/2] \textrm{ (Hz)}$,  then we can identify $X(t)$ by observing $Y$ at a rate $W$ symbols per second. See chapter 2, sampling theorem and Theorem of irrelevance.

Call $P$ the average power (joules per second). The sample power $E_s= P/W$ and the noise symbol power is $N_0$. The signal noise ratio per symbol is defined $\mathrm{SNR} = \frac{P}{BW_0} = E_s/N_0$.

If the complex baseband signal is oversampled $\alpha = f_s/W$, the noise sample power is still $N_0$ while the data sample power is reduced $\alpha$ times, thus $E_s/N_0 = \mathrm{SNR} \times W/f_s$.

In dB, $E_s/N_0 = \mathrm{SNR} + 10\log_{10}{(W/f_s)}$.

The awgn() function add AWG noise by the previously defined $E_s/N_0$.


Finding the minimum phase h[n] and its Z transform Hello, this is one of my homework questions and i have already solved the first question but im having trouble gettin a relation that helps me solve the second one. From the question i understand that R(z)=H(z).H(1/z) but i dont know how to develop this.
This is my response to the 1st question

Thanks

Since this looks like a homework problem I'll help you to solve it on your own, rather than just present you the solution.

The problem here is that you have a mistake in the first part of the solution. The evaluation of the first sum (from $n=-\infty$ to $n=-1$) is wrong. Redo the calculation and you should end up with

$$R(z)=\frac{1}{(1-az^{-1})(1-az)}\tag{1}$$

with some appropriate value for $a$. From $(1)$ it's immediately clear how you can factor $R(z)$ as

$$R(z)=H(z)H(z^{-1})\tag{2}$$

and how you have to choose $H(z)$ to get a minimum-phase systems.


why do low pass and high pass filters generate a phase shift? I just learned following notation for a sine wave:

$A e^{j\phi t}$

A low passfilter and a highpass filter respectively generate a phase shift in the complex plane of the sine wave as follows:

$A e^{-j\phi t}$ and $A e^{+j\phi t}$

I don't understand what phase shifting has to do with letting higher or lower frequencies intact... How do these things( filtering and phase shifting) relate to each other? They just seem like to completely and independent things to me

Because they need a certain time delay in order to "decide" what the output should look like if the output changes a lot. Imagine you see a number series and your brain has to decide what the mean value of these number series is. You will need to look at several numbers to arrive at some reasonable conclusion (let's assume the numbers are represented by a noise distribution with the same mean). The same applies to filters that weigh signal values to arrive at the filtered output signal.

Any filter that does not employ prediction of output will have phase shift. The only way to get rid of phase shift completely (or keep it to fractions of dB) is to use a predictive filter - for example kalman filter. Such a filter, in the simplest form will estimate the derivative of the signal and use that to filter the original signal in a predictive fashion. If you have a model of your signal you can do even better. Since prediction and sensor measurements can be processed in the same time step, from a discrete time point of view you can have zero phase shift if prediction is able to estimate the original signal precisely.


Difference between Fine and Coarse Frequency Offset in OFDM I need to explain for a project the difference between Fine and Coarse Frequency Offset in OFDM signals. I know what the frequency offset is and what it does to the carrier frequencies but I cannot understand where the terms "coarse" and "fine" refer to. Those terms are also used in GNU Radio in the Schmidl and Cox sync block.

Most of Google results take the fine and coarse frequency offset as granted and there is no specific explanation.

Can someone explain this to me?

If in doubt, read the original Schmidl&Cox Paper Robust Frequency and Timing
Synchronization for OFDM. It's actually a good read, and a good paper with respect to demonstration of state of the art, proposed algorithm, theory, and simulation. I highly recommend it.
Anyway, to understand this, you should consider OFDM as a multicarrier system. In fact, the DFT used in OFDM can be understood as Sinc filter bank:

The spacing of the sincs is actually the OFDM bandwidth divided by the size of the DFT. Notice that the important property of the sincs is that they have a zero at every integer multiple of the subcarrier distance – so that one carrier doesn't interfere with the others.
Now, assume your signal only has a small offset, let's say ca ¼ of the subcarrier spacing:

Then, Schmidl and Cox only has to find the offset so that the zeros of the received sincs align with the neighboring peaks of the "ideal" sinc-filterbank-shape.
Now, that's the fine frequency correction.
The coarse frequency correction becomes important when the frequency offset is larger than one subcarrier distance; let's look at something that is shifted by ca 3¼ of a subcarrier distance:

Then, the fine correction first aligns the receive signal so that the zeros/peaks match, but there's a remainder of full subcarrier multiple offset:

The job of the coarse frequency correction is to eliminate the remaining 3 subcarrier bandwidths of offset.


How would one design a (quasi) linear phase adaptive notch filter for a single complex tone? While IIR notch filters are attractive, I need to retain phase linearity at the filter output.  I imagine that it's possible to use a standard IIR notch  filter:
https://www.researchgate.net/publication/224174145_A_Complex_Adaptive_Notch_Filter
followed by some form of phase equalization.
An alternative would be an FIR filter with conjugate symmetric coefficients.  I am unaware of any FIR designs in the literature explicitly for a single complex tone.

Any thoughts on the best way forward?

I might consider doing a complex modulation of the frequency you want to notch down (or up) to 0 Hz.  Then run your favorite deep DC blocker (plus phase adjustment if needed).  Then complex remodulate the result back to the original unshifted spectrum.

DC block: For non real-time, just subtract the entire signal’s average.  Or a long (weighted?) moving average if you can tolerate the latency (Rick Lyons has some notes on this technique).


How to find the coefficients of the following differential equation An arbitrary signal $v(t)$ pass through the following system,

$w'(t) + 5 w(t) = v'''(t) + 320v''(t) + 40 v' (t) + 40v(t)$

How to determine the coefficients of the following differential equation, where the input signal is $w(t)$ and the output signal is $z(t)$,
So that the output $z(t) = v(t)$ (after long time when there is  no transient exist any longer)?

$z'''(t) + d_2 z''(t) + d_1 z'(t) + d_0 z(t) = e_3 w'''(t) + e_2 w''(t) + e_1 w'(t) + e_0 w(t)$

Since this a homework type problem, only a few hints that should help you solve this problem yourself:

determine the transfer function of the first system by transforming the given differential equation and calculating $H_1(s)=W(s)/V(s)$.
do the same with the second system to obtain $H_2(s)=Z(s)/W(s)$.
from $Z(s)=V(s)$, determine $H_2(s)$ in terms of $H_1(s)$.


How to remove the modulation before doing frequency offset estimation? To use DFT/FFT (or maximum likelihood) method to estimate the frequency offset introduced by the channel, we need to remove the modulation on the received data samples in the front. If the unknown data samples belong to a QPSK constellation, how should this modulation removal be done?

For lower-cardinality PSKs like the QPSK, the "classical" way is to
take the signal to the $M$th power, $M$ being the number of constellation points.
From the shape of an $M$-PSK modulation, it's clear that the constellation points are $e^{j2\pi \frac mM},\, m =0,1,\ldots,M-1$, and putting that to $M$th power yields
$${\left(e^{j2\pi \frac mM}\right)}^M=e^{j2\pi m}=e^0=1,$$
due to the $2\pi$-periodicity of the complex exponential.
There you go, undergrad textbook-level data removal :)
Since your frequency offset is but a multiplied $e^{j2\pi f_\text{offset}t}$, you get a constant tone at $M$ times the offset frequency after.
Small problem: you don't have only signal, you do have noise in your received signal, which gets put to the $M$th power, too, and also intermodulated with the signal (you're converting nice additive noise to multiplicative noise). That typically drastically reduces your SNR – you don't get the most efficient frequency estimator that way. I'd advise against this whole approach!
The typical way to counter frequency offset in benign-sized modulation is to do a second-order PLL with loop filter tuned to the symbol rates you're expecting.
Add in that you'll of course have Gray mapping and forward error correcting codes, and that you can thus use your corrected decisions to refine your frequency error estimate, and you'll quickly see how undesirable "deleting" the info through nonlinear means (which are the only means to do that) is. For BPSK, the squaring approach is commonly done (because you'd want to do something similar to recover the timing, not only the frequency error, anyways), as it's low-complexity.
A second-order PLL works relatively well (even abused for QAM) until your constellation becomes large, where the probability that your PLL either "swallows" a legit symbol transition or misinterprets noise become problematic.
As everywhere, the synchronization of larger PSKs quickly becomes very system-specific and usually reflects the assumptions made on the channel. For example, if you have a satellite channel with substantial Doppler, and Doppler rate, but with low phase pertubations other than that, periodically sending a synchronization-aiding words from a smaller PSK might be a very good tradeoff in terms of complexity, robustness and data rate loss.


Why do we observe band when recording CRT monitors or Flouroscent lamp via video recorder? I understand that this has something to do with difference between refresh rate of monitor (or line frequency used by the ballast in case of Fluoroscent lamp) and frame capture rate of video recorder.

When learning about DSP in the past, I vaguely remember that a certain difference in frame capture rates of video recorder against the refresh rate of screen will cause the movement of bands upwards or downwards (assuming screen drawn left to right, top to bottom) ?

Please help me understand the phenomenon (the name, if there is one for this) or direct me to appropriate text on this.

Thank you.

This effect is most obvious in light generators that brighten and dim cyclicly over time, as in an AC current powered photon generator (bulb), or the phosphors on an old CRT periodically hit by a refresh beam.

In some ways this is similar to audio beating between two different but close audio frequency signals or pitched notes.  If the frequencies are slightly different (including the fundamental of one to a strong harmonic of another), then the phase will drift in and out over time.  In audio, this creates changing constructive or destructive interference.  In video recording this might capture the brightest portion of the beam or lamp when in phase, or not capture it when out of phase, or a capture over time while the phase is continually changing in between, resulting in rolling of the brightest band.

Note that the NTSC frame refresh rate is 29.97 Hz, which will beat against half an AC line frequency of 30 Hz.


Oversampling while maintaining noise PSD I'm designing a digital PID for control and stabilization of laser amplitude. The noise PSD of the laser intensity is around -140dB/Hz, and I want to servo that between DC and 100kHz. I think I've found a solution on the input side, but would like to make sure my understanding of noise spectral density vs SNR is correct.

I'm thinking about using an ADC at 5MSPS with a SNR of 100dB. My understanding is that this means that the power spectral density for a full-scale signal will have a noise floor of -163 dBFS/Hz, calculated form the SNR and the bandwidth. Since I only want to servo to 100kHz, I can use a digital FIR filter to maintain that noise floor so that I have a healthy amount of SNR between the noise that I want to servo and the intrinsic noise of the ADC.

What I'm a bit confused about is what the FIR filter does to the noise spectrum. I'm pretty sure that the filter will simply attenuate everything past the passband, and maintain the 167 dBFS/Hz - the SNR increases, but the effective bandwidth decreases. Is this actually correct?

One important point you may be missing is assuming that the noise of the ADC is truly white over your sampling bandwidth- it isn't! You will have spurs from signal correlations and ADC non-linearities that will limit the noise floor you can achieve. This will be specified by the SFDR (spurious-free-dynamic-range) of the ADC for guidance as to what you may be able to achieve. Since the location of the spurs will be difficult to predict, you should assume that the strongest spur could be within your bandwidth of interest, in which case no linear filtering will reduce the power of the spur which exists at a single point in frequency.

That said, the filter will reduce white noise components, and the best metric is the EQNBW (equivalent noise bandwidth) for your filter; which is the equivalent bandwidth of a brickwall filter that would remove the same amount of noise. In this case the SNR improvement (processing gain) is 10Log10(2EQNBW/fs) assuming a low pass filter where EQNBW is the bandwidth from DC to the equivalent cutoff frequency.


Ideal Discrete Time Hilbert Transformer with constant group delay Let's assume we have an ideal discrete time Hilbert Trasnformation system (90-degree phase shifter) with a frequency response over one period:

$$H(e^{jw}) =
\left\{
 \begin{array}{ll}
  -j  & \mbox{if}\;\;\;\;\;\; 0 < w < \pi \\
  \;\;j & \mbox{if}\;\; -\pi < w < 0
 \end{array}
\right.$$

I want to define an ideal frequency response $H_d(e^{jw})$ of an ideal discrete time Hilbert Trasnformation system that has a non-zero constant group delay.

I consider that $$H_d(e^{jw})= [1-2u(w)]e^{j(\frac{\pi}{2}-τ\cdot w)}$$ since $$\left | H_d(e^{jw}) \right | = 1 \; \; \; \forall w$$ and
$$\measuredangle H_d(e^{jw})  =
\left\{
 \begin{array}{ll}
  \frac{\pi}{2}-τ\cdot w  & \mbox{if} \;-\pi < w < 0 \\
  -\frac{\pi}{2}-τ\cdot w & \mbox{if}\;\;\;\;\;\;  0 < w < \pi
 \end{array}
\right.$$

So, it seems an acceptable answer. However, I didn't make any calculations to find it. Just guessed it. So my question is: which is the approach to such a problem is order to define $H_d(e^{jw})$?
Thanks in advance!

There's not much need for complicated calculations. If you have the zero-delay frequency response $H(e^{j\omega})$ as defined in the first equation of your question, the frequency response with a linear-phase (constant delay) is given by

$$H_d(e^{j\omega})=H(e^{j\omega})e^{-j\omega\tau}=e^{-j\left(\frac{\pi}{2}\textrm{sgn}(\omega)+\omega\tau\right)},\qquad \omega\in(-\pi,\pi)\tag{1}$$

where $\tau$ is the group delay in samples.


What is phase congruency? I am beginner in Image processing.I am studying the importance of phase in signal.Can anybody explain what is phase congruency ?

The more that sinewaves that have a zero crossing in the same direction at the same point add up in magnitude, the greater the contrast at that point.  e.g. a possible edge.  Phase congruency is another name (actually a superset) for the amount of sinewaves that line up with zero crossings at (nearly) the same point or points.

Whereas a random distribution of phases is more likely to result in something that a human would see as uniform noise rather than an edge.


What are some efficient algorithms for real-time calculation of eigenvalues? I need to calculate eigenvalues of a square matrix for a real-time implementation using microcontrollers or dSPACE. Obviously, I can find the eigenvalues of a 2by2 matrix analytically. I am wondering if there are any efficient numerical methods/algorithms for calculating eigenvalues of square matrices of dimension greater than 2?

Any help, including introducing a paper, is appreciated.

papers?

why not books

Golub, Gene H., and Charles F. Van Loan. Matrix computations. Vol. 3. JHU Press, 2012.

and if you need something more basic like understanding the difference between the Holder p-norm of a matrix and the Frobenious norm.

Stewart, Gilbert W. Introduction to matrix computations. Elsevier, 1973.

Although there are probably some exceptions, the symmetric square matrix is the most common in application in signal processing.

Typically, you calculate the eigenvectors and the eigenvalues follow.  I'm not aware of any algorithms that just calculate the eigenvalues.

The eigenvalues are sometimes called the matrix spectrum.

It would actaully be hard to not find a numerical linear algebra library that didn't compute eigenvector and SVD of a square matrix.

A popular book for  those who don't want to wade in deeply is

Press, William H., Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press, 2007.


How do I examine if the signal can be z-transformed? It's given signal

$x[n]=\sin(\frac{2 \pi}{N} m n) u[n]$

where $u[n]$ is the unit step function. Can I calculate Z-Transform?

$\mathcal{Z}$ transform exists when
$$ \sum_{n=-\infty}^{\infty} x[n]z^{-n} < \infty $$

Now I'm confused I would say it doesn't exist because it has an infinite number of elements and its sum is infinite but then I check table for transform and transform from sigma function exist? It just jumps into infinity. Can someone tell me where I am making mistakes..

You can use the basic properties of Z-transform to find the result:

Given, $$ x[n] = \sin(\frac{2 \pi}{N} m n) u[n] = \sin(\omega_m n) u[n]  $$

Decompose the sine wave into complex exponentials:

$$ x[n] = \frac{1}{2 j} \left( e^{j \omega_m n} - e^{-j \omega_m n} \right)   u[n] \tag{1}$$

And observe the following :

1 :
$$ U(z) = \mathcal{Z} \{ u[n] \} = \frac{1}{1-z^{-1}} ~~~~,~~~ |z| > 1 $$
2 :

if  $x[n] \longleftrightarrow X(z)$ then
$$ e^{j \omega_m n} ~ x[n]  \longleftrightarrow X( z ~e^{ -j \omega_m} ) ~~~,~~~\text{ ROC the same} $$

Then apply these two observations to Eq (1) :

$$ \begin{align}
X(z) &= \mathcal{Z} \{ x[n] \} = \mathcal{Z} \{ \frac{1}{2 j} \left( e^{j \omega_m} - e^{-j \omega_m} \right)   u[n] \} \\
&= \frac{1}{2 j} \left( \mathcal{Z} \{  e^{j \omega_m} u[n] \} - \mathcal{Z} \{ e^{-j \omega_m} u[n] \} \right)  \\
&= \frac{1}{2 j} \left(  U( z e^{-j \omega_m} ) - U( z e^{j \omega_m} ) \right)  \\
&= \frac{1}{2 j} \left(  \frac{1}{1-e^{j \omega_m} z^{-1}} - \frac{1}{1-e^{-j \omega_m} z^{-1}}  \right)  \\
\end{align}
$$


Why duplicating the signal length alters magnitude in FFT spectrum? I have a signal vector x=[a b c..] having 10 data values. I replicated this signal and created another vector y=[a b c.. a b c..] totaling to 20 data values.

But the Magnitude FFT spectrum results of signal x and y are different i.e. abs(fft(x)) and abs(fft(y)) are different. Any reasons what is happening.

This is due to Parsevals Theorem. The energy before and after DFT is equal (or depending on the normalization factor, at least proportional).

So, if you double the input energy by duplicating the input, you also get a different output energy. In your special case, you have a two-times repetition of the input, this corresponds to a discrete spectrum, where only every second bin is non-zero. Hence, the value on these bins needs to increase to fulfill the energy constraint.

Edit - More explanation on the discrete spectrum:
For the continuous-time Fourier Transform, every periodic time domain function of period $T$ has a discrete spectrum, where the spectrum consists only of discrete lines (Dirac impulses) with distance $1/T$.

A similar property holds for the DFT: if you input a periodic time function, you get out a discrete spectrum (i.e. it only every $n$th frequency bin is non-zero). Namely, if your time-domain input function has period of $N/L$, your spectrum is discrete, i.e. every $L$th bin is non-zero ($N$ is the DFT length).

In your case, your input function has period $N/2$, i.e. only every 2nd frequency bin is non-zero. Now, the double input energy (you have input 2 periods instead of one period) needs to be represented on these bins, leading to higher amplitudes compared to a single period.


Sampling Theorem illustration Can someone please explain the illustration (figure 1.19) at the bottom of this image? It looks like there are four sampling points but I don't understand what the different curves represent. Actually except from the dashed curve I can't even distinguish the different curves from each other. Can someone please clarify this so that the illustration makes sense?

(Image source: Digital Signal Processing Principles Algorithms and Applications)

This plot depicts how to convert your digital signal back to the analog one, using $\mathrm{sinc}$ functions. The nice property of these functions used in this process, is that maximum of each function always occurs at minimums of the other, shifted function:

Now the process of the D/A conversion that is depicted on your plot is basically taking a $\mathrm{sinc}$ function, and multiplying it by one sample that is at the maximum of $\mathrm{sinc}$. You do it for each and every sample of your digital signal, and then sum all of the outputs. For example, for the above plot, the outcome would be following:

Now maybe some more complicated example of D/A conversion. On the plot below, you can see the actual samples (in blue) and "analog" $\mathrm{sinc}$ interpolated signal (in red).

Plotting each function separately you can observe what will become an analog signal later:

I hope that it is more intuitive for you now. In order to complete this example, here is some MATLAB code I used:

clc, clear all, close all
fs = 32;
dt = 1/fs;

% Generate a samples of the signal
ts = linspace(0, 1, 10);
s = 2*sin(2*pi*ts) + sin(2*pi*3*ts);

td = linspace(0, (length(s)-1)*dt,     length(s));
t  = linspace(0, (length(s)-1)*dt, 100*length(s));

sincM = repmat( t, length(td), 1 ) - repmat( td', 1, length(t) );
SS =  s* sinc( sincM/dt );

% Sampled signal
stem(td, s, 'b', 'linewidth', 2)
hold on
% Reconstructed signal
plot(t, SS, 'r', 'linewidth', 2)
grid on

legend({'sampled signal', 'interpolated signal'})

figure
stem(td, s, 'b', 'linewidth', 2)
hold on
grid on
% Each sinc function
basis = sinc(sincM/dt)';
sbasis = bsxfun(@times, basis, s);
plot(t, sbasis)


Layman Description of the Kalman Filter I want to know about Kalman Filter but i tried searching different links including Electrical Engineering StackExchange but the information available there was hardly digestible.

All I am able to understand is that Kalman Filter is an estimator / predictor of future state of any system ,based on a model that is formed by previous measurements / readings.

So please kindly explain in simple words, what is Kalman Filter? What are its applications in Signal / Image processing and how is it different from other filters like Low pass, high pass etc?

KF is actually a mixture of a deterministic state propagator and a statistical estimator.

Despite it's name including the term filter, Kalman filter is not a simple frequency selective one. It's indeed a statistical recursive estimator of a state of a (linear) dynamic system. Yet on a broader sense it's called as a filter as it will separate a desired signal (the state) from corrupted measurements; hence filtering.

A linear dynamical system (LDS) is a physical construct whose behaviour is described by (linear) differential equations;

$$ \dot{X} = A X + Bu + W $$

where the vector $X$ is the state of the system, the matrix $A$ is the state propagation information (derived based on the differential equation of the LDS) , $u$ is a vector of deterministic inputs and $W$ is the state noise; depicting the uncertainity in the LDS propagation model $A$.

KF takes advantage of this deterministic behaviour of a LDS before estimating the current state $X$ purely from noisy measurements.

$$ Y = C X + Du + V $$

where $Y$ is the observation about the system and $V$ is the measurement noise on it.

Given a previous estimate of the state, it first deterministically propagates this previous estimate into current estimate through the matrix $A$ (process equation), and then performs a new finer estimate on the difference of the current measurement and propagated estimate using the measurment equation.


How does using DCT help in image compression I once learnt that the Huffman Encoding is used to rid of redundant information, thus reducing file size. I am not sure how using a DCT, or FFT on an image can help in reducing file size. It does not work on redundancy, then how does it reduce the total amount of data?

If you take an image as is (especially a photograph, as opposed to line art, a document scan...), there's very little redundancy in it for Huffman encoding to do its work - because the image contains many small fluctuations and details which are not perceptible to the human eye but nevertheless present in the data. DCT followed by coefficient truncation/quantization eliminates the barely perceptible details and creates redundancy for Huffman compression to work on.


The Parameter Devil — How to set them when no validation against groundtruth is possible Question:

I want to kick up a discussion on how people set algorithmic parameters when no validation against groundtruth is possible (maybe because groundtruth just cannot be obtained or is very hard/tedious to obtain).

I have read numerous papers and implemented the underlying algorithms wherein --- a set of parameters are said to have been set "empirically" --- and often i found that these are the ones that affect the generality of the algorithm (even though the theory underlying the method is elegant, enticing and sound).

I would appreciate it if you could share your thoughts. And, there is no right or wrong answer for this question. I just want to know, how everybody else deals with this.

Background/Source-of-question:

I am a computer scientist working in the areas of image analysis, computer vision, and machine learning and this question has been on the back of my mind for a while as i have faced this dilemma time and again whenever i design a new algorithm and i found myself spending a considerable amount of time tuning the parameters.

Also, I think, my question here is more general to any area where in computational algorithms are heavily involved, and i want to invite the thoughts of people from all concerned areas.

I wanted to give you some concrete example, just so it helps you think:

--- Take the case of feature detection (lets say circular blobs or salient points). You run some filters (needs parameters) at different scales (scale parameters) and probably threshold the response (threshold parameter). Its usually not possible to get a groundtruth to validate against and thereby automatically tune your parameters in such scenarios.

--- Take any computational framework which involves a lot of signal processing components. There are always parameters to tune and usually there is no groundtruth and when you tune them subjectively on a small random subset of your dataset you will someday encounter the case to which it doesnt generalize.

This parameter devil is more troublesome when you are setting parameters for some intermediate steps in your algorithm.

And I often found, it is not possible to cast the problem of finding good values for these parameters as an optimization problem with an objective function of which you can take a derivative and thereby use standard optimization algorithms to find good values.

Also, in many scenarios exposing these parameters to an end-user is not an option, as we often develop applications/software for non-computational end-users (lets say biologists, doctors) and they usually go clueless when you ask them to tune it unless its very intuitive (like approx object size).

Please share your thoughts.

Assuming that there is a ground truth, (at least theoretically) one of the possible ways to overcome the "tediousness" problem is a "bootstrap" ground truth creation. If you already have a decent algorithm that does the job in about, say 80%-90% of the cases, you can run your algorithm on a large set of instances and ask a user to mark only the mistakes. This approach has its own flaws, such as bias towards your algorithm.

However, there are some cases in which there is no ground truth at all, only different system trade-offs. For instance, an image processing system is required to output a sharp, color-accurate, non-noisy image. Obviously, you can't have them all at the same time. In such case you should use objective metrics that can be calculated on the result of your system. (See Imatest, DXO analyzer for image processing).

Once you have these, there are methods of multi-objective optimization that can create a mapping from the trade-offs (which are clear to the user) to the intrinsic parameters.

In any case, you should never give the user a parameter that he cannot understand. If all fails, simply hard-code the parameter.


how to find the inverse response of a system For the system described by the differential equation below find its inverse zero-state unit step response $$\dfrac{d^2y(t)}{dt^2}-2\dfrac{dy(t)}{dt}-8y(t)=\dfrac{d^2x(t)}{dt^2}-2\dfrac{dx(t)}{dt}-3x(t)$$
so that the cascade system of its inverse system and this system is the identity system.
I am not sure how to tackle this type of problem. My first thought is to switch $x$ and $y$ and find the response of that system using laplace transform, which I presume would be the inverse system. Is this reasoning correct? if not, how does one tackle this problem? I'm not asking for a solution (although one is welcome), but for a general approach to solve problems like this. Thanks!

Dan's answer -- to compute $H(s)$ as normal, and then compute $1/H(s)$ -- is equivalent to your suggestion of swapping the $x$ and $y$ (or doing it in one step by solving for $H_{yx}(s) = \frac{X(s)}{Y(s)}$).
In the Laplace domain it's justified by noting that
$$\frac{H(s)}{H(s)} = H(s)H_{yx}(s) = 1$$
In theory this means that a system followed by its inverse system has an output equal to its input.  In practice you need to make sure that $H_{xy}$ is realizable, and you have to recognize that if you're actually chaining real systems, then your $H(s)$ is just a model that may not capture all of your real-world system's real-world dynamics.  So the actual end-to-end equation is something like
$$\hat X(s) = \left(H(s) X(s) + N(s)\right) \hat H_{yx}(s)$$
where $\hat X$ denotes that it's just an estimate, $N(s)$ is whatever noise is injected in your measurement of $Y$ and your implementation of $\hat H_{yx}$, and $\hat H_{yx}$ emphasizes that it's a guess at your system model, not the real thing.
It's also useful to note that doing this in the Laplace domain immediately lets you know whether what you're trying is practical.  If $H(s) = 1 / (\tau_0 s + 1)$, for instance, you'd immediately see that $1 / H(s)$ is, strictly, unrealizable (because you can't have a true naked differentiator).
However, as long as noise is negligible you could approach it arbitrarily closely by setting $\hat H_{xy}(s) = \frac{\tau_0 s + 1}{\tau_1 s + 1}$, which basically says you're approximating a derivative as $\frac{s}{\tau_1 s + 1}$ -- then you could use my end-end equation above, along with your knowledge of the nature of $X(s)$ and an estimate of the noise you're going to inject at $N(s)$ to estimate your end-end system fidelity.


Teager's Energy Operator Does anyone know how to find information on the derivation (rather than the definition and application) of Teager's Energy Operator, both for the continuous and sampled case.
I have Kaiser's early papers "On a simple Algorithm to calculate the 'Energy' of a signal" and "Some useful properties of Teager's Energy operator" but these seem to stop short of giving the derivation.

This paper does a good job of explaining it, based on the energy of a mass-spring system.

R.C. Guido, "Enhancing Teager Energy Operator based on a novel and appealing concept: signal mass", J. Franklin Institute, 2019. https://doi.org/10.1016/j.jfranklin.2018.12.007


Answering Machine vs. Human -- Neural Network Features Selection I've been tasked with creating an artificial neural network that can classify a telephone call as either answered by a human or an answering machine.
My knowledge of audio processing is, mildly put, limited. I'm quite new to machine learning as well, so feel free to question everything.  Here's what my run-of-the-mill backpropagation network looks like so far:
1. 37 input units;
2. One layer of 10 hidden units;
3. One output unit.

An audio file comes in, and I extract the data like so [X, fs] = wavread(filename). Then, the 37 features I obtain are the following:
1. The length of X;
2 - 21. In 5% increments, the mean of X in each interval (I know);
22. Spread;
23. The mean of short-time energy;
24. Zero crossing rate;
25 - 37. The mean of each in a set of 13 MFCCs.

Now, since I've obviously jumbled it all together, my question is the following: What do you think would be a good set of features for the problem at hand?

Thank you for your time

I don't think that you would do this by feature extraction. Instead, I would train a many-to-one LSTM, which classifies the input sequence as -1 or 1. For you, it is easy to collect a huge amount of labeled speech data, just by calling. Note that, the data size is important.

LSTMs are very powerful techniques and should handle this problem gracefully.


Processing and feature extraction from very short signals I recently have encountered very short time signals (from three to 12 samples, typically). They appeared for instance in small time frames in engine simulation that should be extrapolated in a robust manner, or gaz measurements taken for instance at $t_0 = 0$ (seconds),  $t_1 = 5$,  $t_2 = 10$, $t_3 = 20$, $t_4 = 50$, $t_5 = 100$. The idea behind this time sampling is that the signal should rise or fall, and may reach some steady-state. Such data cumulates caveats:

Very short,
Sometimes corrupted by one missing data, or one potential outlier,
Uneven sampling.

Traditional tools (Fourier, correlation, fitting) appear of limited use, and are probably unusable. Yet I would like to:

describe those signals with few robust parameters,
estimate offsets and delays between two signals,
cluster them with shape similarity.

I for instance would like to detect "V-shaped" short signals that first go down then up,  obtain a location estimator for the wedge, and separate them from monotonous or down-then-flat signals.

So far, I am performing tedious 3- or 4-point linear or parabolic fits, estimate residues and dispersion to select "a best model", with standard quadratic or robust quantile estimation. My questions are:

Is there recommended literature on very-short signal processing?
Would users share their best-practice experience on such signals, that ought to be processed?

Though problem! However, signal processing might have a tool at hand: It's called compressive sensing, and reduces the number of samples you need to sub-Nyquist-rate levels.

It's a bit nonsensical to derive the math behind compressive sensing for signals that are a sum of sinusoids (that being the only thing that I could do without literature) if that doesn't describe your signals overly well. However, the idea is:

If you want to know about everything happening with less than a given frequency over a specific time, you'll need all of the samples at a constant rate (typically, twice the highest frequency for real-valued samples). Then you can, for example, do a base transform to get the samples from time domain to frequency domain using the DFT.

For the following, assume that you could, if you had the "full, long, evenly-sampled" signal $\mathbf S \in \mathbb R^N$, define a base transform $\mathbf T$ to transform the time samples vector to a result vector giving you the information you want $\mathbf Y = \mathbb R^N$. That base transform would then simply be a matrix:

$$ \mathbf Y = \mathbf S^T \mathbf T$$

If what you want is sparse, ie. the dimensionality of the result vector would be much smaller than the observed samples vector, the transform matrix wouldn't be a square full-rank base transform, but a matrix that would have $N$ rows, but fewer columns.

The point here is that most of the values in $\mathbf T$ would be pretty small -- setting those to zero will incur an error, that "feels" like noise, but given a nonzero upper limit for allowable distortion, you can find a sparse matrix that still can allow for reliable conclusions about $\mathbf Y$.


How to extract a smooth contour from a set of points in 3D? I am trying to segment a 3D volume. The outcome of all my volume segmentation algorithms is a set of candidate points in 3D space. Now I need to smooth this point cloud and fit a closed surface it. My question is in two parts:

1) How can one smooth a set of rough points to get a set of soft points when demonstrated as a surface?

2) How can I fit a surface to these points to get a closed smooth shape?

From your description it sounds like surface reconstruction methods such as Poisson Reconstruction would fit your description. These algorithms already have the necessary smoother built-in, mainly in the form of a regularizer. To have slightly more accuracy, in general I prefer the following algorithm:

  SSD: Smooth Signed Distance Surface Reconstruction, F. Calakli and G. Taubin, Computer Graphics Forum, Vol. 30, No. 7, 2011, http://mesh.brown.edu/ssd/

There is a nice Adaptive Solvers repository containing many popular methods under a single roof:
https://github.com/mkazhdan/PoissonRecon

I strongly recommend to take a look at this. Here is how you would like to use it:

Export your point set and the surface normals to a 3D file : sample.ply
Run (make sure to play with the depth value):

SSDRecon --in sample.ply --out sample_reconstruction.ply --depth 10 --density

This will generate a reconstruction but if you have an open surface the reconstruction will contain undesired extrapollations. So move to (3). If the surface is perfectly closed, then the output will already be good in that stage.

If you have an open surface then run:

SurfaceTrimmer --in sample_reconstruction.ply --out sample_reconstruction_trimmed.ply --trim 7

Again, make sure to play with the trimming threshold. This will give you very good reconstructions such as (taken from here):

Note that for all these algorithms to work, surface normals are essential. Luckily, they can be computed from the data using local plane fits. Moreover, one also needs to ensure the correct signs of the normal directions to properly orient the surface. In the case that a single camera captures the data, the canonical choice is to flip all the surface normals so that they make an acute angle with the viewing direction.


File size adjustment using FFT I have been given the task of doing audio compression using FFT.
The assignment is below.

Below is my MATLAB code:

% part 1
[y, f] = audioread('file.wav');
y = y(1:1024*161);
% p = audioplayer(y,f);
% play(p)
Ya = fft(y);
Y = fftshift(abs(Ya));
N = length(Y);
w = -f/2+f/N:f/N:f/2;
figure(1)
plot(w,Y)
% part 2
Y1 = Ya((N/4):(3*N/4));
% part 3
y = ifft(Y1);
player = audioplayer(y,f);
play(player)

Now here is my problem. I was able to do part 1 and part 2 perfectly with no problem, but for part 3 the problem begins. I get no sound at all!!. Now in the prompt it asks to do FFTSHIFT before the IFFT, to me that does not matter, either how I get now sound. can somebody help me on how to fix this please?

The algorithm sets the high frequency DFT coefficients to zero hence effectively providing a lowpass filtering on the time domain signal. Your problem was on the line of creating Y1, which I replaced with the corrected form. I've changed the code a little bit also.

% part 1:
[y, f] = audioread('C:\path to your wav file');
y = y(1:1024*64);                % select a portion of it.

Ya = fftshift(fft(y));           % Compute DFT Ya[k] and SHIFT it
N = length(Ya);                  % Length of DFT Ya[k]

figure,plot(-f/2+f/N:f/N:f/2 , abs(Ya))  % plot the DFT magnitude.

% part 2:
Yr = zeros(1,N);                % Yr[k] of length N!
Yr((N/4)+1:(3*N/4)) = Ya((N/4)+1:(3*N/4));  % Assign nonzero coefficients

% part 3:
yr = real(ifft(fftshift(Yr)));  % RECONSTRUCT the filtered signal.
                                % Take REAL part for convenience.
figure,plot(-f/2+f/N:f/N:f/2, abs((Yr)));

% part 4:
sound(yr,f,16);                 % Listen to the reconstrcusted audio


Kalman Filter implementation is too slow I have implemented with a simple code a Kalman Filter for time domain, based on these:

KG = error_est / (error_est+error_measurment)
estimate = estimate_prev + KG ( measurment -  estimate_prev)
error_est = (1 - KG) * error_est_prev

The simple code to implement is :
error_est = 1.0
error_measurment = 1.0
estimates = []
start = 50.0

# loop over data array of numbers such as [1,2,3,4,..]

for index in range(len(data)):

    measurment = data[index]
    prev_estimate = estimates[index-1] if index>0 else start

    KG = error_est / (error_est+error_measurment)
    estimate = prev_estimate + KG * ( measurment -  prev_estimate)
    estimates.append(estimate)
    error_est = (1.0 - KG) * error_est

Now I have 3 questions,

The filter is extremely slow, which mean for array of [1...50] , I will get a result of [1...25], and generally it is much slower than moving average, but i read that it's advantage is that it is fast. What am i doing wrong?

I read that Kalman is used to estimate next result, but how is it possible that a filter ( based on the past) can be used to estimate the next result ?

Are there any parameters you can add to tweak it (such as window in MA) ?

Maybe what is extremely slow is python, if you are comparing with a rolling average function implemented in some library.Specially if the dimension of the observations large, each step it will invert a matrix, $O(n^3)$, it seems you have only one observation in your example. Even for similar implementations Kalman filter will be slower than moving average.
The Kalman filter uses the model of the system, from one state it can get a reasonable estimate, for instance if you are tracking a device with GPS data, you may create a model that assumes that the acceleration has limited variance, from the previous measurements of position you can have an estimate of velocity (part of the state), and if from there you can predict where the device will be.
The things you can tweak in Kalman filter is the covariances of the process and the covariances of the observation, this will tell how much the filter should rely in the estimate vs. the observation. For instance if you don't trust your model you can set the model covariance $Q$ large, if your observations are not that accurate you acn increase the observation covariance $R$.


Identify Where Singing Starts in a Voice Only Recording I have a little karaoke-style app where a user sings 4 lines of a song, with a one second gap between each line.  There is no backing music, so it's voice only, hopefully making the problem easier to solve.

I am looking for the most robust way to detect exactly where in my recording the user starts and ends singing line 1, starts and ends singing line 2, etc.

I have cobbled together a simple-minded algorithm that works when there is very little background noise in the recording (like when does that happen?), but it falls to pieces in the presence of the smallest noise.

Can anybody point me towards something more robust?

If the background noise is white-ish, you could measure spectral flatness and consider it to be voice when the amplitude is above some threshold and the spectral flatness is below some threshold.

Basically you just take an FFT of a chunk of the signal, then divide the geometric mean of the magnitude of the spectrum by the arithmetic mean.

You could also use a band-pass filter to emphasize only the frequency regions where the human voice usually sits (as simple as setting the unwanted regions of the FFT to 0 before measuring the spectral flatness)


Removing periodic spike noise from ECG signal The signal shown in the following figure is collected from a ECG sensor. The spike noise that is observed with a periodicity of 30 seconds was traced to the periodic blip of the LED as it draws current.

The following is one instance of the spike noise after zooming in,

Frequency domain filters are not effective in this case because the frequency domain characteristics of the LED blip artifact is similar to the heart beat signal(QRS).

Is there any effective way to remove this artifact..? I was think of using the blip as a mother wavelet and performing a multi-scale correlation with the ECG signal to detect the spikes. Is this the procedure for wavelet denoising? Is there any other effective way of denoising for this case?

In the first place I would suggest you to look for some electrical methods of avoiding those spikes... If it's impossible to do so, then you must find a close enough mathematical model of your noise generation process from which you may decide to choose an adequate method of noise removal...Provided you can do so...


What fps can be considered as a hardline for real-time performance? Is there any academic paper that describes this? I am trying to design a deep learning based inferencing solution for security applications. My deployed program achieves an FPS of 15fps for classification. Can it be considered real-time?

There are as many ad hoc "definitions" of "real-time" are there are people or places you hear or read about what it supposedly means. Some people claim there is a "hardline" such as 10uS response time, but I believe there is no academic paper that makes that mistake.

Many people implicitly have an informal mental model that considers a system as being "real-time"
•   if, or to the extent that,
•   its actions are manifest to action observers with latencies that can be related (usually in an application-specific way) to the actions’ perceived currencies—
•   i.e., in a time frame that those actions’ latencies and predictability of latencies have acceptably satisfactory value to their observers (e.g., system users). Think of high frequency trading.

The magnitudes of the time frames are application- and situation-specific—e.g., microseconds to megaseconds—not related to if, or to the extent, that a system is a real-time one.

So a system "operates in real-time" if its actions have low enough latency to satisfy your needs. If 15 FPS is a satisfactory enough rate, then the system is real-time as far as your needs are concerned; otherwise it is not--that's all that matters to you.


numpy.linalg.lstsq underdetermined case I would like to understand what I am doing wrong here.
I am trying to perform polynomial regression by minimizing the least squares, ||Au-y||^2, where y is the given data and A is the matrix where the i-th line holds [1, x_i, (x_i)^2, ... (x_i)^n-1].

If n is larger than the number of data points, the problem is underdetermined, and I expect the numpy.linalg.lstsq() routine to give any of the infinitely possible solutions. But, as you can see, I don't get a solution at all.

import matplotlib.pyplot as plt
import numpy as np

x = np.array([-6 ,1, 2, 3, 4])     # x data
y = np.array([2, -3, 4, 20, -10])  # y data

A = []
n = 50      # polynomial degree
for i in range(0,n):   # create A matrix of proper form: i-th line is [1, x_i, x_i**2, ...]
    A.append(x**i)
A=np.array(A)
A=A.T

u=np.linalg.lstsq(A,y, rcond=None)[0]  # solve underdetermined problem

x_test=np.linspace(-6, 5, 100)   # create more x values for plotting
B=[]
for i in range(0,n):   # same as before, their power matrix
    B.append(x_test**i)
B=np.array(B)
B=B.T

plt.plot(x_test, B@u)   # plot polynom
plt.scatter(x,y)   # plot data
plt.show()

print(A@u-y)    # this should be zero vector??
```

I ran your code and got the below plot as well as an almost zero vector of:

[-1.19209290e-07, -2.42083020e-09, -8.94069672e-08,  2.23517418e-07, -6.33299351e-08]

This means numpy successfully curve-fit a line to your data.

An order of 6 works a lot better, though:

[-2.79396772e-09 -2.04281037e-14 -1.27897692e-13  4.33431069e-12 -2.27373675e-11]

numpy won't give you an infinite number of solutions.  It will give you an optimally curve-fit solution to your data for the given polynomial order.


Comparison of size of an entity in different images with different resolutions? Can we compare the size of two objects appearing in different images with different resolutions?

You must have some reference between the images,
either some other similar elements or the camera parameters,

If you know the different resolutions that it is easy to calculate:

relative size = resolution ratio * object_size_ratio

otherwise you must have some other spatial reference between the images.


Identify random repetitive patterns Forgive me if it’s too basic, I finish engineering a while ago.
Given any time series, not periodic, I would like to find any repetitive pattern that is distinct (by some given measurement) and is unknown.

For example given this :
 I would like to find number 1 and 2( get their index on time).

Is FFT going to help anyway?
Without knowing or hard programming an algorithm, is it possible or do I need a serious machine learning?

(Are there practical books about finding patterns with code and not too much math?)
Thanks.

Finding repeating but not periodically repeating patterns of an unknown template which you expect the algorithm to identify is a hard problem.  Nature typically doesn't shut off all the other patterns that may be present for our convienance either.

Signals have sources or generating mechanisms. It would help if you had some idea how your patterns are generated.

The concept of noise can be helpful because what isn't noise might very well be signal.

I would approach a problem as a time series that corresponds to an ODE and try to learn the ODE.

One way to do this would be to construct a phase space representation, that is something like $x(t)$ and $y(t)=dx(t)/dt$. If the noise is low using central differences can give a good approximation to the derivative of your time series. A phase space plot might identify orbits and limit cycles produced by your pattern.

These ideas are common in chaos theory and I would search Google accordingly

As an example https://en.m.wikipedia.org/wiki/Van_der_Pol_oscillator


What is the relation between hilb and Hilbert transform in MATLAB? In MATLAB, there is a command called hilb(n) which it says gives a Hilbert matrix of order n.  Is this matrix the one used for Hilbert transformation?

No, the Hilbert matrix has nothing to do with the Hilbert transform, other than that both are named after the famous German mathematician David Hilbert.


How to filter ECG and detect R peaks I am trying to extract R peak from raw ECG data and some samples are seemed to be interfered by EMG.
I used the lib provided by biosppy with python, biosppy.signal.ecg, it seems that the lib does,

Default filter is firwin with filtfilt, mode is bandpass with 3-45Hz, the filter order is 0.3 * sample_rate, in my case, sample_rate is 500Hz, and order is 150
The R peak extract algorithm is P.S. Hamilton, "Open Source ECG Analysis Software Documentation", E.P.Limited, 2002

The result seems not stable when handling some samples interfered by EMG,

Between 2 and 8, the R peaks seem to be not correct.
I'm not familiar with DSP and what can I do to get the raw ECG from the composed signal?

What I tried(I read some papers for ECG filtering) and the issue encountered,

Baseline Wander, some papers suggest to use firwin with highpass filter to remove this kind of noisy, and the frequency can be 0.5Hz, I tried this both with firwin and butter, the filtered signal was more smooth than before, but the EMG interfered still exists
Powerline Interference, the bandstop frequency is 50Hz or 60Hz in most cases, since the firwin using frequency from 0.5Hz to 45Hz, this kind of noisy also filtered
EMG Noisy, some paper mention to use moving average to remove this kind of noisy since this noisy is complex and may overlap with the raw ECG, some suggest to use a 8 window moving average. But in my case, the noisy last a long time with many samples, even the window set to 150, the EMG also exists

What I want is to get the stale R peaks without the interference? if the EMG is not easy to remove, can I remove the invalid R peaks and do interpolation after R peaks detected, such as,

If the R peaks detected is,

[439, 433, 433, 245, 193, 150, 414, 307, 135, 300, 432, 430, 429]

The first 3 and last 3 indices seem to be correct with the correct raw PQRST wave, but other indices are the invalid ones interfered by EMG, remove these indices, and insert 3 new indices with the value 433?

Thanks for your kindly help.

You should use py-ecg-detectors
Siply install by doing

pip install py-ecg-detectors

Then you can use for instance the well known Pan Tompkins algorithm to find the R-peaks

Here I used an ECG recording from the Apnea-ECG Database

from ecgdetectors import Detectors
import matplotlib.pyplot as plt
import pandas as pd
fs=100 # sample freq

heartbeat = pd.read_csv("/myworkdir/mitbih-database/100.csv")

detectors = Detectors(fs)

r_peaks_pan = detectors.pan_tompkins_detector(heartbeat.iloc[:,2][0:1000])
r_peaks_pan= np.asarray(r_peaks_pan)

plt.plot(heartbeat.iloc[:,2][0:1000])
plt.plot(r_peaks_pan,heartbeat.iloc[:,2][0:1000][r_peaks_pan], 'ro')

If you want you can also add a R-peak correction algorithm:

def R_correction(signal, peaks):

    num_peak=peaks.shape[0]
    peaks_corrected_list=list()
    for index in range(num_peak):
        i=peaks[index]
        cnt=i
        if cnt-1<0:
            break
        if signal[cnt]<signal[cnt-1]:
            while signal[cnt]<signal[cnt-1]:
                cnt-=1
                if cnt<0:
                    break
        elif signal[cnt]<signal[cnt+1]:
            while signal[cnt]<signal[cnt+1]:
                cnt+=1
                if cnt<0:
                    break
        peaks_corrected_list.append(cnt)
    peaks_corrected=np.asarray(peaks_corrected_list)
    return peaks_corrected

corrected_R_peak=R_correction(heartbeat.iloc[:,2][0:1000],r_peaks_pan)
plt.plot(heartbeat.iloc[:,2][0:1000])
plt.plot(corrected_R_peak,heartbeat.iloc[:,2][0:1000][corrected_R_peak], 'ro')


How could I approach determining if this 2D system represented as a 2D summation formula is linear? I have a given 2D system:
$$y(m,n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} x(k_1,k_2)$$
My usual approach to determining if a system is linear is to test if it is homogeneous and additive. First I would input an arbitrary point: (0,0), (1,1), (0,1), (1,0), calculate its output. Then I would scale the input and scale the output and see if the scaling amount matches. Next I would test the additive property in a similar manner. This is fine for a simple 2D system formula, for example: $y(m,n) = x(m,n) + 3$
With a 2D system represented as summation, I am not sure how to apply this test.
Searching the web, it seems this system looks like a 2D convolution:
$$y[i,j]= \sum_{m=-\infty}^{\infty} \sum_{n=-\infty}^{\infty} h[m,n]⋅x[i−m,j−n] $$
However, I am still stuck thinking of how to prove this system is linear. I would assume the process is fundamentally the same (prove it is homogeneous and additive or not). I am just not sure of the steps.

I'm actually starting this out not knowing the answer, although I have my suspicions.  I think it'll be a linear system, but I'm not sure.  Here's your starting point: you define the system $y = h(x)$ as:
$$y = h(x) \mid y(m,n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} x(k_1,k_2)$$
So let
$$\begin{aligned}
y_1(m,n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} A\,x_1(k_1,k_2) = h\left(A\, x_1 \right) = A h(x_1) \\
y_2(m,n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} B\,x_2(k_1,k_2) = h\left(A\, x_2 \right) = B h(x_2)
\end{aligned}$$
Then
$$
y_1(m,n) + y_2(m, n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} A\,x_1(k_1,k_2) + \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} B\,x_2(k_1,k_2)
$$
Because the summations have identical limits, we can bring everything under one pair of summations:
$$
y_1(m,n) + y_2(m, n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} \left ( A\,x_1(k_1,k_2) + B\,x_2(k_1,k_2) \right)
$$
but this is just $y_1 + y_2 = h(A\,x_1 + B\,x_2)$.  You already know that $y_1 = A h(x_1)$ and $y_2 = B h(x_2)$, so linearity is proved.
(Which is what I thought -- I was going to just dismiss your original problem as "summations are linear, go away!" -- then I realized I wasn't sure).


How to distinguish between recorded and a live vocal signals? Is it possible to detect if a sample of human voice by smartphone/computer(44Khz) is sourced from another electrical device or originated  from a human being

I want to build a model in matlab such that, samples a vocal signal and decides whether it is a real human being talks a live or recorded, is it even possible?

On a smartphone, flat out no. The microphones are so poor that you're lucky to get any approximation of the actual audio.

On high-quality hardware, perhaps. Barely. If you're lucky. Check for MP3 compression artifacts. If the source is not compressed with a lossy codec, you're out of luck.


What factors should I consider in choosing an edge detection algorithm? I've learned about a number of edge detection algorithms, including algorithms like Sobel, Laplacian, and Canny methods. It seems to me the most popular edge detector is a Canny edge detector, but is there cases where this isn't the optimal algorithm to use? How can I decide which algorithm to use? Thanks!

There are lots of edge detection possibilities, but the 3 examples you mention happen to fall in 3 distinct categories.

Sobel

This approximates a first order derivative. Gives extrema at the gradient positions, 0 where no gradient is present. In 1D, it is = $\left[ \begin{array}{ccc} -1 & 0 & 1 \end{array} \right]$

smooth edge => local minimum or maximum, depending on the signal going up or down.
1 pixel line => 0 at the line itself, with local extrema (of different sign) right next to it. In 1D, it is = $\left[ \begin{array}{ccc} 1 & -2 & 1 \end{array} \right]$

There are other alternatives to Sobel, which have +/- the same characteristics. On the Roberts Cross page on wikipedia you can find a comparison of a few of them.

Laplace

This approximates a second order derivative. Gives 0 at the gradient positions and also 0 where no gradient is present. It gives extrema where a (longer) gradient starts or stops.

smooth edge => 0 along the edge, local extrema at the start/stop of the edge.
1 pixel line => a "double" extremum at the line, with "normal" extrema with a different sign right next to it

The effect of these 2 on different types of edges can be best viewed visually:

Canny

This is not a simple operator, but is a multi-step approach, which uses Sobel as one of the steps. Where Sobel and Laplace give you a grayscale / floating point result, which you need to threshold yourself, the Canny algorithm has smart thresholding as one of its steps, so you just get a binary yes/no result. Also, on a smooth edge, you will likely find just 1 line somewhere in the middle of the gradient.


Upsampling with time offsets Suppose I have done 4x oversampling for a continuous time signal, but the successive sampling times have a linearly increasing offset. Specifically, the samples with indices {4k; k=0, 1, 2,...} are correctly sampled at times {4kTs}, but the samples at indices {4k+1; k= 0, 1, 2,..} are sampled at {(4k+1)Ts+e}, the samples at indices {4k+2; k= 0, 1, 2,..} are sampled at {(4k+2)Ts+2e)}, and the samples with indices{4k+3; k= 0, 1, 2,..} are sampled at {(4k+3)Ts+3e)}. How does the spectrum of this oversampled signal look ?

Answer: You will see residual images of $X(f)$ at multiples $f_s$, $2f_s$ and $3f_s$, and distorted image of $X(f)$ at non-zero multiples of $4f_s$, when sampling in the manner you explained. Depending on value $e$, the size of residual will change. I have explained how in detail below.

Ideally, sampling at $4f_s$ would have completely cancelled those images of $X(f)$ at multiples of $f_s$, $2f_s$ and $3f_s$ and you would've seen images of $X(f)$ only at multiples of $4f_s$ and the magnitude scaled by $4f_s$.

(Explanation seems long only because I have included lots of pictures to show. Please follow through.)

I would like to give you an intuition on how to visualize sampling at any rate $f_s$. You probably have a pretty good idea about that. But then I would like to show pictorially what happens when you sample at $2f_s$ and then you can extend the idea to $4f_s$.

  Sampling a bandlimited signal $x(t)$ at sampling rate $f_s$ :

When you sample $x(t)$ at sampling rate $f_s$, you are basically, multiplying $x(t)$ with a periodic pulse train with period $T_s = \frac{1}{f_s}$ in time domain. Hence, in frequency domain you see a convolution of $X(f)$ with Fourier representation of that periodic pulse train.

Since the pulse train is periodic, it's Fourier representation will be obtained by computing Fourier Series. The pulse train can be represented by: $$\sum^{\infty}_{k=-\infty} \delta(t - kT_s)$$

 and it's Fourier transform as: $$\frac{1}{T_s}\sum^{\infty}_{k=-\infty} \delta(f - kf_s)$$

Notice that the magnitude of Fourier domain representation of the sampling pulse train has scaled up magnitude of $\frac{1}{T_s} = f_s$.

Assume, frequency representation of $x(t)$ as below :

Now, as described above, Sampling is nothing but multiplication of $x(t)$ with a $T_s$ periodic pulse train in time domain and hence a convolution of $X(f)$ with Fourier transform of the pulse train in frequency domain. Mathematically,
$$x(kT_s) = x(t).\sum^{\infty}_{k=-\infty}\delta(t-kT_s)$$
$$X_{sampled}(f) = X(f) * f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$

  Now, lets explore what happens when we sample at $2f_s$ or when the sampling pulse train becomes $\frac{T_s}{2}$ periodic :

The sampling pulse train becomes as shown below :

The sampling pulse train can be represented by a sum of two $T_s$ periodic pulse trains:

One pulse train as if $T_s$ periodic and centered at 0.
Second pulse train as if $T_s$ periodic but shifted by $\frac{Ts}{2}$.

Hence, mathematically it will be as follows:
$$\sum^{\infty}_{k=-\infty}\delta(t-k\frac{T_s}{2}) = \sum^{\infty}_{k=-\infty}\delta(t-kT_s) + \sum^{\infty}_{k=-\infty}\delta(t-kT_s - \frac{T_s}{2})$$

Hence, the Fourier Transform of the sampling pulse train will also sum of these 2 trains, as convolution is a linear operation. Also, use the time shift property of Fourier transform to get the result as follows:

$$\mathcal F \{ \sum^{\infty}_{k=-\infty}\delta(t-kT_s) + \sum^{\infty}_{k=-\infty}\delta(t-kT_s - \frac{T_s}{2}) \}$$
$$= f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s) + f_s.e^{-j\pi \frac{f}{f_s}} \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$

$$= f_s. \sum^{\infty}_{k=-\infty}\delta(f-kf_s) + f_s.(cos(\pi\frac{f}{f_s})-\mathbb i.sin(\pi \frac{f}{f_s})). \sum^{\infty}_{k=-\infty}\delta(f-kf_s)$$

Notice that the sum is evaluated only at integral multiples of $f_s$, because of the $\delta(f - kf_s)$. What this means is that $sin(\pi \frac{f}{f_s})$ will always be $0$, so, no imaginary images of $X(f)$ will be seen, and $cos(\pi \frac{f}{f_s})$ will be $1$ at even multiples of $f_s$ and $-1$ at odd multiples of $f_s$. Pictorially, the fourier transform of pulse train which is used to sample at $2f_s$ will look like following :

So, the even multiples of $f_s$ will be doubled in magnitude to $2f_s$ and odd multiples of $f_s$ will cancel out each other to cancel the images of $X(f)$.
This is the reason you see images of $X(f)$ only at multiples of $2f_s$ when sampling at double the rate, because images of $X(f)$ at odd multiples of $f_s$ cancel each other out.

Now, consider the case which you have explained in your question. When you break your pulse train into 4 pulse trains which are $T_s$ periodic individually, but shifted as below:

$\sum^{\infty}_{k=-\infty}\delta(t-kT_s)$
$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{T_s}{4} - e)$
$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{2T_s}{4} - 2e)$
$\sum^{\infty}_{k=-\infty}\delta(t-kT_s -\frac{3T_s}{4} - 3e)$

Conclusion:

When you check their Fourier transforms, you will find that images at multiples of $f_s$, $2f_s$ and $3f_s$ will not get cancelled completely because the negative impulses (both real and imaginary) are shifted by $e$, $2e$ and $3e$ respectively. And, image at multiples of $4f_s$ will also not be aligned exactly to give a scaling of $4f_s$ but they will be fudged around to give a distorted image of $X(f)$ except at $k=0$, that is the original image of $X(f)$ centered around DC.

Depending upon the value of e, the real and imaginary images of $X(f)$ will have residuals at $f_s$, $2f_s$ and $3f_s$, and images at non-zero multiples of $4f_s$ will be fudged around.


Filter design - Rules of thumb regarding ratio of sampling and stop/pass frequency The issue: I have a signal in the time domain sampled at 51.2KHz. I want to apply a highpass filter with a pass frequency of 20 Hz.

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

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

In the most general case, there isn't really a formula for how long the filter should be. There are papers which show the empirical relationships between the filter transition width, passband ripple, and stopband attenuation. You typically have - for a fixed sampling frequency:

A smaller transition width (from the stopband to the passband), as a percentage of the sampling frequency, requires a longer filter. A 100 Hz transition width for fsamp=5kHz would require a smaller filter than a 100 Hz transition width with fsamp=15kHz - assuming other requirements are equal.
Smaller ripples in the passband, requires a longer filter
More attenuation in the stopband requires a longer filter

That said, for remez type designs there are a couple of published formulas that give you an estimated length of the filter. In Matlab you could use the remezord() function. These formulas work fairly well, but they do run into problems when the transition frequencies are near 0 Hz or fs/2.

There are some published papers that give equations for the length of filter using a Kaiser window for a windowed Sinc design - but you don't have as much control in the specification of the filter - if I recall correctly.


Inverse DTFT of $H_1(\Omega)=\begin{cases} 10,& \frac{\pi}{3} \leq |\Omega| < \pi\\ 0,& 0 \leq |\Omega| < \frac{\pi}{3}\\ \end{cases}$ What is the inverse DTFT of the $2\pi$-periodic extension of following function:

  $$H_1(\Omega)=\begin{cases}
10,& \text{for } \frac{\pi}{3} \leq |\Omega| < \pi\\
0,& \text{for } 0 \leq |\Omega| < \frac{\pi}{3}\\
\end{cases}$$

I have found out using the definition that it is: $$h[n]=-\frac{10\sin(\frac{\pi}{3}n)}{\pi n}$$

But the problem is that this has a DTFT of

$$H_2(\Omega)=\begin{cases}
-10,& \text{for } 0 \leq |\Omega| \leq \frac{\pi}{3}\\
0,& \text{for } \frac{\pi}{3} < |\Omega| < \pi\\
\end{cases}$$

Now, $H_1(\Omega)$ and $H_2(\Omega)$ seem quite different to me. Am I missing something? Have I done something wrong in my calculations?

Let's use $\omega$ instead of $\Omega$ for the definition of the inverse DTFT and state the following:

$$h[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} H(\omega) e^{j \omega n} d\omega $$ Given that $H(\omega) = 1 $ , for $\pi/3 < |\omega| <\pi $ and zero otherwise, the integral becomes:

$$\begin{align}
h[n] &= \frac{1}{2\pi} \int_{-\pi}^{-\pi/3} e^{j \omega n} d\omega + \frac{1}{2\pi} \int_{\pi/3}^{\pi} e^{j \omega n} d\omega \\
&= \frac{1}{2\pi} \left(  \frac{ e^{j \omega n} }{j n} |_{-\pi}^{-\pi/3} + \frac{ e^{j \omega n} }{j n} |_{\pi}^{\pi/3}  \right)\\
&= \frac{1}{2\pi} \left(  \frac{ e^{-j (\pi/3) n} -e^{-j\pi n}  }{j n} + \frac{ e^{j \pi n} - e^{j \pi/3 n} }{j n}  \right)\\
&= \frac{1}{n\pi} \left(  \frac{ e^{-j (\pi/3) n}- e^{j (\pi/3) n} }{2j} + \frac{ e^{j \pi n} - e^{-j\pi n} }{2j}  \right)\\
&= \frac{1}{n\pi} \left(  -\sin(\frac{\pi}{3}n) + \sin(\pi n) \right)\\
&= \frac{\sin(\pi n) }{n\pi} - \frac{ \sin(\frac{\pi}{3}n)}{n\pi}  \\
\end{align}
$$

at this point recognize  that
$$\frac{\sin(\pi n) }{n\pi} = \begin{cases} 1 &,\text{ for } n=0 \\0 &, \text{ for } n \neq 0 \\ \end{cases} = \delta[n]$$

Hence the I-DFT becomes:
$$ \boxed{ h[n] =  \delta[n] - \frac{\sin( \frac{\pi}{3}n) }{\pi n} }$$

Note that your $H(\omega)$ had a gain of $10$ instead of $1$ so multiply this with $10$.


How to extract traffic signs from a photograph? What image analysis techniques can I use to extract the traffic signs from an image such as the one below?

Edit:

After Anisotropic diffusion: The background that i don't want gets cleared a little

After Dilation :

Thresholding after Diffusion : Not able to figure out the best thresholding for this purpose

However i am not able to figure out how to remove the background?

Edit : i just want these parts of my image

Taking another input image :

Applying median filtering and edge detection :

After Bottom hat filtering:

How do i isolate the road signs please help?

Did you try something simple like correlation?

(EDIT). The idea behind correlation is to use a template (in your case a trained road sign sample), and compare it to every position in the test image. The comparison operation I've used to generate the images below is called normalized cross-correlation. Roughly speaking, you standardize (mean=0, standard deviation=1) the pixels in the template and the image part you want to match, multiply them pixel by pixel, and calculate the mean value of the products. This way you get a "match score", i.e. a measure of similarity between the template and the test image at every position in the test image. The position with the best match (highest correlation) is the most probably candidate for the position of the road sign.  (Actually, I've used the Mathematica function CorrelationDistance to generate the image below, which is 1 - (normalized correlation). So the darkest spot in the match image corresponds to the best match).

I don't have any other templates, so I simple cropped the sign from the second picture you posted:

Even though the template is rotated slightly, cross correlation still looks usable

and the best match is found at the right position:

(You'd need multiple scaled versions of each template to detect signs at any size, of course)


Inverse FFT - Synch the Phase Is there any way to synchronize phase of output of inverse DFT in each buffer? When I send the output of inverse DFT to the speaker it sounds nasty.

Of course I know the windowing functions but it just gives me in/decreasing level at the beginning and the end of buffer. It sounds better but with tremolo effect. I’ve also heard about overlapping windowing, I haven’t tried that but I am pretty sure it also won’t give me clean signal. Probably it will give me some phasing/chorusing effect that will appear periodically on each beginning and end of buffer. So I wonder how to repair that?

If you were attempting to modify the contents of a DFT before the inverse DFT, say for filtering, then you need to do overlap-add or overlap-save fast convolution instead of just an Inverse DFT.  You will also need to use a longer DFT and inverse DFT.  If you can determine the length of the impulse response of your filter, then the length of the new DFT/IDFT needs to be extended and zero-padded by the length of that impulse response, so that the filter convolution is equivalent to linear, not circular.

If filtering, the noise you hear is likely the result of doing a circular convolution, instead of zero-padding and overlapping.

If, instead you are synthesizing audio, look up phase vocoder methods.


need general resource about image matching i'm working on image matching mainly SIFT and SURF for my thesis. i read Moravec's and harris corner detectors and understand them quiet well in my own opinion.A COMBINED CORNER AND EDGE DETECTOR[1]

I even understoond harris-lapace detector but had no luck with understanding it's descriptor.Indexing based on scale invariant interest points[2]

i understood SIFT and SURF too.

what i didn't really understand is some of these topics : harris-laplace descriptor , how does hessian matrix help to detect good interest points , RANSAC(have general idea but need details for image matching) , Haar wavlet( also have a general idea, don't know how it's output look like and how i can work with that)

that's pretty much it.
i would really appreciate If anyone can help me with finding general reference for these topics.

(i have sent this question on stackoverflow. they suggest here would be better place to ask it and get proper answer)

Check out Szeliski's book:
http://szeliski.org/Book/drafts/SzeliskiBook_20100903_draft.pdf

There is also an old book on feature detection:
http://www.amazon.com/Feature-Extraction-Processing-Computer-Edition/dp/0123965497

You can read the sections that you care about. Also, I think it is always a good idea to read about scale space theory if you are to work with such algorithms. For example:

http://www.nada.kth.se/~tony/book.html


what is the use of signal processing in CT scans? I'd like to ask about the use of signal processing in CT scans ?
would be great if you suggest Book names or courses about it !
 thanks in advance !

The C stands for 'computer', which should give you a strong hint that the answer is yes. The fundamental problem is the inverting the inherent radon transform that comes out of the measurement setup. There are other signal processing related issues too, like denoising, or image processing ones like feature detection and image segmentation.


Intuitively, why is windowing function a low-pass filter? I'm trying to intuitively understand why a windowing function of a signal is a low-pass filter. I know its DTFT is a sinc function, making it mathematically a low-pass filter, but my logic says truncating a signal with a windowing function would preserve high frequencies better than low frequencies, since high frequencies have a better chance at showing up in the window.

Thanks.

Windowing is not a filter.

Windowing is a multiplication of two signals in time (the input samples with the window function: $x_w[n] = x[n] \times w[n]$ ).

What you get in the frequency domain is a (circular) convolution of the transforms of your signal and the window.

This convolution in the frequency domain can be seen as the spectrum being low-pass-filtered, giving a resulting spectrum with less detail.

[In the case of a filter, you get convolution in time and product in frequency.]


Does it matter if the scaling term is in the DFT versus the inverse DFT? From Wikipedia, the equation for the 1D DFT is

From a separate source, the equation for the 2D DFT is

Notice how the 2D DFT definition features a scaling term while the first definition does not. Is this just an artifact of where the scaling takes place between the DFT and inverse DFT? Does it matter which transform contains this scaling term?

It's mostly a matter of convention.
There are three options. $1/N$ for the inverse transform, $1/N$ for the forward transform and $1/\sqrt{N}$ for both. The first is the one most commonly used, the last one has the advantage that Parseval's theorem holds directly, i.e. $$\sum |x_n|^2 = \sum |X_k|^2$$


Cancelling out known signal with sampling offset I would like to cancel out a known audio signal at my receiver.

Specifically at my transmitter I am playing a song, and at my receiver I am recording it. Both transmit and receive sampling rate is 44.1kHz.

At the receiver I basically want to null out the song so that the only thing that is left is the ambient noise in the room.

When I lined up the received and transmitted signal, I noticed that there is sampling offset so I cannot simply subtract out the transmitted signal.

Are there tools in Matlab that are designed for this task?
And what are the general principles I can use to solve this?

What you are describing is basically an acoustic echo cancellator, as used in e.g. any car with a «hands free» phone solution. Use the (known) dry song track as a reference, and use a stochastic gradient search algorithm (eg NLMS) to online identify the FIR filter which relate your recorded signal to the known reference. As long as the channel (room acoustics) remains relatively stationary (LTI) and the ambient noise is uncorrelated to the reverberated music, you should be able to adopt standard methods to your application.


What is the DFT of a binomial filter? I hope this question is not too simple, I just started learning digital image processing. The 1D binomial filter of size 2 is defined by $B_2 = \frac{1}{4}\begin{bmatrix}1 & 2 & 1\end{bmatrix}$. The DFT should be $\frac{1}{2} + \frac{1}{2}\cos(\omega)$, but I have some trouble deriving this solution. I found it here.

My calculations:

$$X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N}nk} = \frac{1}{4} + \frac{2}{4}e^{-\frac{2\pi i}{3}1k} + \frac{1}{4}e^{-\frac{2\pi i}{3}2k}$$ Then I took the real part $Re(X_k)$, and found after simplifying $\cos^2(\frac{\pi k}{3})\cos(\frac{2\pi k}{3})$, which is incorrect. Only if I use two times the second element in $B_2$, I get the same result.

I'm also trying to understand the reason for applying the DFT in this case. My line of thought is: By taking the Laplace/z-transform of the impulse response, one obtains the system's transfer function. Because we are only interested in the steady-state response, we can ignore the imaginary part and calculate instead the fourier transform. Since the input $B_2$ is real, the real part of the DFT is enough. Is this correct?

I don't really understand the thing about cosines (as in: how is that helpful?) – a DFT is really just a mapping from a complex vector with $N$ elements to a complex vector with $N$ elements; and your calculation seems to be wrong, and I'm not sure where, but doing two of the three elements of the DFT manually might actually be enough to clear things up.

\begin{align}
X_0 &= \sum\limits_{n=0}^{N-1} x_n e^{-i2\pi\frac nN\cdot 0}\\
&= \sum\limits_{n=0}^{N-1} x_n e^0\\
&= \frac 14 (1+2+1)\\
&=1\\[2em]
X_1 &=\sum\limits_{n=0}^{N-1} x_n e^{-i2\pi\frac nN\cdot 1}\\
&=\sum\limits_{n=0}^{N-1} x_n e^{-i2\pi\frac n3}\\
&=\frac14 \left( 1 \cdot e^{-i2\pi\frac 03} + 2 e^{-i2\pi\frac 13} + 1 e^{-i2\pi\frac 23}\right)\\
&= \frac 14\left[1+ 2\cos\left(2\pi\frac13\right)+ i2\sin\left(2\pi\frac13\right)
+ \cos\left(2\pi\frac23\right)+ i\sin\left(2\pi\frac23\right)\right]\\
&\text{bit of basic trigonometry}\\
&=\frac14\left[
1+ 2\cos\left(2\pi\frac13\right)+ i2\sin\left(2\pi\frac13\right)
- \cos\left(2\pi\frac13\right) + i\sin\left(2\pi\frac13\right)
\right]\\
&=\frac14\left[
1+ 1\cos\left(2\pi\frac13\right)+ i3\sin\left(2\pi\frac13\right)
\right]\\
&=\frac14\left[
1-\frac12+ i3\frac{\sqrt3}2
\right]\\
&=\frac18+i\frac{3\sqrt3}8
\end{align}

Note how there's no taking of real parts or anything: the DFT of a sequence is complex in general, and you can't just drop the imaginary part. So, not quite sure why we're talking in control theory terms here, but no, you cannot just drop it. And: the discrete Fourier is not the same as the Laplace transform.


Understanding sine wave generation in Python with linspace I was trying to sample a 12.8 MHz sine wave (78.125 ns) signal at every 160us (micro seconds). Since 160us is multiple of base period 78.125ns(x2048) i expected to get a sample of fixed amplitude but instead what I am seeing is a another periodic sine wave. I don't understand why ?

I am doubting quantization error but shouldn't that generate uniform noise instead of a creating a periodic sine wave.

    import numpy as np
    from matplotlib import pyplot as plt

    fig2 = plt.figure()
    ax2 = fig2.add_subplot(1, 1, 1)

    capture_size1 = 2048
    timestep1 = 160e-6
    freq1 = 12.8e6
    time1 = np.linspace(0, capture_size1 * timestep1, capture_size1)
    w1 = np.sin(time1 * 2 * np.pi * freq1)
    ax2.plot(time1, w1, '.')

    plt.show()

Edit1 :
1. the 12.8 MHZ is intentionally under sampled

Adding the screenshot of the plot with capture_size1 = 2048, the sine wave has proper amplitude of [+1, -1]

Edit2: I tried to increase the precision by using Decimal and i see it is behaving as expected. I expect a straight line as the sampling point is an exact multiple of period.

from decimal import Decimal
from math import pi as mpi
from math import sin as msin
import numpy as np
from matplotlib import pyplot as plt

fig2 = plt.figure()
ax2 = fig2.add_subplot(1, 1, 1)

capture_size1 = 2048
timestep1 = 160e-6
freq1 = 12.8e6
time1 = np.linspace(0, capture_size1 * timestep1, capture_size1)
w1 = np.sin(time1 * 2 * np.pi * freq1)
ax2.plot(time1, w1, '.')

capture_size3 = Decimal(2048 * 16)
timestep3 = Decimal(160e-6)
freq3 = Decimal(12.8e6)
time3 = [Decimal(i) * timestep3 for i in range(capture_size1)]
w3 = [msin(Decimal(i) * timestep3 * Decimal(2) * Decimal(mpi) * freq3) for i in range(capture_size1)]
ax2.plot(time3, w3, '.')

plt.legend(["Actual", "Expected"])
plt.show()

Edit3:
I further did some analysis thanks to the comment by @jithin. Looks like this is an issue with linspace. I tried to generate the time interval by just multiplication as shown in below code and removed the original plot which used linspace(this is crucial) so now i am able to see the values in 1e-9 range as others suggest. So is there indeed issue with linspace ?

from decimal import Decimal
from math import pi as mpi
from math import sin as msin
import numpy as np
from matplotlib import pyplot as plt

fig2 = plt.figure()
ax2 = fig2.add_subplot(1, 1, 1)

capture_size1 = 2048
# timestep1 = 160e-6
# freq1 = 12.8e6
# time1 = np.linspace(0, capture_size1 * timestep1, capture_size1)
# w1 = np.sin(time1 * 2 * np.pi * freq1)
# ax2.plot(time1, w1, '.')

capture_size2 = 2048
timestep2 = 160e-6
freq2 = 12.8e6
time2 = [i * timestep2 for i in range(capture_size2)]
w2 = [np.sin(i * timestep2 * 2 * np.pi * freq2) for i in range(capture_size2)]
ax2.plot(time2, w2, '.')

capture_size3 = Decimal(2048)
timestep3 = Decimal(160e-6)
freq3 = Decimal(12.8e6)
time3 = [Decimal(i) * timestep3 for i in range(capture_size1)]
w3 = [msin(Decimal(i) * timestep3 * Decimal(2) * Decimal(mpi) * freq3) for i in range(capture_size1)]
ax2.plot(time3, w3, '.')

plt.legend(["multiply", "Decimal"], fontsize='xx-large')
plt.show()

The image of the above python code is below

Change the following line :

  time1 = np.linspace(0, capture_size1 * timestep1, capture_size1)

To the following:

  time1 = np.linspace(0, capture_size1 * timestep1, capture_size1, endpoint=false)

You will see correct results. Your original time instances is not what you intend because Python will create 2048 equally spaced point between 0 and 2048*Ts. What you want is equally spaced 2048 points starting from 0 and 160usec apart.

You can also use the following line if you dont want to use 'endpoint=false' :

  time1 = np.linspace(0, (capture_size1-1) * timestep1, capture_size1)

At your current sampling period of $160\mu sec$, it will mean that you are taking 1 sample of the sinusoidal every 2048 periods. There will be aliasing, but you are not bothered about that because you want to see fixed amplitude across discrete time. Basically, you want aliasing to happen.


Find the period of a signal with the DTFT plot I have an exercise and I'm struggling to resolve it. Here it is :

My problem is about the DTFT. I've always been taught that we use DTFT for infinite-lenght signal that are not periodic (if the signal were periodic I would have used DFS)... So, how am I supposed to find the period of x[n] and y[n] ? I must have missed or misunderstood something. Can somebody explain it to me ? (obviously without giving the explicit answer). Thanks !

The discrete-time Fourier transform (DTFT) can be used for general infinite length signals, and it can also be used for periodic signals if we allow Dirac delta impulses in the expression for the DTFT. You have the correspondence

$$\textrm{DTFT}\{e^{jn\omega_0}\}=2\pi\delta(\omega-\omega_0)\tag{1}$$

where it is understood that the expression $(1)$ is valid in the interval $[ -\pi,\pi)$, because the DTFT is always $2\pi$-periodic.

An $N$-periodic signal has Dirac impulses at multiples of $2\pi/N$. Note that not all of these Dirac impulses must be present, some of them can have a weight of zero. I'm sure you can take it from here.


rate 1/128 channel code for 5 bits of information I need to construct channel code for 5 bits of information for a wireless channel. Possible code rate could be as low as 1/128 (very reliable). At decoder I assume to use channel LLRs as decoder input to take into account channel conditions. At a first glance I'd use a repetition code (or spreading sequence) thus there will be 2^5 coded word in total and use ML decoder (correlation with each code word). These seems to be simple in realization. Are there any other types of channel codes for such requirements.

Since you say you want to take the channel conditions into account, one thing that immediately comes to mind is to suggest an adaptive link. Assuming you have access to the channel conditions, you can do two things to adjust the link:

Change the modulation: for example, think about using a higher order scheme like 16 QAM for "clear" channel conditions and using QPSK or BPSK for worse channel conditions.
Change the code rate of the channel code: this can be done by choosing a code with some rate, then depending on the channel you can perform puncturing to remove some parity bits (thus increase rate) when the channel is "clear".


Why is it a bad idea to filter by zeroing out FFT bins? It's very easy to filter a signal by performing an FFT on it, zeroing out some of the bins, and then performing an IFFT.  For instance:

t = linspace(0, 1, 256, endpoint=False)
x = sin(2 * pi * 3 * t) + cos(2 * pi * 100 * t)
X = fft(x)
X[64:192] = 0
y = ifft(X)

The high frequency component is completely removed by this "brickwall" FFT filter.

But I've heard this is not a good method to use.

Why is it generally a bad idea?
Are there circumstances in which it's an ok or good choice?

[as suggested by pichenettes]

Zeroing bins in the frequency domain is the same as multiplying by a rectangular window in the frequency domain.  Multiplying by a window in the frequency domain is the same as circular convolution by the transform of that window in the time domain.  The transform of a rectangular window is the Sinc function ($\sin(\omega t)/\omega t$).  Note that the Sinc function has lots of large ripples and ripples that extend the full width of time domain aperture.  If a time-domain filter that can output all those ripples (ringing) is a "bad idea", then so is zeroing bins.

These ripples will be largest for any spectral content that is "between bins" or non-integer-periodic in the FFT aperture width.  So if your original FFT input data is a window on any data that is somewhat non-periodic in that window (e.g. most non-synchronously sampled "real world" signals), then those particular artifacts will be produced by zero-ing bins.

Another way to look at it is that each FFT result bin represents a certain frequency of sine wave in the time domain.  Thus zeroing a bin will produce the same result as subtracting that sine wave, or, equivalently, adding a sine wave of an exact FFT bin center frequency but with the opposite phase.  Note that if the frequency of some content in the time domain is not purely integer periodic in the FFT width, then trying to cancel it by adding the inverse of an exactly integer periodic sine wave will produce, not silence, but something that looks more like a "beat" note (AM modulated sine wave of a different frequency).  Again, probably not what is wanted.

Conversely, if your original time domain signal is just a few pure unmodulated sinusoids that are all exactly integer periodic in the FFT aperture width, then zero-ing FFT bins will remove the designated ones without artifacts.


what is the difference between YUV format and HSV format? i wonder why HSV and YUV are different formats, even though they both works similarly (based on the brightness), if not what is the difference and i have confusions with YPbPr and YCbCr , can anyone give directions to understand these color formats.

So the best discussion I can point to is in this PDF on Charles Poynton's site.

The PDF goes on to much more detailed explanations of each color space, and how to convert between them. It's a bit longer than I want to type out (or cut-and-paste images of). :-)