What is done to minimize distortion due to the hold operation? A hold operation can be modeled using a step function over one sampling period i.e.

$R(t) = 1/T * (h(t) - h(t-T))$, $h(t)$ the step function

In frequency domain this is equivalent to

$R(jw) = e^{-jwT/2}*sinc(wT/2)$

We see there exists a phase distortion of $T/2$ and a magnitude distortion from the $ sinc(wT/2)$ term.

This means that all the reconstructed signal will be distorted unless operated at very low frequency.

Can someone explain what induces this distortion and what is done to minimize both the phase and magnitude distortion

What you are describing is the distortion introduced by an ideal digital-to-analog converter (DAC) in the analog domain.  Two things are typically done to reduce this distortion:

Analog filtering
Oversampling

As you note in your question, the distortion is modeled, in the frequency domain, by a sinc rolloff.  Increasing the sample rate before converting the digital signal to analog reduces $T$/2, which make the rolloff slower, or "flatter".

Analog filtering can range from simple low-pass filters to compensation filters that counteract the sinc rolloff.  These filters do not, I believe, fix the phase issues, but they can greatly improve the frequency amplitude response.

CD players are a good example of these techniques.  Audio CD's are sampled at 44.1 kHz, but before the DAC the CD players will typically increase the sample rate to a few hundred kHz.  I'm not sure how sophisticated their analog filtering is.  My guess would be that the cheap stuff does simple low-pass filtering and the expensive ones do sinc rolloff compensation.


Matching two video conversations I have two video where two person having conversation. They begin recording different times. I want to match these two videos' frame numbers, i.e. By matching audio.

How can I do that?

The problem you want to solve is called automatic audio alignment, or auto-synchronization of audio. Common algorithms are phase correlation, cross correlation on the raw audio waveform. People also use more "cooked" audio features like MFCCs, say each 10ms frame, and measure frame-wise similarity between the audio segments to be compared.

There are many threads on Stack overflow that show how to do this. For example: find-short-clip-of-audio-within-longer-clip-of-audio, cross-correlation-vs-fft-for-finding-phase-between-2-signals. Also some interesting repos on github, audiosync, videosync


Filtering difference of correlated measurements I have an application with two separated GPS receivers giving live positions and I'm deriving a heading/displacement from the vector between them.  Each set of position measurements is noisy and has dropouts but I do have a quality estimate. So far it's classic Kalman filter.

The relation between the two antennas isn't fixed, both are moving independently and generally one is moving at higher rate than the other.

But obviously some of the error is correlated between the two receivers, so do I filter each antennae location separately or the filter the displacement?
Or is there a cleverer alternative?

edit: As an example suppose you have a GPS on a boat and on a towed sonar array - and you are trying to track the relative position of the sonar receiver.

The general uncertainties on each antennae position are similar. But some errors will correlate between them - eg. if an atmospheric effect moves both readings 10m North this has zero effect on the difference but would effect any filter/average of each position.

Arguing against simply averaging/filtering on the displacement vector is that the motion of each receiver is different, the boat is relatively steady - the sonar bounces around in the waves. So the statistics of each if you have any sort of averaging filter are different.

I have a couple thoughts on your problem:

Regarding the handling of atmospherically-induced biases to the GPS receivers' reported positions, I don't think you need to do anything at all. If there is a truly common bias inherent in both positions, then it will simply cancel when you take the difference of the two positions. Thus, I would just treat each platform independently and subtract the time-aligned positions to yield your displacement vector estimate.

This is essentially the technique of differential GPS (DGPS), which was used successfully to combat the effects of the US government's policy of selective availability (SA) on the publically-available GPS carrier. SA was implemented by pseudorandomly dithering the timebase of the spacecraft; this would result in small position errors reported by GPS receivers.

This type of implementation meant that the SA error was highly correlated with the receiver position, such that two receivers in close proximity would experience approximately the same SA error. In the typical DGPS case, if one receiver's position is known precisely, one can solve for the SA error and resolve the second receiver's absolute position accurately, cancelling out the effects of SA. In 2000, the US government ended the use of selective availability, so this isn't much of an issue any more.

Your case is similar but slightly different, in that you don't know the absolute position of either receiver, but that doesn't seem to be what you care most about. Instead, you only care about the differential displacement between the two. Therefore, you can apply a similar technique; if there is a constant position bias inherent in both, then calculating the difference of the two positions will cancel it out, just as DGPS systems did for SA error.
It's not likely that you can improve the performance of the GPS receivers' reported positions by filtering them further. Dynamics tracking systems are a key part of any GPS receiver, and a lot of work has gone into their development. Especially if the receivers that you are using are of a recent vintage, I would hazard that any additional filtering that you would apply would be superfluous.

That's not to say that you would never want to process the reported positions further. There are some exceptions to my above statement:

If you have data from other instruments to fuse with the GPS readings, that can yield a gain. Sensors used for dead reckoning (e.g. gyroscopes or accelerometers) can yield measurements that can help to improve tracking of platform dynamics in periods where the GPS receiver outputs are unreliable or inaccurate. However, in time periods where the receiver measurements are assumed to be accurate, an optimal tracker would rely proportionally more upon the GPS measurements in lieu of the higher-variance dead reckoning sensors.
If you have some specific knowledge of your platform's dynamics that you know the GPS receiver isn't taking into account, then you may be able to obtain some gains. If you can specifically configure your receiver for maritime operation, you may be able to obtain better performance, as its tracker can make certain assumptions about your platform's position (e.g. zero altitude relative to sea level, no vertical velocity or acceleration). Failing that, it's possible that you could take the receivers' position outputs and apply some filtering to improve things, but your results will be very specific to the equipment that you use and whatever processing is applied by their internal trackers.


Design of Reed-Solomon Code where $n < q-1$ I'm trying to figure out how to design a Reed-Solomon generator polynomial where $n$ divides $q-1$ (I believe it's denoted by $n|q-1$) but where $n < q-1$. For example, lets say we have a $q=16$ and we want to design a codeword with length 3 and corrects 1 error. so we have the following:

$n = 3$

$l = \frac{q-1}{n} = \frac{15}{3} = 5$

$t=1$ (for correcting 1 error) therefore $\delta = 2t+1$

from here we can get our roots of unity as $\{\alpha^5, \alpha^{10}, \alpha^{15}=\alpha \}$. We then need to choose $\delta -1$ consecutive powers of $\alpha$, but i thought we chose these powers of $\alpha$ from the roots of unity, but there aren't any consecutive roots.

I can solve this for $n=15$ since I will have consecutive roots, however since $n=3$ in this case, i'm not sure how to proceed.

You have the right idea but are missing some of the fine details.  A cyclic Reed-Solomon code of length $n$, where $n$ is a proper divisor of $q-1$, is a BCH code and the key property is that if a sequence of $2t$ consecutive $\mathbf{n}$-th roots of unity are roots of the generator polynomial, then the minimum distance is at least $2t+1$. Here, $n = 3$ and $\alpha^5$, $\alpha^{10}$ and $\alpha^{15}=1$ are the three cube roots of unity. So, choose two of them to be roots of $g(x)$ (it doesn't matter which ones because you will always get two consecutive roots in this simple case!) and you are done. If you choose $\beta = \alpha^5$ and $\beta^2 = \alpha^{10}$ as the roots of $g(x)$, then $g(x) = x^2+x+1$, and so all the codewords on this $[3,1,d=3]$ Reed-Solomon code are of the form $\gamma +\gamma x + \gamma x^2, \gamma \in$ GF$(2^4)$, or equivalently $[\gamma, \gamma, \gamma]$,   that is, we have a repetition code!

It is also the case here that the cube roots of unity are elements of the subfield GF$(2^2)$ but this is not the case in general. Take $n=5$ which is also a divisor of $q-1=15$ in which case $\alpha^3, \alpha^6, \alpha^9, \alpha^{12}$ and $\alpha^{15}=1$ are the five $5$-th roots of unity which all belong to GF$(2^4)$.


Fourier transform of even/odd parts of a complex signal Why does Oppenheim state the following properties:

\begin{align}
\mathcal F\big\{x_e (t) \big\} &= \Re\big\{ X(j\omega) \big\}\\
\mathcal F\big\{x_o (t) \big\} &= j \Im\big\{ X(j\omega) \big\}
\end{align}

where $x_e (t)$ and $x_o (t)$ denote the even and odd parts of a signal $x(t)$, respectively, only for a real signal $x(t)$?
I don't understand how those two properties could break down if $x(t)$ were complex in general.

Defining
$$ x_e(t) = \frac{x(t) + x(-t)}{2} $$
$$ x_o(t) = \frac{x(t) - x(-t)}{2} $$
We can take the Fourier transform of each:
$$
FT\{x_e(t)\} = \int_{-\infty}^{+\infty} x_e(t) (\cos(2\pi\omega t) - j \sin(2\pi\omega t)) dt
$$
$$
FT\{x_o(t)\} = \int_{-\infty}^{+\infty} x_o(t) (\cos(2\pi\omega t) - j \sin(2\pi\omega t)) dt
$$

So if $x_e(t) : \mathbb{R} \mapsto \mathbb{R} $ then the $\sin$ term disappears and we are left with just a real-valued function.

So if $x_o(t) : \mathbb{R} \mapsto \mathbb{R} $ then the $\cos$ term disappears and we are left with a purely imaginary-valued function.

Now, if instead $x_e(t), x_o(t) : \mathbb{R} \mapsto \mathbb{C} $ then the term cancellations still happen, but we are still left with a complex value for each integral because both functions are now complex-valued.

As a result, the right-hand sides of your quoted equations do not hold in general.


Fourier transform of the sampled signal I want to calculate Fourier transform of the sampled signal in two ways. Let $$s(t) = \sum_{k = -\infty}^{\infty}\delta(t - kT)$$And $z(t) = x(t)s(t)$. So we have $$z(t) = \sum_{k = -\infty}^{\infty}x(kT)\delta(t - kT)$$Directly computing Fourier transform leads to $$\mathcal{F}(z(t)) = \sum_{k = -\infty}^{\infty}x(kT)e^{-j\omega kT} $$ On the other hand, if we use $\mathcal{F}(x(t)s(t)) = \frac{1}{2\pi}X(j\omega)\star S(j\omega)$ we have $$\mathcal{F}(z(t)) = \frac{1}{T}\sum_{k = -\infty}^{\infty}X(j(\omega - k\frac{2\pi}{T}))$$ How these are equivalent to each other mathematically? I mean how we can prove $$\frac{1}{T}\sum_{k = -\infty}^{\infty}X(j(\omega - k\frac{2\pi}{T})) = \sum_{k = -\infty}^{\infty}x(kT)e^{-j\omega kT}$$ And is there any intuition for this connection in terms of other subjects like sampling theorem or DTFT?

The equality
$$\frac{1}{T}\sum_{k = -\infty}^{\infty}X(j(\omega - k\frac{2\pi}{T})) = \sum_{k = -\infty}^{\infty}x(kT)e^{-j\omega kT}\tag{1}$$
is an instance of Poisson's sum formula. The term on the right-hand side of $(1)$ is just the Fourier series representation of the periodic function on the left-hand side of $(1)$.
The samples $x(kT)$ of the time domain signal are basically the Fourier series coefficients of the periodic spectrum of the sampled signal.
The significance of $(1)$ is that it shows that the discrete-time Fourier transform (DTFT) $X_d(e^{j\Omega})$ of the sequence $x_d[k]=x(kT)$ equals a periodized version of the Fourier transform of the corresponding continuous-time signal $x(t)$:
$$X_d(e^{j\Omega})=\sum_{k=-\infty}^{\infty}x_d[k]e^{-jk\Omega}=\frac{1}{T}\sum_{k=-\infty}^{\infty}X\left(\frac{j(\Omega-2\pi k)}{T}\right),\qquad\Omega=\omega T\tag{2}$$
Clearly, if $X(j\omega)$ is band-limited with a maximum frequency $\omega_c<\pi/T$, then the shifted spectra in the sum on the right-hand side of $(2)$ won't overlap, i.e., there is no aliasing and the signal $x(t)$ can be reconstructed without error from its samples $x(kT)$. That means that the basic form of the sampling theorem is implicit in Eq. $(2)$.


Moving average filter using DF-I and DF-II I understand that the moving average filter is simply the average of a number of points from the input signal to produce the output signal. If x[] is the input and y[] is the output and M is the number of points, then the filter is just $$y[i]=\frac{1}{M}\sum_{j=0}^{M-1}x\left[ i+j \right]$$
This can be done through many ways, but I want to ask about using Direct Form I and Direct Form II. How do I follow the DF-I and DF-II structures to create a moving average filter? Don't I need a two separate multipliers like $$b_{M}$$ and $$a_{N}$$ to use the DF-I and DF-II, which moving average filters don't necessarily have? How about if I want to use DF-I and DF-II to create an autoregressive filter? Isn't it even more incompatible with DF-I and DF-II?

To start, notice that the diagrams are in the $\mathcal{Z}$ domain, multiplying by $z^{-1}$ is the delay of one sample.
DF-I
Writing the equation for $y_n$ of the DF-I
$$y_n = \sum_{k=0}^{M} b_k x_{n-k} - \sum_{l=1}^{N} a_{l} y_{n-l}$$
Comparing with the equation of the moving average
$$y_n = \sum_{k=0}^{M} \frac{1}{M+1} x_{n-k}$$
we see that we could realise the moving average by settings $b_k=1/(M+1)$, and $N$ set to zero.
DF-II
For the DF-II realization you the sequences $a$ and $b$ must have the same size, so what we can do is to set $a_l=0$, and $b_k = 1/(M+1)$.
If you choose to compute a moving sum then $b_k=1$ and you can save the multiplications, the moving average still can be computed by a gain, thus, performing one multiplication per sample instead of $M+1$ multiplications per sample.


Finding system stability LTI system with impulse response:

$h(t) = e^t \sin(-3t) \text{u}(-t)$

I know the rule and formula , but I am lost in how to do this because of the absolute. Please help.

Sorry if I posted in the wrong section.

Thank you.

A system is stable if $ \int_{-\infty}^{\infty}  |h(t) |dt < \infty$

In your case $ \int_{-\infty}^{\infty}  |h(t) |dt = \int_{-\infty}^{\infty}  |e^t \sin(-3t)u(-t)| \ dt$

Since $\sin(-3t) = -\sin(3t)$,

$ = \int_{-\infty}^{\infty} | e^t \sin(3t) u(-t)| \ dt$

$ = \int_{-\infty}^{0}  e^t |\sin(3t)| \ dt$ , $u(-t) = 1$ if $ -\infty < t < 0$ else $0$

$ < \int_{-\infty}^{0}  e^t \ dt < \infty$

So the system is stable.


Pixel Binning: Effect on SNR for Hardware vs. Software Binning Is there a difference in SNR (or any other measure of image quality) between using 2x2 pixel binning in the camera and shrinking the full resolution image in software (cv::INTER_AREA) to the same size (w/2, h/2) ?

I imagine that SNR is higher with binning if pixel binning is done "in hardware" but I'm not sure how binning is actually implemented or if it does make a difference, to what extent the difference is measurable/negligible.

If you have access to the raw pixel data, you can implement binning in software too. They should be identical upto the step-size of the ADC, by which the quantized analog sum of pixels may differ from sum of quantized pixels.

Hardware binning has the advantage that it has to transfer less number of pixels, hence can operate at high fps compared to raw data transfer.


why if we pad any signal with enough zeros we can get the same result if we pad any signal with enough zeros we can get the same result as linear
convolution whilst using the fft function we compute circular convolution.
Why?

Hint:  try computing the linear convolution of $~[1\quad 1\quad  1]~$ and
$~[0\quad 1 \quad 1]~$ by hand and the circular convolution of the zero-padded
versions
$~[1\quad 1\quad  1 \quad 0\quad 0\quad 0]~$ and $~[0\quad 1 \quad 1\quad 0\quad 0\quad 0]~$ also by hand (or via FFTs if you are desperate enough).


Any meaning to the magnitude of a FFT when audio is recorded using a PC? I understand the theoretical meaning of the magnitude of a FFT from a mathematical sense. But what about when the sound audio is recorded using a PC microphone using a package like PyAudio, and FFT is calculated using Python. It seems the magnitude is unitless and the range depends on factors like sound card, number of bits used in encoding, etc.

Is it fair to say that the magnitude in this case has no meaning in an absolute sense but can be used to compare two signals?

Actually, to say the either the sound values or the DFT magnitutes are unitless (or meaningless) is quite misleading.  The sound values actually represent changes in the sound pressure level at the microphone.  This pressure has physical units.  You can think of the microphone and A/D conversion as a "transducer" which converts the units into something else.  Ideally, this conversion is linear.  So the units in the sound file do have meaning, but most of the time it is not worthwhile to calibrate them to the actual physical units.

The DFT magnitudes are then in the same units, after you normalize by 1/N, where N is the sample count of your DFT.  The bit depth has to do with precision, not scale.

Your conclusion is still correct, but now better founded.  You can compare two signals on the resultant scale.

Ced


y(t) a of an integrator circuit I have a signal $$ x(t)= \frac{1}{T} e^{-\frac{t}{T}} u(t) - \frac{1}{T} e^{\frac{t}{T}} u(-t) $$

and I know that it transits in a integrator circuit and I have to find y(t) in time and frequency domain. From theory on my book i know that

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

So I started to calculate this integral with my $x(t)$:

$$ \frac{1}{-RCT} \left[ \int_{0}^{+\infty} e^{\frac{-\tau}{T}} d\tau - \int_{—\infty}^{0} e^{\frac{\tau}{T}}d\tau \right] $$

but this gave me

$$ \frac{1}{-RCT} \left[  \dfrac{e^{-\infty } -1}{-\frac{1}{T}} - \dfrac{-e^{-\infty } +1}{\frac{1}{T}} \right] $$

and this gave me 0.

The correct result should be $$ -e^{ -\frac{|t|}{T} } $$

Thank you so much

The input and the integrator output are shown below:

So, you should consider two cases: $t<0$ and $t>0$.

In the first case ($t<0$), you have
$$\begin{align}
y(t) &= \int_{-\infty}^t x(\tau)d\tau = -\frac{1}{T}\int_{-\infty}^t e^{\frac{\tau}{T}}u(-\tau)d\tau \\
&= -\frac{1}{T}\int_{-\infty}^t e^{\frac{\tau}{T}}d\tau = -\frac{1}{T}Te^{\frac{\tau}{T}}\Big]_{-\infty}^t \\
&= -e^{\frac{\tau}{T}}\Big]_{-\infty}^{t} = -(e^{\frac{t}{T}} - 0) = -e^{\frac{t}{T}}
\end{align}$$
again, for $t<0$. So now you have
$$y(t) = \left\{\begin{array}{ll}
-e^{\frac{t}{T}}, & t < 0 \\
? \quad \:\:\:, & t > 0
\end{array}\right.$$
and you're looking for the other branch of the output. Can you show that it is $$y(t) = -e^{-\frac{t}{T}}, \quad t > 0?$$

If you do, then you'll get
$$y(t) = \left\{\begin{array}{ll}
-e^{\frac{t}{T}}, & t < 0 \\
-e^{-\frac{t}{T}},  & t > 0
\end{array}\right.= -e^{-\frac{|t|}{T}}$$ for all $t$.


Difference between order and length of a filter Let, the channel be modeled as a FIR model with order L. Mathematically, $y_n = \sum_{k=1}^L h_k u_{n-k} \tag{1}$ where $u$ is the data input to the channel and h is the channel coefficient. In statistics, the model y is known as the Moving Average process with order L. In signal processing, L is known as the length. Consider an example where there are 3 coefficients viz. $h_1,h_2,h_3$

Confusion 1: My confusion is, what is the order for this example? There are 3 coefficients  so the order should be 3. I am confused between these terms - order, length and delays.

Confusion 2: Can somebody please clarify what is the order for this example? Please correct me where I have made any mistakes.

In Matlab I have implemented the equation(1) as follows:

N = 1000;
u = rand(1,N);
for n = 3 : N
    y(n) = h(1)*u(n) + h(2)*u(n-1) + h(3)*u(n-2);
end

In your example order is $2$ and length is $3$.

Its has more to do with convention. Consider a polynomial in $x$:
$$ a_n x^n + a_{n-1} x^{n-1} + ... a_0 x^0 $$
The order of this polynomial is $n$, whereas, number of coefficients is $n+1$. In your case an FIR filter of order $L$, is a polynomial you write as:
$$ h_L \delta(n-L) + h_{L-1} \delta(n-L+1) + ... h_0 \delta(n)$$
The $z$-transform of which is:
$$ h_L z^{-L} + h_{L-1} z^{-L+1} + ...+ h_1 z^{-1} + h_0 z^{0}$$
Note this transform, which you might not have studied yet, just places the shift in $\delta(n)$ on power of $z$, e.g., $z^{0}$ and $\delta(n)$ correspond to no-delay, just scaling the input signal. $z^{-1}$ and $\delta(n-1)$ correspond to order $1$ and one delay.

I think real confusion comes from MATLAB, as it does not allow $0$ as first index value. You cannot write x(0) in MATLAB. The first index is $1$. Therefore, we have to start from x(1).


Fourier Transform Problem - absolute value, time-saving tricks, etc I am given the following signal:

$$[e^{-at}cos(w_{o}t)]u(t),\ a>0$$

Then I am told to find the Fourier Transform, which tells me I need an answer of the form:
$$X(jw)=\int_{-\infty}^\infty \! x(t)e^{-jwt} \, \mathrm{d}t.$$

I know I can reset the bounds of my integral with the unit step function, so my equation becomes
$$X(jw)=\int_{0}^\infty \! e^{-at}cos(w_ot)e^{-jwt} \, \mathrm{d}t.$$

From here, can I essentially solve this out and get a correct answer, keeping $w$ and $w_o$ as separate variables?

I also know I can solve it by using the relation
$$x(t)=a(t)b(t)\xrightarrow{\mathscr{F}} X(jw)=\frac{1}{2\pi}A(jw)*B(jw)$$

So, essentially, I can figure out a transform for each part and convolve to find my answer?

From the book examples, it seems $cos(w_ot)$ can be broken down into $\pi\delta(w-w_o)+\pi\delta(w+w_o).$ If we convolve this with the result from the transform of $e^{-at}u(t)$, a correct answer should be obtained.

It is known that $X(jw)$ when $x(t)=e^{-at}$ is $\frac{1}{(a+jw)}$

Therefore, since impulses sift through the other function in a convolution to get the nonzero values, is the following a correct result? $$\frac{1}{2\pi}(\pi)(\frac{1}{(a+j(w-w_o))}+\frac{1}{(a+j(w+w_o))})$$

$$=\frac{1}{2}(\frac{1}{(a+j(w-w_o))}+\frac{1}{(a+j(w+w_o))})$$

Thank you, sorry for the long question!

As you suggested, you could simply solve the integral using $\cos(\omega_0t)=(e^{j\omega_0t}+e^{-j\omega_0t})/2$. But as you've also noted, with this type of problems there is usually some smart way using known transform pairs. And what you suggested appears to me a very sane approach: convolve the known transform of $\cos(\omega_0t)$ with the known transform of $e^{-at}u(t)$ and you're done. I would recommend to you to cross-check your result by solving the integral, which in this case is also quite straight-forward. But I can assure you that your result looks good.


How to calculate $E_b/N_0$ for a 8-PSK? How to calculate $E_b/N_0$ for a 8PSK? My BER is given. But I don't understand how to calculate $E_b/N_0$ as it is inside a Q function.

Assuming that you have the correct formula for the BER (which is probably based on the assumption that a Gray code is used), then - as mentioned by MBaz in a comment - you would have to invert the $Q$ function, which is a problem for which no analytic solution exists. You can use numerical methods implemented in Matlab and probably in other software packages, but you can also just estimate the required value of $E_b/N_0$ from a graph of the bit error rate. I would assume that the latter is what you're supposed to do if that's a homework problem.

Simple approximations (actually bounds) of the $Q$ functions exist, but inverting them will very likely result in very inaccurate estimates of $E_b/N_0$.


Interpolation of estimated frequency between fft bins: what about phase? I am developing an additive re-synthesier, which involves identifying partials from FFT frames.

Partial identification requires amplitude peaks.  The simplest way to do this is to just look at the amplitude of each bin, but this of course can miss partials that happen to lay between bins, and there energy is spread over the two adjacent bins.

I understand there are various interpolation methods (of varying complexity) to estimate the actual amplitude peak and frequency... however, if taking such an approach, how to estimate the phase offset?  Is it even possible?

If you do an fftshift before the FFT (to center the phase reference in the data frame), you can use the same interpolation methods on the real (even) and imaginary (positive odd) components in the complex FFT result, and thus get interpolated phase estimates (referenced to the original data frame’s center) between FFT result bins.


Upsampling and Downsampling using IFFT and FFT (DFT Based Resampling) Assume in OFDM , there are N=64 subcarriers. I would like to upsample my signal by Factor 8 after IFFT and pass it through channel then downsample it by factor 8 before FFT. I can use upsample, downsample and interpolating functions. However I was told that there is an alternate to use (8*64) point IFFT and (8*64) points FFT. I used to think that since my symbol size is originally N then I can zeros at the end of my signal (64*8-64) that its length is 512 then I can apply IFFT and FFT correctly and take the first 64 samples and discard the rest. However after reading different post here (like this) enter link description here about adding zeros between samples or re-scaling the signal to get the correct result. I am wondering if anyone can explain to me about right way of upsampling and downsamping my signal using IFFT and FFT and also the reason for that. In different posts here people mention to add zeros between samples but they never said why.

Given a time domain signal $x_1[n]$ of length $N$ and its $N$-point DFT $X_1[k]$, If you pad $M$ zeros to the end of $x_1[n]$ and then take ($N+M$) point-DFT of it, this will give you an interpolated set of frequency samples, $X_2[k]$; i.e., a denser set of uniform samples of the DTFT, $X(e^{j\omega})$, than was previously available from $X_1[k]$.
Given a time domain signal $x_1[n]$ of length $N$ and its $N$-point DFT $X_1[k]$, If you carefully stuff (insert) $M$ zeros to the central portion of the DFT $X_1[k]$ (by moving those right half of the available $X_1[k]$ samples to the right by the amount of $M$) and then take ($N+M$) point-inverse DFT of the resulting $X_2[k]$, this will give you an interpolated set of time samples, $x_2[n]$; Note that you should scale the amplitudes in accordance with the interpoaltion amount.

Note that implementing 1-) is straightforward, while implementing 2-) requires careful selection of the correct indices at the central portion for the original length $N$ being even and odd, the details of which is provided by some examples here dsp.se,  including the link you have provided.


PSD - Spectrum analysis Could someone tell me where it is that I might be going wrong with the following problem?
Ok, so I'm computing a spectrogram, and, trying to return the PDS of the signal (to show where the signal has the most energies) which I will then use in plotting the decibels.
In my approach, I multiply each of the frames by the Hanning window, perform FFT of the resulting blocks and then only take the first half of the result (resulting vector).
Finally, I take each of the blocks (computed by the STFT) and perform PSD.
The problem that I'm having is this:
The chart below illustrates the results that I'm getting using my own C++ code:

Whereas in Python, using matplotlib if I replicate this, then I get the following:

I can't understand where the "noise" is coming from. I can only think that it could be some periodic noise which could be due to little discontinuities at the start and at the end of each frame. But I don't know where I have gone wrong.
The NFFT = 256 hop size = 128
Thanks

I'm a little confused by this statement given the previous one indicates you've already taken the FFT.

  Finally, I take each of the blocks (computed by the STFT) and perform PSD.

Generally the following is used; Seperate into frames (overlap if desired) -> apply window -> fft -> transform and scale into dB (if desired).

The scaling for a power spectrum or power spectral density can be a bit tricky at first, but this paper goes into some great detail about how to set it up.


Ways to analyze electrical signals without FFT? I am collecting a 15Hz sine wave at 1500Hz. Therefore, every 100 points of data collected should contain a single complete sine wave.  If I sample for 10 seconds, I will have 15000 points of data. If I perform an FFT on this data I will have very high frequency resolution.  However, If I wanted to divide these 15000 points into 100 point sections and average them together to get a mean signal of 100 points, and performed an FFT on that, my frequency resolution would be very low.  So, by averaging my data I am losing frequency resolution.

The reason I would like to get an average is because high amplitude noise is periodically injected into the sine wave.  I would like to throw away the sections containing noise and only have the frequency information of the average of clean data.

I'm fairly new to DSP so I'm not quite sure how to approach this problem.  Any thoughts would be greatly appreciated.

Please be specific about what kind of analysis you want to perform. My reply is based on the assumption that you want to measure the exact frequency of the input (since you are talking about frequency resolution) - that is known to be near 15 Hz.

Averaging is not the only option to remove the noise. You can filter your input signal by a low-pass filter with a cutoff-frequency near the frequency of interest. If your noise is very spiky (just a glitchy sample with an abnormally high amplitude), there are non-linear techniques which might be even more efficient at removing spurious transient events without much impact on the overall shape of the signal - for example, median filtering. The danger of averaging is that the averaged segments will have to be perfectly synchronized with the period of your input signal - if they are not, the phase differences will cause cancellations and your signal will lose amplitude as segments are averaged.

Even with averaging and/or in a situation with access to a limited number of samples, there are better techniques than the FFT. For example, accurately determining the frequency of a sinusoidal signal can be seen as a model fitting problem and can be solved by least square techniques. You could also use a PLL to try to lock into the input signal, and thus get the phase/frequency information - if this is what you are interested in knowing...


anti-alias filter of square pulses I'm trying to determine whether or not an anti-alias filter is needed for sampling square waves. The goal is to sample square wave pulses from a video detector with an ADC, do some time-domain digital processing on it, and reconstruct it with a DAC.

I do understand that signals with frequencies above the nyquist rate will alias into the "wrong" frequency bin. I guess the best analogy is the camera taking a picture of a car wheel turning at a certain rate, and at certain speeds it looks like the wheel stops or starts spinning backwards (a "misinterpretation" caused by the aliased frequencies). A square wave is made up of an infinite amount of odd harmonics but....

From an ADC perspective, it is just taking a sample of the voltage in time. I fail to see how a "misinterpretation" could be made since there is no "turning car wheel" to take pictures of at the wrong time. Do the harmonics alias in such a way that the wave shape is preserved?

In my mind, adding a filter to the signal will modify the shape of the original signal, getting rid of those upper harmonics. Depending on the filter design, it could add overshoot, ripple, and/or rise/fall time changes. So wouldn't the best representation of the pulse be obtained by direct sampling?

From an ADC perspective, it is just taking a sample of the voltage in time. I fail to see how a "misinterpretation" could be made since there is no "turning car wheel" to take pictures of at the wrong time. Do the harmonics alias in such a way that the wave shape is preserved?

You can reason this out yourself, in the time domain.  Consider a square wave with values -1 and 1; and sample it at exactly five times its period.  You'll get something like {+1, +1, -1, -1, -1, +1, +1 ...}.

If you just have to think about this in the frequency domain, you can, after a great deal of work, demonstrate to yourself that not only do the harmonics not alias in a way that preserve the shape, but in fact alias in a way that does not preserve the shape -- you end up with edges aligned with the ADC sample points.

  In my mind, adding a filter to the signal will modify the shape of the original signal, getting rid of those upper harmonics. Depending on the filter design, it could add overshoot, ripple, and/or rise/fall time changes. So wouldn't the best representation of the pulse be obtained by direct sampling?

That depends on what you're trying to do.  If your goal is to accurately capture the timing of the edges, to something less than the ADC sampling interval, then you need to round those edges out -- because a signal that goes {+1, +1, -1, -1} has less information about the location of the edge than one that goes, e.g. {+1, 0.75, -0.25, -1} or {+1, 0, -1, -1}.

And finally:

  I'm trying to determine whether or not an anti-alias filter is needed for sampling square waves.

That depends entirely on your problem at hand.  If you need every bit of information contained in the wave (and not just its timing), and if it truly has energy content out to infinity (which is physically impossible) then you need to sample infinitely fast.  If you know it's bandlimited, then you know the sampling rate you need to use.  If you know that, for example, the edges are super-sharp but what you really care about is their timing, then you can filter and acquire with an ADC, and infer the actual position of the edges from the way that the measured "square-ish" wave transitions.


How to recognize colors from the color chart? I am developing a vision application that is using the color chart below and camera to extract color of each patch in the chart. In order to do that, I have to first detect the chart area in an image and match the area with existing template which contains locations of patches.

My question: I need to detect the color patch area by recognizing four corners of the chart. You can see the corners that the small inset rectangles indicate in the image below. I know one way to go about it is to let user to select those corners by clicking it. But is there any way to automatically detect four corners from the image.

You can house hough lines to detect the color area

First detect all lines
Lines with longest length would be the outer most ones
Pick the top horizontal line and do template match for the squares.
Once you find the squares, traverse downwards toll you reach the other set of squares.
Mark all points where the squares were found
You now have an enclosing area for the color chart.


Speed of car with sound I've got this task:

A traffic policeman standing on the highway notices a change in tone of exactly one major third (f2/f1 = 4/5).
What speed can he infer? Speed of sound c = 340 m/s

I've got this attempt: 340 m/s * (1+√5/4) = 765m/S = 2754 km/h
Can someone tell me if this is correct or if I'm completely wrong?

Hints

Define the frequency of the car as $f_0$, the heard frequency when it's approaching as $f1$ and the heard frequency when it's departing as $f_2$
Write the Doppler formula for both $f_1$ and $f_2$ using a receiver velocity of 0 (since the policeman is stationary with respect to the air) and a car velocity of $v_s$. Pay attention to the sign of the source velocity.
Write the ratio of both formulas and equal it to a major third. $f_0$ should drop out and you are left with a single equation for $v_s$
Solve for $v_s$


Is there a simpler way to calculate the amplitude response of the following filter I would like to calculate the amplitude response $|H(z)|$, $z=e^{j\omega}$, of the following filter:
$$H(z)=\frac{\frac{b}{2}+z^{-2}}{2+bz^{-2}}$$
and I would like to avoid using Euler's formula and complex numbers cartesian form as much as possible.

The given transfer function has a specific form. With the denominator polynomial $A(z)=2+bz^{-2}$, you can write the transfer function as

$$H(z)=\frac{z^{-2}}{2}\frac{A\left(\frac{1}{z}\right)}{A(z)}\tag{1}$$

Since for $z=e^{j\omega}$ (i.e., on the unit circle, where we evaluate the frequency response), we have $|A(z)|=|A(z^{-1})|$, the magnitude of $H(e^{j\omega})$ is given by

$$\left|H(e^{j\omega})\right|=\frac12\tag{2}$$

The given filter is a (scaled) all-pass filter. It's easy to recognize all-pass filters by noting that their numerator and denominator polynomials have mirrored coefficients.


Audio Watermarking - What is typically embedded? Just starting research in audio watermarking and have read many papers regarding audio watermarking but I am still a little confused on what exactly gets embedded into the signal. For example, one paper I read embedded an image or written signature into the audio.

My question is what typically gets embedded in the signal using audio watermarking? Is it just random data, a picture, someone's name? Thank you!

Research papers or textbooks won't delve into the specifics - they just assume you have a stream of binary data that needs to be reliably and undetectably embedded into the audio file.

What this bitstream is actually going to be is application dependent. It could be:

Digital rights management applications: A short serial number identifying the individual transaction in which the audio file has been sold. Given the small size (a hundred bits) of this information, it is typically repeated over and over the file to maximize the probability of recovering it.
Metadata management: Metadata such as artist name, song title, or even lyrics.
Proper steganography: Any secret data you want to share with the recipient of the audio file - without an eavesdroppers suspecting that you are exchanging anything other than music. This can be a text-file, an image - any binary data.
Informed source separation: Any information about the individual sources the audio file is made of, that would subsequently help in recovering these individual sources.
Second-screen interactive TV applications: These application allow a device such as a smartphone or a toy to display secondary information synchronized with a TV show. To allow the synchronization, the audio signal of the show contains a regular watermark containing a show identifier and time stamp.


Fitting a grid to points in a 2d space I have this problem where I have multiple points in a 2D space. A subset of these points makes up a rectangular grid. Some points of the grid might be missing. The shape and number of points the grid consists of are known. Also, the other points not making up the grid are few. The general location and rotation of the grid are also known. But we don't know the exact location or rotation. What would be a good, and most important, efficient way to find this grid?

A perfect grid:

The grid can appear as shown below with a missing point:

Or with for example a rotation and being split in two as shown below:

In the actual problem, the grid is larger but always rectangular, also we can assume the gird is located around the center of the data.

The points come from an image where there is monospaced text. The missing points and noise are due to glare in this image. The points are the centers of the
(by preprocessing) found characters (so not exact). The text in the image always consists of the same number of characters and lines. The idea is reconstructing this grid to be able to find the lost characters. We do not know what points are characters and what points are not, but I am aware this could be approximated by analyzing the "blobs" we find in the image.

You could represent the grid with 5 parameters:

(x0, y0) to represent the offset from your image origin to the grid
origin
(xs, ys) to represent a multiplication factor from your image
plane to your grid plane. If your grid is spaced apart every 100
pixels on your image plane, these values might be 100,100
theta to
represent the rotation of the grid

You could fit the grid to point in image by varying these 5 parameters and optimize these parameters by calculating a cost function at each step.

An example cost function would be the total Euclidean distance between known point and grid points. The grid parameters that best overlap your given point should represent an optimal solution in most cases, especially in cases where you have several "correct points" and few of the "outlier points".

For efficiency you could ensure that your x0 value is [0,xs] and your y0 is [0,ys]. Since it sounds like you have a guess for the size of your grid, you could input a decent guess for xs and ys and also constrain these values from changing too much.


Isolating Bat Species without prior training base I am attempting a project and I have been given Bat calls (in .wav) which are not ultra-sound navigation calls and I'm looking to isolate the different species of bats and monitor their activity.

The problem I am having is that I do not have a prior training set to be able to feature matching on the samples to try and identify them. I have been doing part research and I have looked at the following Link which gives details of the species of bats, their rep rate, FM/CF and Peak power.

From the link above, I am assuming that I can have some kind of array which contains the species of the bats, as well as their frequency (range) or "Peak Power" and then using an FFT algorithm, convert the samples that I have obtained and then match the output of the FFT against the pre-set array and therefore enable a platform to monitor their activity.

Here are the 2 questions:

1) Based on the above, does would this be an ideal solution in order to isolate the bats and therefore monitor their activity or would I need to implement further algorithms in order to isolate the bats calls?

2) I have only ever had experience of using the Cooley-Tukey FFT algorithm, would this be sufficient enough for this type of project?

Thank you for taking time to read this.

1) Based on the above, does would this be an ideal solution in order
  to isolate the bats and therefore monitor their activity or would I
  need to implement further algorithms in order to isolate the bats
  calls?

The short answer is no, this would be no where near ideal. This type of work strongly overlaps with speech/voice recognition, and the algorithms (or more so, the features) involved are not trivial.

There are certainly frequency features that one can extract from audio signals of species to try and cluster them, and the DFT would certainly be the transform of choice for this type of raw feature. However the DFT will only give you frequency scores, although it where one would start. Other frequency metrics and features can be more complex. (Look at cepstral analysis, MFCC, etc).

  2) I have only ever had experience of using the Cooley-Tukey FFT
  algorithm, would this be sufficient enough for this type of project?

Remember that the DFT is a Transform, whilst the FFT is an efficient algorithm for this transform. Furthermore, the Cooley-Tukey Algorithm is one way of implementing an FFT. In your work, you will use canned libraries that have the FFT. You don't need to implement an FFT in the 21st century, just as you would not be expected to implement the square root function.

I strongly suspect that you will need to compute the Short-Time-Fourier-Transform (STFT) for this type of project. In anticipation of your next needs, please see the wiki link on the STFT, as well as this post here. The STFT allows you to analyse signals in the dual time-frequency domain. You will most likely want to start your feature extraction in that space.

  let's assume I have computed the STFT and the MFCC and extracted the
  features, how would I therefore identify the species of the bats?
  Would I therefore need some kind of training in order to complete this
  task?

Once you have decided on what sets of features you need to use, you will of course need a learning algorithm to discriminate them properly. Since you do not have a training set, there are two options. The first obvious is to attain a training set, and/or meticulously label your data based on what you know to be true. Once you have your labels, then you can use supervised training techniques like logistic regression or a simple perceptron to come up with separating hyperplanes in feature space that disjoint your data. The second option however, would be for you to use unsupervised learning techniques (such as K-means for example).

If you use something like K-means, you will need to know apriori how many different bat species you are trying to discriminate against. This is the value of your hyper-parameter $k$. The output of K-means will be $k$ vectors, each of them being in $N$-dimensional points in feature space, representing your cluster centroids. Now, every time you get a new audio sample, you:

Compute $N$ features and place this into a feature vector.
Compute the L2 vector norm, (aka, the Euclidean distance) between the given feature vector, and the $k$ centroid vectors you pre-computed from using $k$-means.
Pick the $i$th centroid out of the $k$ cetroids that corresponds to the closest distance to your point.

Mathematically, you want to declare a bat of type $i$, among $k$ type bats by:

$$
{arg\,min}_{i \in (1, 2, ...k) }  \sqrt{\sum_{n = 1}^{N} (\mathbf{x} - \mathbf{f_i})^2 }
$$

where $\mathbf{x}$ is an $N$ x $1$ feature vector of your unclassified audio sample, and $\mathbf{f_i}$ is the $i$th $N$ x $1$ cluster centroid found by $k$-means.


How to realize Poles and zeros at infinity??especially through transfer function? I have a question regarding the poles and zeros at infinity
I often read here in DSP SE and also in some textbooks about poles and zeros at infinity

This question also answers somehow (but not in much easy/simple way)
 Pole/Zero existence at infinity

I am confused about how to realize the poles and zeros at infinity especially from observing/looking at system transfer function??

It's actually quite straightforward: positive powers of $s$ (or, in discrete-time, $z$), correspond to poles at infinity. Negative powers give you zeros at infinity.

Let's look at some examples. In continuous time, an ideal differentiator has the transfer function

$$H(s)=s\tag{1}$$

Clearly, $\lim_{s\to\infty}H(s)=\infty$, hence you have a pole at infinity (and a zero at $s=0$).

The transfer function of an ideal integrator is

$$H(s)=\frac{1}{s}\tag{2}$$

Here we have $\lim_{s\to\infty}H(s)=0$, i.e., a zero at infinity (and a pole at $s=0$).

In discrete-time, a delay of one sample corresponds to

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

and you get $\lim_{z\to\infty}H(z)=0$.

A (non-causal) advance by one sample is represented by the transfer function

$$H(z)=z\tag{4}$$

In this case you have a pole at infinity: $\lim_{z\to\infty}H(z)=\infty$

Note that in continuous time as well as in discrete time, systems with poles at infinity cannot be causal and stable.


Difference between Gabor filtering and Discrete Wavelet Transform Both Gabor filtering and discrete wavelet transform (DWT) analyze the image in both spatial and frequency domains, unlike Fourier transform which analyzes the image only in the frequency domain. What is the difference between DWT and Gabor filtering?

Per se, a  Gabor filter  in image processing is one linear filter at a certain scale and 2D frequency used for orientation filtering and texture analysis.
It would be easier to compare Gabor representations and discrete wavelet transforms. Both are related to a linear decomposition of possibility multidimensional data at different scales with somehow oscillating functions. Main differences are:

Discrete wavelets: critical or non-redundant scheme (stable, invertible), discretize some families of wavelets, more scale-based  though wavelets range from weakly oscillating (Haar) to oscillation (Shannon), often applied separably across dimensions so weakly directional, a lot of statistical and approximation properties (moments, regularity). Their redundant counterparts: discretized continuous wavelets, shift-invariant or stationary wavelets.

Gabor transforms: (highly) redundant decomposition, mostly one shape: a modulated Gaussian, more frequency-based though computed at a couple of scale (Gaussian spread), inherently non-separable or directional. Their non-redundant counterparts: modulated or orthogonal lapped transforms, Malvar wavelets.


How to Shift Carrier Frequency in QAM Signal? If I wanted to apply a 100Hz frequency shift to a QAM signal with a 500MHz carrier frequency, I would

Demodulate it at the carrier frequency
Apply a frequency shift to the original signal
Modulate it again at the new carrier frequency

Is there a way to apply a frequency shift to the IQ values of a QAM signal without demodulating it?

Perhaps I reply too late but never the less ..

I am not sure how this problem is actually intended to be understood and the info on electronics.SE does not help much. Of questions some could be:

Why demodulate at the carrier? Should it be understood in the way that you want to convert the information at the carrier to baseband in-phase ($I$) and quadrature-phase ($Q$) signals?
What do you mean by the original signal? Is that an RF signal or the $I$/$Q$ signal components?
Why do you want to modulate (I assume you just mean frequency translate) the signal to the new carrier frequency?
What is actually the objective of the question?

Without this information I am forced to guess a bit.

As I understand it you have a QAM (for the following, the modulation is not relevant) signal at a carrier frequency $f_c=$500$\:$MHz which can be expressed as:

$$
x(t) = I(t)\cdot \cos(2\pi f_c t) - Q(t)\cdot \sin(2\pi f_c t)
$$

where the in-phase signal $I(t)$ occupies the frequency band $-B \leq f \leq B$ and the quadrature-phase signal $Q(t)$ also occupies the frequency band $-B \leq f \leq B$.

To simplify things I have omitted an initial constant phase in $\cos()$/$\sin()$ - this can easily be included though and in fact any constant phase can be included in $I$ and $Q$ without loss of generality.

To me it seems as if you want to force an offset by saying that the carrier is $f_c=f_0+f_\Delta$ and then you would like to know how "new" $I$/$Q$ signals referred to just the carrier frequency $f_0$ looks like. This means that we can rewrite the original signal as:

$$
x(t) = I(t)\cdot \cos(2\pi (f_0+f_\Delta) t) - Q(t)\cdot \sin(2\pi (f_0+f_\Delta) t)
$$

which we need to rewrite as:

$$
x(t) = I_\mathrm{new}(t)\cdot \cos(2\pi f_0 t) - Q_\mathrm{new}(t)\cdot \sin(2\pi f_0 t)
$$

This modification be done by applying standard trigonometric identities (see here).  What we need is:

$$
\cos(\alpha + \beta) = \cos(\alpha)\cdot\cos(\beta) - \sin(\alpha)\cdot\sin(\beta)
$$

$$
\sin(\alpha + \beta) = \sin(\alpha)\cdot\cos(\beta) + \cos(\alpha)\cdot\sin(\beta)
$$

Using these identities we can then determine the "new" in-phase and quadrature-phase signal components as:

$$
I_\mathrm{new}(t) = I(t) \cdot \cos(2\pi f_\Delta t) - Q(t) \cdot \sin(2\pi f_\Delta t)
$$

$$
Q_\mathrm{new}(t) = I(t) \cdot \sin(2\pi f_\Delta t) + Q(t) \cdot \cos(2\pi f_\Delta t)
$$

The "new" in-phase and quadrature-phase baseband signals occupy the frequency band $-B + \vert f_\Delta \vert \leq f \leq B + \vert f_\Delta \vert$.  The bandwidth increases no matter if we change $f_\Delta$ in the positive or negative direction.  Should it be needed the opposite transformation can also be done such that move from known "new" signals $I_\mathrm{new}$/$Q_\mathrm{new}$ to the $I$/$Q$ signals.

By use of this approach you can receive the signal at the carrier frequency $f_c+f_\Delta = 500\:$MHz, perform quadrature down-conversion as normal, and the apply the above technique to modify the effective carrier frequency by using e.g. $f_\Delta = 100\:$Hz or whatever.  As shown above you do not need to demodulate the signal and the conversion can be applied to any signal and be handled at baseband.  Depending on the actual application, which is not known from the question, it may not be allowed to use the above. You can extract the original $I$/$Q$ signals directly from $x(t)$ if this signal is quadrature down-converted with the carrier frequency $f_c$.  If you quadrature down-convert $x(t)$ with a locally generated frequency of $f_0$ you extract $I_\mathrm{new}(t)$/$Q_\mathrm{new}(t)$.  If you have $I_\mathrm{new}(t)$/$Q_\mathrm{new}(t)$ you can reconstruct $I(t)$/$Q/t)$ from the above equations (solving two equations with two unknowns).

One thing to remember, though, is that the requirement to sampling frequency increases as the bandwidth of the "new" signal is wider than the original signal.  In case of a $100\:$Hz offset this is likely not a big issue but it obviously depends on the bandwidth $B$ of the original baseband signal.

If the problem is more in the direction of RX a signal at one carrier and relaying to TX the same signal at another frequency then some other technique should be applied.  An obvious technique is to quadrature down-convert and then just quadrature up-convert to the desired offset carrier.  Also a translational loop combined with amplitude information can be applied. The translational loop filters the phase part of the signal.  If further filtering is needed depends entirely on the required spec.


Question about poles and zeros in AR filter For AR HP filter

  zeros, to the right of the imaginary axis
  poles outside the unit-circle
  zeros on the real axis
  poles, to the left of the imaginary axis

Apparently, the right answer is (4). Why?

We did learn in class some properties of the zeros and poles, and I probably need to relate it to the fact that it's a AR HP filter.

I'd be glad for your help.

Thanks!

EDIT:
I know that AR is an all-poles filter so it must be (1) or (4). What is the meaning of the position of the pole with relation to the imaginary axis?

EDIT 2:
Is it because the large $\omega$'s (frequencies as angles of the unit circle) are at the left side of the imaginary axis and so we get an HP filter?

An AR process is defined in a way that it is inherently causal. It depends on past and present values, but not on the future. Thus, its ROC must be:

In the continuous case, $s$ such that $|s|>\sigma$ for some $\sigma\in\mathbb{R}$.
In the discrete case, $z$ such that $|z|>r$ for some $r\in\mathbb{R}$.

For an AR filter, all the zeros are at the origin. This could make the 3rd option you wrote correct, as the origin is indeed on the real line. The first one is definitely incorrect.

Regarding the poles, it depends on the coefficient of the process, $a$. If $|a|<1$, the filter will be stable and all the poles will be on the left half plane of the $s$-plane, or inside the unit circle in the $z$-domain. If $|a|>1$, the contrary will happen.

So the answer would be the 4th, if the filter is required to be stable (which makes sense). If the filter is unstable, then the 2nd one would be the one to chose.


Estimating the time offset (delay) between two audio signals in real-time I'm looking for an algorithm which can accurately estimate the time-offset/delay between two audio signals in real-time. One audio signal is a loudspeaker signal and the other signal is the microphone signal.

Audio is processed in blocks of 10ms (sample frequency 16kHz).

There is light office noise embedded in the signals
The algorithm must be able to estimate delays up to 500ms.
The delay is usually constant for a period of time (with some variability -10ms to 10ms), but occasionally the delay "jumps". Those sudden jumps can be in the order of 50 to 80 ms. The algorithm must be able to handle that.

I've read that GCC PHAT is one algorithm which can be used, but I'm not sure how it would work in a real-time environment.

There are a variety of techniques available. See this thesis for a nice write-up.

See the list of algorithms:

In particular, note that they recommend pre-whitening the signals before correlation --- to ensure that the cross-correlation has a nice, narrow peak for better "localization".


Resample with unknown rate change I am doing digital signal processing for a software defined radio. I have a signal coming in packets of data (several samples) at an average sample rate of approximately some known value FSin. I am given a request for packets to be sent out periodically, with an average rate out of approximately FSout. If I knew these rates exactly I could resample my data by a ratio R=FSout/FSin and use a FIFO to make sure data always has somewhere to go when I get a packet and there's always data to send when I'm asked to send a packet. However, the clocks of the radio and my computer are not perfectly synchronized, so there is no exact known ratio for R. What would your solution be?

This is what I am thinking: I will have to adjust the ratio by which I resample over time to keep the size of the FIFO some target, on average. If my ratio goes too high, my buffer begins to grow and the latency gets large or it overflows, if I am using a circular buffer. If my ratio goes too high, my buffer begins to shrink and eventually underflows. I will periodically find the error between the number of elements in my buffer and the target size. The integral of this error over time, I, will be used to determine R by R=e^-I. The integrator will be set with its initial value I=-ln(R). If the buffer size drifts too high, I goes towards zero and R drops, shrinking the buffer. If it drifts too low, the opposite occurs. What do you think? Should I expect anything weird (oscillation through my feedback path, etc.)?

I have done this before and it works as long as you tune your parameters correctly. There should be a gain factor on the integrator which is < 1.0.

If your application is audio, one thing to watch for is audible pitch changes during periods when the ratio is changing. You can add a proportional term to your loop filter to tune the transient response if this is an issue.


How to add a frequency offset to a buffer of MPSK signal for simulation purposes? I need to simulate some channel effect to a buffer of MPSK signal. How to add a frequency offset to a buffer of baseband MPSK signal?

Multiply the time domain samples by a time-varying complex exponential to effect a (cyclic) shift in the frequency domain.
$$y[n] = x[n]\cdot e^{j\pi \frac{f_{offset}}{F_s/2}n}$$
Where $f_{offset} \in \left(-\dfrac{F_s}{2}, \dfrac{F_s}{2}\right)$.
This is a result of the Shift Theorem for the DFT.


How to make series circle shift invariant I have a series of numbers. They may be integer, floating, or binary. Length of series is constant. I need criteria which is circle shift invariant. For more information assume series S1, S2, S3, S4 and S5:

S1 = {A, B, C, D, E, F}
S2 = {F, A, B, C, D, E}
S3 = {D, E, F, A, B, C}
S4 = {A, C, D, F, E, B}
S5 = {A, B, C, D, E, G}

where A, B, C, E, F, and G are numbers. S2 and S3 is resulted by shift circle of S1. So the criteria should return same values for these three series. S4 have same number of S1, S2, and S3 but order of numbers are not the same. S5 and S1 are similar so I want to criteria of S1 and S5 be similar.

I want to use them as descriptors in image processing. In other words, I want to make a feature (like LBP, or HOG) rotation invariant. For example mean and variance can be used but I need more. Are there any suggestion?

A possible way to approach the problem would be to extract the maximum values of the discrete circular convolution of each pre-normalized sequence with an NXN matrix formed by the elements of an orthonormal base of the Nx1 sequence space.

The output of this step would be a Nx1 vector that would contain a measure of the similitude of each sequence with the each base sequence. If base sequences are adequately chosen, similar sequences should have similar output vectors, and an operator like mean square difference between outputs should tell how similar are two input sequences.


How do you normalize a sample-by-sample processed signal? I am writing an audio effect based on a second-order filter that adds two samples. The resulting output is clipped. So I multiplied the signal by 0.5, which results in an output that is too low. Since I am process sample-by-sample, I don’t see how you can apply normalization because the peak-values in the normalization window will vary too much.
How do you do this? Calculate the right audible half volumes in dB and add those?

Calculate the right audible half volumes in dB and add those?

Can't do that: dB is a logarithmic thing. Adding logarithms means multiplying the values. You want to build a linear filter, not a self-intermodulation generator, don't you?
If your filter coefficients are really just [1 1], well, then multiplying by 0.5 exactly retains your original signal power for low-frequency signals. Since this is a low-pass filter, of course higher frequencies than DC are attenuated, which is not a fault of the 0.5, but exactly the result of your filtering.
A usual way to avoid clipping that only happens rarely  would be a soft-clip function that lets through any "benign" amplitudes and starts saturating "softly" towards the maximum amplitudes.
If you want to do something more "time-smooth", that would for example be an adaptive gain control (AGC), which reduces gain gradually as your signal works towards higher peaks. However, that's not a "sample-by-sample" thing, but something with at the very least memory (and often even implemented in a somewhat acausal manner).


How do you calculate RMSE in degrees for a DOA estimator system with multiple incoming signals? Let's say you have a system with 1 transmitter and 5 receivers. We'll use the MUSIC algorithm to determine DOA of incoming signals.

If we have 1 signal (at 0° for example), it's easy to calculate RMSE vs SNR. Simply take the peak of the MUSIC spectrum (an angle) and find the difference between the actual target angle.

However, what if we have 2 signals? One at 0° and a second at 10° for example. You can no longer look at the MUSIC spectrum and simply pick the 2 greatest points, since the absolute greatest points may be at 0.0° and 0.1°, even if there's a clear (lower) peak at 10°.

What would be the/an appropriate way to measure the error when multiple signals are present? Or perhaps another metric (besides RMSE) is needed, such as probability of observing 2 peaks, etc.

The way to do it is to find peaks and not maximum values. Like the function findpeaks in MATLAB (local peaks) you can define resolution (distance between peaks and more...)


Constant Magnitude Chirp across audio band (Parseval's Theorem?) I want to excite a loudspeaker with a constant magnitude chirp that sweeps across all frequencies in the audio band. I downloaded a test file from this Audio Check - High Definition Audio Test Files to start with.
I plot the log magnitude of the file with the following code which I use for other applications as well.
[X Fs nbits] = wavread('audiocheck.net_hdchirp_96k_-3dBFS_log.wav');
X = X(:,1);
fft_prep = fftshift(fft(X));
fft_mag = real(fft_prep).^2 + imag(fft_prep).^2;
pos_fft=fft_mag(ceil(length(fft_mag)/2+1:length(fft_mag)));
db_fft=20*log10(pos_fft);
figure(1);
plot(linspace(0,48000,length(db_fft)), db_fft);

I see this graph:

I was hoping to see a constant magnitude response. A friend of mine mentioned something about how Parseval's Theorem states that there will be different energy in different frequencies. I found this link that talks about Parseval's Theorem and also says in the section "Spectral features and energy of chirps":

Important and desired spectral properties of excitation signals to improve the quality of
wideband measurement are: flat amplitude spectrum with minimal fluctuation (ripple) together with the absence of overshoots inside the generated (excitation) bandwidth Bexc =ffin - fst; steep drop-down of the amplitude spectrum outside the bandwidth Bexc; maximal energy-efficiency, i.e., the ratio between the energy lying within the generated
(excitation) bandwidth Bexc and total energy of the signal.

To avoid worrying about a chirp someone else made. I am also trying to generate a chirp in Matlab and I am seeing the following response with:
t = 0:0.001:10;      % 10 seconds @ 1kHz sample rate
fo = 10; f1 = 400;   % Start at 10Hz, go up to 400Hz
X = chirp(t,fo,10,f1,'logarithmic');
fft_prep = fftshift(fft(X));
fft_mag = abs(fft_prep);
pos_fft = fft_mag(ceil(length(fft_mag)/2)+1:length(fft_mag));
db_fft = 20*log10(pos_fft);
figure(1);
plot(linspace(fo,f1,length(db_fft)), db_fft);

How can I generate a relatively constant magnitude sweep across the audio range in Matlab like I see in figure 8 above? Any advice is welcome and I appreciate the help.

If you change the "X = chirp(t,fo,10,f1,'logarithmic');" to "X = chirp(t,fo,10,f1,'linear');" you will get the response that you want.

To have the energy be equal it must spend an equal amount of time at each frequency.  With the linear sweep it does that.  With the logarithmic sweep it spends more time at the lower frequencies and much less time at the higher frequencies, resulting in unequal energy.


Discrete Wavelet Transform (DWT) Filter Bank I have some stumbling block in my thesis writing.

Do we use the same filter pair while implementing DWT filter bank with downsampling of the filter output, or filters do change also from level to level?
Though on the Wiki page and in the book "Biosignal and Medical Image Processing, 3rd edition" by J. L. Semmlow and B. Griffel there are same labelling of filters at each decomposition level on the filter bank schemes:

,
the description in Wikipedia states: "The filter output of the low-pass filter g in the diagram above is then subsampled by 2 and further processed by passing it again through a NEW low- pass filter g and a high- pass filter  h with half the cut-off frequency of the previous one.."(c)

This resulted in big confusion for me. Logically, as I have been thinking before, for DWT we "travel" in frequency domain by signal downsampling and leaving the filters untouched (while e.g. in Stationary Wavelet Transform it is achieved through upsampling the filters themselves without signal modification).
The word "new" in Wiki frustrated me a lot.

Please, help me to resolve this issue.

With respect,

-Andrey.

In the DWT scheme, whether it is the classical $2$-band or the $M$-band wavelet setting, the very same analysis filter bank (lowpass/highpass + subsampling) is used at each level. Under this condition, one can derive the cascade algorithm that provides  the spectrum of the scaling function:
$$\Phi(\omega)= \prod_{k=1}^\infty \frac {1} {\sqrt 2} H\left( \frac {\omega} {2^k}\right) \Phi^{(\infty)}(0)$$
where iterated half-cut-off frequencies are apparent.

However, this is especially important when one wants to address the properties of the underlying wavelet, or cascade many wavelet levels. In practice, you can easily choose a different set of perfect reconstruction filter bank at each level, for several reasons:

actual filters down the scales tend to become larger, sometimes too much with respect to the signal features, and the signal length.  When a convolutive filter becomes longer than the data itself, you run into troubles, especially in regard of boundary handling (symmetry, padding). The consequence can be high details coefficients caused mostly by wrapping around border artifacts.
filtered/down-sampled data don't have the same characteristics at each scale, so different filters can do a better job
do we really need multicale theory after all on 3 levels? In practice, it can be useful to use longer filters on the first levels, for better frequency separation, and shorter ones on lower levels, for the multiscale effect.

This is not mainstream, as one has to choose filters wisely, sometimes without many hint. However, there were works with DWT or GenLOT filters on the first scales and DCT or Haar on the lower resolutions.  JPEG XR standard somehow behaves in that philosophy. One of the papers I have met on this direction is: On the frame bounds of iterated filter banks, 2009.


reference parseval theorem limited signal some days ago I asked here parseval for a continuos but limited signal  if the Parseval can be applied for limited signal.

Can you recommend me a book or a paper that I can use as reference for this?

Can someone show me the demonstration considering a limited domain

Thank you

I don't think you will find a reference for this. Your are simply stating that "integrating over a whole bunch of zeros results in 0". In other words, if
$$x(t) = 0, x\not\in [0,T]$$
then
$$ \int_{-\infty}^{+\infty}x(t) dt = \int_{-\infty}^{0}x(t) dt + \int_{0}^{T}x(t) dt + \int_{T}^{+\infty}x(t)dt = \int_{0}^{T}x(t) dt $$
since
$$ \int_{-\infty}^{0}x(t) dt = \int_{-\infty}^{0}0 dt=0$$ and

$$\int_{T}^{+\infty}x(t)dt = \int_{T}^{+\infty}0dt = 0$$

If that doesn't work, than I don't know what would


Rapid changes in signal correspond to high frequencies - proof If I consider a generic aperiodic signal $x(t)$, how can I prove that rapid changes in signal correspond to high frequencies?

Thank you for your time.

One might argue that the Dirac delta $\delta(t)$ , is a signal that exhibits maximum change because all frequencies are uniformly increased, not just the high frequencies


How to connect $\omega$ in time domain to physical frequency? I'm a noobie to DSP, and recently I'm learning DSP for an urgent assignment. There's a problem that has been confused me for days. I hope anybody could shed some light on it, although it might seem to be silly.

The assignment asks me to filter a signal so that it's in $20-800$Hz. After I spent several days to learn DSP, I finally get some idea about frequency domain, DTFT, FIR, IIR, etc. But I still have no idea about how to proceed. Furthermore, I have a huge confusion about frequency: in the frequency domain, we always talk about the radian $\omega$ which is always less than $\pi$, how can I connect this $\omega$ to the physical frequency? Is it right to use $\omega/2\pi$?

There are several "frequency" variables in signal processing, the most common all being related by constant factors. I will use the notation from Digital Signal Processing by Proakis & Manolakis:

$F$: "Ordinary" frequency in units of Hz.
$\Omega$: Angular frqeuency, in units of radians/second; $\Omega = 2\pi F$
$f$: Normalized or "digital frequency". $f = F/F_s$, where $F_s$ is the sample rate in Hz. This is a formally unitless quantity, but sometimes its units are listed as being cycles/sample.  Sampling produces aliasing, which means that the spectrum of a sampled signal will be periodic with period 1 cycle/sample. For this reason DSP engineers only look at the spectrum of the signal in the principal period of $-1/2 \leq f < 1/2$.
$\omega$: Normalized radian frequency. $\omega = \Omega/F_s = 2\pi F/F_s$. Sometimes its units are listed as being radians/sample. Because of aliasing, it is only necessary to study the spectrum of a signal from $-\pi \leq \omega < \pi$ in the digital domain.

Other notions of the concept of "frequency" exist (such as instantaneous frequency, and a notion of an "instantaneous spectrum"), but these are more advanced concepts that you probably don't need to worry about yet.


Discrete implementation of the PI controller I have been implementing discrete PI controller in the incremental (velocity) form in C++.
I have been looking for the anti-windup mechanism. One idea which I have can be described by
following pseudocode:
    integral_gain = proportional_gain*sampling_period/(2*integral_time_constant);

    // e(k)
    error = reference_value - actual_value;

    // dup(k)= kp*[e(k) - e(k-1)]
    proportional_increment = proportional_gain * (error - error_previous);
    // trapezoidal integration rule i.e. dui(k) = ki*[e(k) + e(k-1)]
    integral_increment = integral_gain * (error + error_previous);
    // du(k) = dup(k) + dui(k)
    action_increment = proportional_increment + integral_increment;
    // u(k) = u(k-1) + du(k)
    tmp = action_previous + action_increment;

    // anti-windup
    if (tmp > action_max) {
      action = action_max;
    } else if (tmp < action_min) {
      action = action_min;
    } else {
      action = tmp;
    }

    // u(k-1) = u(k)
    action_previous = action;
    // e(k-1) = e(k)
    error_previous = error;

In math:
Let $x_k$ be the "action".  Further, define an anti-windup function as
$$f_{aw}(x) = \begin{cases}
x_{max} & x > x_{max} \\
x_{min} & x < x_{min} \\
x & \mathrm {otherwise}
 \end{cases} $$
Then calculate $x_k$ from:
$$
\Delta p = k_p (e_k - e_{k-1}) \\
\Delta i = k_i (e_k + e_{k-1}) \\
x_k = f_{aw} \left ( x_{k-1} + \Delta p + \Delta i \right)
$$
I am not really sure whether the anti-windup mechanism mentioned above will be appropriate despite the fact that the incremental form is used. Can anybody help?

Define "appropriate".
Yes, the output of the controller will be limited to between your action_max and action_min -- this is good.  Since your output is also your integrator state, you'll be effecting anti-windup.  If I'm not horribly mistaken, if your output is saturated, as soon as the error term starts pointing in the opposite direction of the limit it's at, then it'll back away from that limit.
So you have an anti-windup method that is a good candidate.
However, anti-windup is a nonlinear control law, which is appropriate to use when you're dealing with some plant nonlinearity or some optimality criterion that is not least-square-error.  Usually you're dealing with actuator saturation, and usually you just want a stable system and you're willing to let it take a bit extra-long to settle when it has to move a long way.  In that case then, yes, this control law may be "appropriate".
But ultimately, the "appropriateness" of a controller depends on your plant, how it's being used, and if you're working in a team, whether your colleagues agree that it's appropriate*.  It's an opinion call.  So you need to determine if it's appropriate by your lights, and if you're in a team, you need to convince the team that you're right.
* Because the best control rule ever isn't any good at all if your boss won't let you implement it.


calculating similarity scores between similar and different classes? I have images for 99 person for each person i have 5 images .
i have this function

[similarity] = ssim(img1, img2, K, window, L);where K,window ,L are previously determined with their values.

now , i need to use this function to measure intra-class similarities and inter-class similarities .

i mean i need to measure the similaity within and outside the class and put the results in 2 matrices , one contains the results of similarities within the same classes , and the other contains the results of similarities among the determined class and all other classes.

i wrote this in matlab :

%%%%for intraclass

K = [0.05 0.05];
window = ones(8);
L = 100;
img1=imread('E:\matlab\images\0001hv1.bmp');
img2=imread('E:\matlab\images\0001hv2.bmp');
img3=imread('E:\matlab\images\0001hv3.bmp');
img4=imread('E:\matlab\images\0001hv4.bmp');
img5=imread('E:\matlab\images\0001hv5.bmp');
[s1] = ssim(img1, img2, K, window, L);
[s2] = ssim(img1, img3, K, window, L);
[s2] = ssim(img1, img4, K, window, L);
[s3] = ssim(img1, img4, K, window, L);
[s4] = ssim(img1, img5, K, window, L);
intra_class=[s1 s2 s3 s4];
histintra_class=hist(intra_class);
plot(histintra_class);title('intra_class');

%%%%for interclass

K = [0.05 0.05];
window = ones(8);
L = 100;
img1=imread('E:\matlab\images\0001hv1.bmp');
img2=imread('E:\matlab\images\0002hv1.bmp');
img3=imread('E:\matlab\images\0003hv1.bmp');
img4=imread('E:\matlab\images\0004hv1.bmp');
img5=imread('E:\matlab\images\0005hv1.bmp');
[s1] = ssim(img1, img2, K, window, L);
[s2] = ssim(img1, img3, K, window, L);
[s3] = ssim(img1, img4, K, window, L);
[s4] = ssim(img1, img5, K, window, L);
inter_class=[s1 s2 s3 s4];
histinter_class=hist(inter_class);
plot(histinter_class);title)'interclass');

now how can i implement this on all my images automatically , i mean how can i read all images and finding intra and inter class similarities and get automatically the 2 matrices which are inter class and intra class

regards

From a general perspective, you are trying to take SSIM - which is image difference measure. But it is probably a bad idea if you are looking for measure the intraclass, distance.

You can look at distance measures like Euclidean distance, Mahalonbis distance (if you have co-variance metrics) or others.

As a very simple example : To compute Euclidean distance $L_2$, which is computed between the centroid of the Clusters - as $ Distance = \sqrt{(x_a- x_b)^2 + (y_a- y_b)^2 }$

NOTE: here the dimension are NOT image dimension but are rather the dimension of the feature vector.

  I WILL EXPAND THIS ANSWER WITH YOUR SPECIFIC NEEDS.


Symbol rate and bit rate Nyquist's formula for maximum channel capacity (noiseless channel): $C=2B \log_2(M)$

Shannon's formula (noisy channel): $C=B \log(1+S/N)$.

I need to distinguish between the symbol rate and the bit rate.
Ragarding the symbol rate, do the two formulas says that one cycle of a signal (i.e 1 Hz) cannot carry no more than 2 symbol/cycle?
Regarding the bit rate, Nyquist formula doesn't put a limit on the number of bits that can be put in each symbol, while Shannon formula takes that consideration.

Is my understanding correct?

What Nyquist says is that $$R_p = 2B\,\text{ Bd},$$ where $R_p$ is the number of pulses per second transmitted, $B$ is the available bandwidth, and the units are Bauds or symbols/sec ($\text{Bd}$). Here it is assumed that the pulses being used are sinc pulses; otherwise, the pulse rate will decrease. In this sense, $2B$ is an upper bound on the pulses per second that can be transmitted. It is correct to say that you can transmit two pulses for every Hz of bandwidth available.

Now, each pulse can carry several bits, by means of its amplitude. If you allow, say, two amplitudes (for example, -1 and 1), then each pulse carries one bit and $R_b=R_p$, where $R_b$ is the bit rate. In general, if you allow $M=2^k$ amplitudes, you'll be able to transmit at a bit rate $R_b=kR_p$.

Shannon capacity does constrain the number of bits that you can transmit per second; furthermore, the usual formulation assumes sinc pulses.


simplest Notch filter in the world What do you think of this notch / band reject filter for a very very narrow band ?

Original signal x[n]
--> FFT
--> then I just put "0" in the FFT bin of the frequency I need to reject (notch filter)
--> inverse FFT
--> output signal y[n].

How is called such a filter ? Does it really work, or is it too simple to be true ?

The frequency response of a single FFT bin filter looks like a Sinc function, which has a massive amount of overshoot or ripple at frequencies between FFT bins.  So your filter is only useful if you can strictly guarantee that the input to the FFT only contains unmodulated frequencies that are strictly and exactly periodic in the FFT aperture length (e.g. not a window on longer data of unknown periodicity or a non-periodic stream).

A notch filter is equivalent to adding the inverse of the equivalent bandpass filter to the signal.  If the bandpass has ripples, they will get added to the result of the notch, thus potentially adding a lot of noise.

So you end up with a not-so-simple (as it requires pairs of FFT computations) that likely produces terrible results.


Linear Time-Invariant system frequency response and impulse response given input and output in time domain I am having difficulty with this question specifically with the $\sin()$ part. So here goes nothing:

Consider a Linear Time-Invariant system with input $x[n]$ and output $y[n]$. When the input to the system is
$$x[n]=5 \sin(0.4\pi n)/(\pi n) + 10 \cos(0.5 \pi n)$$
and corresponding output is $$y[n] = 10 \sin(0.3 \pi (n-10))/(\pi (n-10))$$

Determine the frequency response $H(e^{jw})$ and the impulse response $h[n]$ for the Linear Time-Invariant system.

Thanks in advance.

A classical question on signals and systems, for which I provide you one method of solution.

Now given an LTI system, and its input $x[n]$ and output $y[n]$, the fundamental approach to determine the impulse response $h[n]$ involves the utilisation of a frequency domain approach to assert $H(e^{jw}) = Y(e^{jw})/X(e^{jw})$

In this problem however, I would choose a different approach which involves more of an analysis by inspection rather than pure mathematical operations, because of the fact that the input $x[n]$ involves a cosine term whose DTFT is includes impulses and will be problematic in the division process.

The method proceeds by drawing the DTFTs of input and output:
The input $x[n]$ involves two terms one of which is $$5 \sin(0.4 \pi n)/ \pi n$$ whose DTFT is a zero phase ideal lowpass filter of cutoff frequency $W_c = 0.4 \pi$ with a gain of 5.

And the second term of input is $$10 \cos (0.5 \pi n)$$ whose DTFT is modeled as two impulses at $$ 5 \pi \delta(\omega + 0.5 \pi) +  5 \pi \delta(\omega - 0.5 \pi) $$ Because of linearity property of DTFT, these two DTFT are added to represent the DTFT of the input $x[n]$.

Then, look at the DTFT of the output $$y[n] = \frac {10 \sin(0.3 \pi (n-10) )}{\pi (n-10)} $$ which is again an ideal low pass filter of gain = 10 and $W_c = 0.3 \pi$ with a phase term of $e^{-10 \omega}$ hence $Y(e^{j\omega}) =10 e^{-j10 \omega}$ when $|\omega| < 0.3 \pi$ and zero elsewhere.

Now we want to find a frequency response $H(e^{j\omega})$ such that when it is multiplied by the DTFT $X(e^{j\omega})$ of the input $x[n]$, it should produce the DTFT $Y(e^{j\omega})$ of the output $y[n]$.

The observation of the input and output DTFTs yields the conclusion that one such $H(e^{j \omega})$ is a linear phase ideal lowpass filter of gain = 2, $W_c = 0.3\pi$ and phase term of $e^{-j10 \omega}$, whose time domain correspondence is found by inspection as $$h[n] = \frac {2 \sin(0.3 \pi (n-10))}{\pi (n-10)}$$

The justification of this solution can be made by observing the effect of this filter on its input: it will completely supress any signal above the frequency $\omega = 0.3 \pi$ and it will multiply any signal below that frequency by 2, while adding a phase shift of -10 radians which reflects itself as a group delay of 10 samples in (n-10) term.


Calculating the output of a pole eigen signal in a difference equation Let an IAR system be defined by the following difference equation:
$$y[n]-\frac{1}{4} y[n-2]=x[n]+3x[n-1]$$
and an input signal $x[n]=(-0.5)^n$.
The transfer function is defined as $H^z(z)=\frac{1+3z^{-1}}{1-0.25 z^{-2}}$. The value $z=-0.5$ is a pole of the system, hence, it is not in the ROC. We cannot use $y[n]=H^z(-0.5)x[n]$ to solve, as we would for any other eigen-signal.
Is there any other way to solve this? What can we say about the output? Is it infinity? is it undefined?

If the input is $x[n]=\left(-\frac12\right)^n$ for all $n$, then the output is infinite, because the convolution sum doesn't converge. Note that also the input grows without bounds for negative $n$.
If, on the other hand, the input starts at $n=0$, i.e., $x[n]=\left(-\frac12\right)^nu[n]$ then you can use the $\mathcal{Z}$-transform to compute the output:
$$Y(z)=X(z)H(z)\tag{1}$$
With $$X(z)=\frac{1}{1+\frac12 z^{-1}}\tag{2}$$
you obtain
$$\begin{align}Y(z)&=\frac{1+3z^{-1}}{\left(1+\frac12 z^{-1}\right)^2\left(1-\frac12 z^{-1}\right)}\\&=\frac{a}{\left(1+\frac12 z^{-1}\right)^2}+\frac{b
}{1+\frac12 z^{-1}}+\frac{c}{1-\frac12 z^{-1}}\tag{3}\end{align}$$
with some real-valued constants $a$, $b$, and $c$.
Inverse $\mathcal{Z}$-transform of $(3)$ gives
$$y[n]=-2a\;n\left(-\frac12\right)^n+b\left(-\frac12\right)^n+c\left(\frac12\right)^n\tag{4}$$

EDIT: First, to clarify, I use a basic example that shows that a double pole in the $\mathcal{Z}$-transform is no problem. Assume a system with impulse response $h[n]=(1/2)^nu[n]$ and an input signal $x[n]=(1/2)^nu[n]$. Computing the output in the time domain is straightforward:
$$\begin{align}y[n]&=u[n]\sum_{k=0}^n\left(\frac12\right)^k\left(\frac12\right)^{n-k}\\&=u[n]\left(\frac12\right)^n\sum_{k=0}^n1\\&=(n+1)\left(\frac12\right)^nu[n]\tag{5}\end{align}$$
The $\mathcal{Z}$-transform of the output sequence is simply given by the multiplication of the $\mathcal{Z}$-transforms $H(z)$ and $X(z)$:
$$Y(z)=\frac{1}{\left(1-\frac12z^{-1}\right)^2}\tag{6}$$
The inverse $\mathcal{Z}$-transform of $(6)$ equals the result $(5)$, which is very easy to confirm.

The probable cause of the OP's confusion could be the fact that for an input signal $x[n]=z_0^n$, $-\infty<n<\infty$, the output is given by
$$y[n]=H(z_0)z_0^n$$
if and only if $z_0$ is inside the ROC of $H(z)$, because otherwise the convolution sum doesn't converge and the expression $H(z_0)$ is meaningless. However, in the computation of the outputs $y[n]$ in the examples above we never evaluate $H(z)$ for any $z$ outside the ROC, so the situation is entirely different.


Applying filter in scipy.signal: Use lfilter or filtfilt? I saw in a SO thread a suggestion to use filtfilt which performs backwards/forwards filtering instead of lfilter.

What is the motivation for using one against the other technique?

filtfilt is zero-phase filtering, which doesn't shift the signal as it filters.  Since the phase is zero at all frequencies, it is also linear-phase.  Filtering backwards in time requires you to predict the future, so it can't be used in "online" real-life applications, only for offline processing of recordings of signals.
lfilter is causal forward-in-time filtering only, similar to a real-life electronic filter.  It can't be zero-phase.  It can be linear-phase (symmetrical FIR), but usually isn't.  Usually it adds different amounts of delay at different frequencies.

An example and image should make it obvious.  Although the magnitude of the frequency response of the filters is identical (top left and top right), the zero-phase lowpass lines up with the original signal, just without high frequency content, while the minimum phase filtering delays the signal in a causal way:

from __future__ import division, print_function
import numpy as np
from numpy.random import randn
from numpy.fft import rfft
from scipy import signal
import matplotlib.pyplot as plt

b, a = signal.butter(4, 0.03, analog=False)

# Show that frequency response is the same
impulse = np.zeros(1000)
impulse[500] = 1

# Applies filter forward and backward in time
imp_ff = signal.filtfilt(b, a, impulse)

# Applies filter forward in time twice (for same frequency response)
imp_lf = signal.lfilter(b, a, signal.lfilter(b, a, impulse))

plt.subplot(2, 2, 1)
plt.semilogx(20*np.log10(np.abs(rfft(imp_lf))))
plt.ylim(-100, 20)
plt.grid(True, which='both')
plt.title('lfilter')

plt.subplot(2, 2, 2)
plt.semilogx(20*np.log10(np.abs(rfft(imp_ff))))
plt.ylim(-100, 20)
plt.grid(True, which='both')
plt.title('filtfilt')

sig = np.cumsum(randn(800))  # Brownian noise
sig_ff = signal.filtfilt(b, a, sig)
sig_lf = signal.lfilter(b, a, signal.lfilter(b, a, sig))
plt.subplot(2, 1, 2)
plt.plot(sig, color='silver', label='Original')
plt.plot(sig_ff, color='#3465a4', label='filtfilt')
plt.plot(sig_lf, color='#cc0000', label='lfilter')
plt.grid(True, which='both')
plt.legend(loc="best")


Histogram Equalization I am reading the etext on Digital Image Processing by Gonzalez and under the section on Histogram Equalization, an equation which is apparently basic is shown. I have no clue what the equation means and where it is from. I hope someone can point me to it in wiki or something. I have been googling for awhile but I have yet to get any meaningful results.

Assuming that in your question $T(r)$ denotes a map $:\mathbb{R}\to[0,1]$ and assuming that it is continuous and strictly increasing, then it is obviously an isomorphism among $r\in[0,1]$ onto its image.
As a consequence you get the following result: let $r$ be any value in $[0,1)$, and let $a=T(r)$, then for any $\rho\in[r,r+dr]$ it is $\alpha = T(\rho) \in [a,a+da]$, having denoted $a+da = T(r+dr)$. Therefore the probability of $r$ and $a$ to belong their infinitesimal intervals (which is actually the definition of PDF) is the same. hence as $dr\to 0$, also $da\to 0$ and you can write
$$
p_r(r) dr = p_a(a) da
$$
which yields the result


Detecting a Signal Inside a Time Series Using MATLAB's $ \tt xcorr() $ I would like to find the exact replica of signal A in time series B. I am working with the simple toy example below, in which A is located starting at the 4th position in B.

I thought that using xcorr would achieve this but I am clearly failing somewhere. xcorr doesn't seem to be the right function despite the MATLAB help notes on xcorr and comments such as:

  Cross-correlation is usually the simplest way to determine the time lag between two signals. The position of peak value indicates the time offset at which the two signals are the most similar.

A = [1 2 1];
B = [0 9 8 1 2 1 7 3 1 0];
[c, lags] = xcorr(A,B)
lags(c==max(c))
c =
Columns 1 through 12
  0.0000    1.0000    5.0000   14.0000   18.0000   11.0000    6.0000   12.0000   26.0000   26.0000    9.0000    0.0000
Columns 13 through 19
  0.0000    0.0000    0.0000   -0.0000   -0.0000   -0.0000   -0.0000
lags =
  -9    -8    -7    -6    -5    -4    -3    -2    -1     0     1     2     3     4     5     6     7     8     9
ans =
  -1     0

The max correlation values are located at lags of -1 and 0, which doesn't match the offset required to get A to match B. This leads me to ask:

When is xcorr useful for finding the exact signal in another signal? 2. And what function should I be using to find the location of A in B?

Ideally I would like to test this method on more complicated signals as well as look for methods of detecting "similar" signals inside other (longer) signals.

The function xcorr calculates the correlation of 2 signals.
The correlation is known to be a good (The MLE) for delay estimation under Gaussian Noise.

Yet, as can be seen in your data you're not using it in the cases it meant to be used.
If we assume you have a model of a known signal with Additive White Gaussian Noise (AWGN or any other Additive White Noise), than the SNR is so terrible you can not expect it perform well.

If your signal model is additive deterministic signal (The signal [0 9 8 0 0 0 7 3 1 0] in that example) then you should look for a different approach (Though xcorr can be part of it).

Any tool you use, you must be aware for the model it was designed for.


How identify echo patterns in audio samples How do I estimate if in recorded sound samples is present an echo caused from walls reflection ?
(maybe there are multiples echoes: delayed, attenuated versions of the signal itself).

I'm recording several sound samples in my room and I suppose that in the recorded audio samples is present an echo image of the original sound.
How do I verify if this is true and in case, how can I calculate the samples delay (ms echo reflection time from the walls) ?

I plot the cross-correlation of the samples with itself (xcorr in MATLAB, about 40000 samples, ~1 second) but I don't know which pattern to look for to check if one or more echoes are present.

If I can understand the presence of the max peak in the middle of the graph, what are all the others equally spaced peaks on the both sides of the xcorr plot ?
Eventually, how do I estimate the echo parameters from the recorderd samples (ms delay or samples delay) ?

The side peaks could be the result of periodic structure in your input signal (say around on the order of 8-10Hz, plus a 2nd harmonic).

If you are creating the sound samples, try playing pseudo-random noise, which has good autocorrelation properties for a long enough sample, notably tight peaks & aperiodic.

2nd John's answer. You can do a back-of-the-envelope estimate of the max echo lag (assume a few reflections) using the speed of sound & room size.


Applying a window to input signal I am sorry if I say any thing silly .Please forgive me.I want to apply a window to input signals.I have two signals each of which contains zeros (462 samples) and non zeros data (38 samples).I want my window to take only the non zeros data (or data above certain threshold).My ultimate goal is to reduce the processing time and resources by taking only non zero data .I will use this data in finding cross correlation and finally the delay between those signals.Can any one comment How I can do that (some hints about matlab code).
Thanks in advance and I am sorry for any inconvenience!

% copy and paste into matlab command line

% say you vector is this:

a = randn(100, 1);

% you take out these small values:
threshold = 0.3;
a = a(abs(a) > threshold); % or try a(abs(a)<threshold) = [];


How to derive the stationary Kalman filter predictor? In its chapter on Kalman filters, my DSP book states, seemingly out of the blue, that the stationary Kalman filter for a system

$$\begin{cases} x(t+1) &= Ax(t) + w(t) \\
y(t) &= Cx(t) + v(t)
\end{cases}$$

has the predictor

$$\hat{x}(t+1|t) = (A-A\bar{K}C)\hat{x}(t|t-1) + A\bar{K}y(t)$$

and stationary state vector covariance and Kalman gain

$$\bar{P} = A\bar{P}A^T - A\bar{P}C^T ( C\bar{P}C^T + R )^{-1}C\bar{P}A^T + Q$$
$$\bar{K} = \bar{P}C^T(C\bar{P}C^T+R)^{-1}$$

where $Q$ and $R$ denote the covariances of the input noise $w$ and measurement noise $v$, respectively.

I can't see how to arrive at this from the minimum variance predictor. Could someone explain it to me, or point me to a resource that derives the expression? This is the time-variant minimum-variance filter, which I can derive:

$$\hat{x}(t+1|t) = (A-K(t)C)\hat{x}(t|t-1) + K(t)y(t)$$
$$P(t+1|t) = A\left(P(t|t-1)-P(t|t-1)C^T(CP(t|t-1)C^T + R)^{-1}CP(t|t-1)\right)A^T+Q$$
$$K(t) = AP(t|t-1)C^T(CP(t|t-1)C^T+R)^{-1}$$

I'm just unsure about how to go from here to the stationary filter above.

Update: I can see that substituting $\bar{P} = P(t+1|t) = P(t|t-1)$ and $K(t)=A\bar{K}$ into the time-variant filter results in the stationary filter, but why multiply with $A$? Is this just a symptom of an unfortunate choice of notation, meaning that either $K$ or $\bar{K}$ doesn't really denote the Kalman gain?

Your derivations are correct.

$\bar P = P(t|t-1)$ and $K(t) = A \bar K$

Is this your confusion:

Why didn't they have the term $t|t-1$ in the Kalman Gain and Covariance Matrix Expressions?
How can this be "stationary" when your derivation shows that it is time varying?

Bad choice of notation on the book's part

Let's look at the expression:
$ \bar P = A\bar PA^T - A\bar P C^T(C\bar P C^T + R)^{-1}C\bar PA^T + Q $.
The fact that $\bar P $ is a function of itself shows a recursive relationship. In other words, it uses its past values. So, it is NOT the same for all time instants -  It changes at every iteration.

Misunderstanding of the word "stationary".

When the author of the book said "stationary" he/she did not mean that $P$ and $K$ have the same value at all time instants. Instead, the author wanted to emphasise that the expressions for those to values are the same for all statistical realisations.  Stationarity is a statistical concept which means that the Statistics of the system is same all the time. Look at the expressions for $\bar P$ and $\bar K$ again. They depend only on
\

Previous values of themselves
Transition matrices $A$ and $C$ which are deterministic and in your case time-invariant ($A$  and $C$ are the same at all times)
$Q$ and $R$ which are the noise covariance matrices. These 2 matrices describe the statistics of the noises and are the same in all realisations and time instances.

The Kalman gain, $K$, and the state covariance Matrix $P$ will have the same value for all realisations of this random process. (Side Note: None of these 2 terms depend on the measurements, $y$. So they can be computed beforehand.)

Conclusion:

The "time-variant" equations you derived were equivalent to the ones in the book. Besides, the notational differences, there was a slight misunderstanding on your part regarding what changes and what doesn't.


Transfer function: from poles and zeros to polynomial coefficients I understand how this transfer function was solved except I don't know how to get the 0.87.
Any help would be appreciated.

In short: $(0.937)^2=0.878$

A bit more detail:

$$\begin{align}(z-re^{j\theta})(z-re^{-j\theta})&=z^2-r\left(e^{j\theta}+e^{-j\theta}\right)z+r^2\\&=z^2-2r\cos(\theta)z+r^2\end{align}$$


Image Segmentation applications I wanna experiment k-means clustering method on different kind of images, so I am trying to find different kind of images used in image segmentation such as MRI images.I want to gather some more categories.
Any suggestion would be gratefully appreciated.

I would search for "segmentation benchmark dataset" or so. One of the first results seems to be useful for your task:

https://www.eecs.berkeley.edu/Research/Projects/CS/vision/grouping/segbench/BSDS300/html/dataset/images.html


Find identical audio signal I have a long audio file (in wav) and a short audio file (in wav) that has been copied from the long file. Now I need to find the correct location where the short signal has been copied from. I tried to compare the sample data but for some reason (possibly resampling or introduction of short fades) it doesn't match up. What is a good and fast way to solve this?

If you are trying to find a signal inside another signal, Correlation is a useful technique.  You will get a large peak anywhere the signals match.  The easiest way to do it is to flip the shorter signal, such that samples 0 & N-1 are swapped, 1 & N-2, etc, then convolve the two signals.  This should attenuate any effects of noise or edge condition resulting from resampling or fades, unless the phase data got mangled, in which case you’ll probably need something more sophisticated.