Change phase of the output of a buffer I have a question about changing the phase of a buffer.

Question - How to change the output phase of a buffer by 90 degrees, 180 degrees, and 270 degrees?

I know to change the phase by 180 degrees, we just invert the signal using an inverter. But what about changing it by 90 degrees?

Thank you for your help!

If you are referring to a buffer as in a logic gate, this question may be more appropriate on the electronics stack exchange site but here are some suggestions.

If it is a logic gate, you must be referring to a fixed frequency clock, and shifting the phase for this is quite trivial. Here are some approaches to doing this:

For 180° simply invert with a not gate (inverter).

For 90° and 270° you can double the frequency with a frequency doubler and then invert the 2x clock and then divide by 2 with a flip flop (since it's output changes state on each rising edge). You need to gate the reset on the divide by twos with your original signal to ensure a consistent phase state (either 90° or 270°)

For 90° you can phase-lock a like frequency VCXO using an XOR gate as the phase detector which will have a lock point at 90°, invert this to get 270°.

For 90° you can implement a delay lock loop (similar to above but with a voltage variable delay line and compare one clock output to the delayed one with the xor gate, integrate it's output and use that to control the delay line for a first order DLL which will lock to 90°)

For non precise (temperature sensitive) approximate phase shifts, cascade buffers with RC filters in between.

For variable phase shifts there are plenty of delay lines on the market: https://www.digikey.com/products/en/integrated-circuits-ics/clock-timing-delay-lines/688


Which window function to use for QPSK signal? I am trying to Receive a QPSK signal and for Sampling I need to recover clock. Now QPSK Signal is generated at 2.4GBaud. It is than sampled at 80GHz.Than further it is downsampled to 10GS/s (Symbols per second).
Note: This signal have only one polarization.

Now I have Inphase and Quadrature Phase component of the signal(10GS/s). This signal is ofcourse in time-domain. Now I want to window the whole signal (the entire length in time). So which window function is suitable ?
Am I missing any details?

For purpose of recovering the clock, no window function is needed. Also to clarify, your signal is downsampled to 10GS/s which is samples per second, not "Symbols per second" as you specified the symbol rate to be 2.4 GBuad which is 2.4 G Symbols/second.

For QPSK clock recovery, I recommend considering the Gardner Loop in which case you would filter and downsample to 2 samples per symbol (4.8 GS/s). You can see the link below for more information on the Gardner Loop implementation and considerations.  I have seen very efficient implementations for closing the timing loop using polyphase filters; let me know if you would like more information on that (and update your title to not mention windowing) and I can add that here.

You will likely also have a carrier offset that you will need to track and remove. The Gardner is fairly insensitive to carrier offset, so once timing is acquired you can downsample to one sample per symbol (at the best decision locations) and use those complex samples to determine the carrier offset and correct for that in a carrier tracking loop (such as a digital Costas Loop).

See: Gardner Timing Recovery for Repeated Sybmols


Conditions on $v \in \mathbb{C}^{n}$ such that the complex to real DFT is sensible I'm using FFTW to compute some DFTs, and as a memory space optimization, you can give it a $v \in \mathbb{R}^{n}$ and it will return a $\hat{v} \in \mathbb{C}^{\lfloor n/2 \rfloor + 1}$ via the call fftw_execute_dft_r2c. (To those without context, it's just claiming that since $\overline{\hat{v}_{k}} = \hat{v}_{n-k}$, then you just don't need to care about the top half of the complex vector.)

In addition, you can take $\hat{v}$ back to $\mathbb{R}^{n}$ via fftw_execute_dft_c2r, but of course no one does that, they want to apply some filter to $\hat{v}$ and then call fftw_execute_dft_c2r.

What is the largest class of operations on $\hat{v}$ such that fftw_execute_dft_c2r make sense?

(Clearly, pointwise multiplication by the DFT of another real function is allowed, but it's not clear to me that every bandpass filter produces a sensible output and many other operations seem suspect.)

Obviously any operation is allowed that keeps the symmetry of the frequency domain data intact, i.e. any operation that doesn't add an imaginary part to the time domain sequence. And this is the case for all standard filtering operations unless you want to generate a complex valued (time-domain) signal by some operation. This could be the case if you want to generate an analytic signal, but in that case you actually make two signals out of one (the real part and the imaginary part in the time domain). So for normal low pass, band pass, high pass etc. filtering you don't need to worry. As long as you do things right, the symmetry of the time domain sequence should not be changed by these operations.


Calculating original signal from Discerete Fourier Transform I am trying to calculate the original equation using a DFT. I start with a equation, generate values from this equation and then get the dft of these values. The aim is to generate the original equation using the dft values.

Consider the following equation
$$f(t) = 1 + 0.75\cos(2\pi t/4) + 0.2\cos(4\pi t/4) + 0.5\cos(6\pi t/4)$$

original signal for $t=[0,1,2,3]\Rightarrow [2.45,0.80,-0.05,0.80]$

DFT: $[4,2.5,0.8,2.5]$

From DFT (i.e $[4,2.5,0.8,2.5]$) how do I retrieve the coefficient $0.75, 0.2$ and $0.5$ present in the original equation?

Cheers,
Ambi.

At first glance your first paragraph seemed wildly confusing to me. (The notion of "calculating" an equation makes no sense.) But now I see what your professor is asking of you. He wants you to determine the coefficients of an equation. He gave you a continuous $f(t)$ equation of a signal containing a DC component (1cos(0$\pi$$t$/4) and three sinusoids. And the 2nd and 3rd sinusoids have frequencies that are integer multiples of the 1st sinusoid.

So, your prof wants you to select the values and length ($N$) of a $t$ time sequence that will make those sinusoids reside exactly at an $N$-point DFT's bin centers. When that happens the bin centers' frequencies will give you the frequencies of the three sinusoids in your $f(t)$ equation. The peak values of the non-zero DFT spectral components (DC plus the three sinusoids) will enable you to establish the [1, 0.75, 0.2, 0.5] peak amplitude coefficients in your $f(t)$ equation.  Your prof wants you to learn how a sinewave whose peak amplitude is $A$ produces a DFT spectral component magnitude value of $M$ given the $N$ length of your DFT.  Your job is to see how $A$, $M$, and $N$ are related to each other. This is a good homework problem.

Your $t$ = [0,1,2,3] sequence's period is one which means your sample rate is 1 sample/second. And as Ollie said, that's too slow. Set your $t$ sequence to [0, 0.5, 1, 1.5, 2, ..., 6.5, 7, 7.5], compute a 16-point DFT, compute the spectral magnitudes, and see what happens
.


Blackman-Tukey Autopower equation I am looking over Blackman-Tukey method for Autopower.

It uses the DFT after a window is applied on a autocorrelation.

$$
Power Spectrum = \frac1{2\pi} \ \sum_{(k=-(N-1))}^{N-1}  w[k] R[k] e^{-i\omega k}
$$

Where does the $ 1/(2\pi) $ come from?
It's not from the Autocorrelation $R[k]$ or window $w[k]$ and it's not an inverse transform, it's forwards, so why require it here?
$e^{-i\omega k}$ why is this not the DFT exponential term of $e^{-i 2\pi kn/N }$
We seem to be missing an $n$?
Is this autocorrelation method 'biased' by division of $N$ (length of raw original signal) after the autocorrelation calculation?
I assume it must be to give a Power spectrum because if there was no
division of N it would be giving Energy spectrum.

The equation looks reasonable to me: The scaling by $\frac{1}{2\pi}$ is to have the result in units of normalized frequency $f$ instead normalized angular frequency $\omega$. What may be confusing is using the index $k$ from the autocorreation instead of $n$ since it would be the time domain variable for $R[k]$ and $w[k]$, while $k$ is often associated with the frequency index in the DFT. The OP is seeking the frequency variable to be in the exponent, in this case suggesting $n$, but instead the unit used is actually $\omega$ representing the fundamental normalized angular frequency (which for a DFT over $N$ samples would be given as the expected $2\pi n/N$ as in the OP's 2nd question).

In this case, there are $2N-1$ total samples consistent with a linear correlation over $N$ samples apparently given as:

$$R[k] = \sum_{n=0}^{N-1}x[n]x^*[n-k], \space\space\space k \in [-N-1, \ldots N-1]\tag{1}\label{1}$$

The Discrete Fourier Transform of the windowed autocorrelation function could then be (using a frequency index $\ell$ and assuming there are also $2N-1$ frequency samples):

$$X(\ell) = \sum_{k=N-1}^{N-1}w[k]R[k]e^{-j2\pi k \ell /(2N-1)}, \space\space\space \ell \in [-N-1, \ldots N-1]\tag{2}\label{2}$$

This will grow by a factor of $2N-1$, which we could normalize if desired by dividing by $2N-1$. As in any DFT, this is just a scaling and its used would be dictated by defining the units of the output.

Finally the normalized radian frequency (the fundamental frequency when $\ell = 1$) in this case is $\omega = 2\pi \ell/(2N-1)$, if substitute this in $\ref{2}$ we get a function of $\omega$:

$$X(\omega) = \sum_{k=N-1}^{N-1}w[k]R[k]e^{-j\omega k}\tag{3}\label{3}$$

Or if instead if we want a function of $f$, given $\omega = 2\pi f$:

$$X(f) = \frac{1}{2\pi}\sum_{k=N-1}^{N-1}w[k]R[k]e^{-j\omega k}\tag{4}\label{4}$$

I was not able to locate an actual copy of Blackman and Tukey's 1958 paper (which predated the FFT) but did find this additional variation from "Noninvasive Instrumentation in Medical Diagnosis" by Robert B. Northrop. In this variant, I would think that either the exponent should be divided by $(N+1)$ or the summation and number of output frequencies should be done over $N$ samples:

The OP asked in comments if this is best practice for estimating power spectral density. The power spectral density is given as the Fourier Transform of the auto-correlation function, and doing a linear correlation (which can be done by zero-padding the time domain sequence and doing a circular correlation) and windowing the resulting time domain sequence prior to computing a DFT are indeed best practice for estimating the Fourier Transform. The linear correlation avoids time domain aliasing, and the windowing reduces spectral leakage error. So this is standard practice.

One precaution not yet mentioned with windowing and spectral estimation: the equivalent noise bandwidth of a rectangular window is 1 DFT bin, which means when the spectrum in generally equally distributed the total power is equivalent to being distributed neatly in each DFT bin without double counting - and we can then use Parseval's Theorem to measure the total power (or the power across a band of frequencies) by summing the power in each bin. Windowing increases the equivalent noise bandwidth of each DFT bin (as given by the kernel, or discrete-time Fourier Transform, of the window), so that if we added up the power in each bin double counting would result and the total power would be over-estimated. The window itself also removes power from the signal since we are decreasing the signal over much of the time interval, resulting in a net loss (the loss from the window can be predicted based on the sum-square of the window samples, and the net loss is this loss plus the gain from the increase in kernel bandwidth). A single tone that is not spread out over multiple bins would only be effected by the power loss but not the gain from double counting, so the SNR is effected. This is called the processing gain of the window as detailed by fred harris in his classic paper "On The Use of Windowing in the Fourier Transform" and further detailed in this other posting:

How to calculate resolution of DFT with Hamming/Hann window?


Time recovery algorithm and a symbol with samples I am working on my graduation thesis. The topic is a receiver implementation in FPGA. I am doing an implementation of symbol timing recovery. For this part I have decided to choose the Gardner algorithm. I have found here discussions about Gardner, Mueller&Muller and other algorithms. All the posts were useful, but I still have the following questions.

I have noticed that my implementation of Gardner algorithm doesn’t work as intended in case of a big initial phase offset. Could you please give a suggestion, what I can do in such case?

I've already started a simulation in MATLAB with generating a signal with N symbols. Each symbol has to have 8 samples (as given by my supervisor).
I did the first part:
Input_signal = randi([0 15],2000,1); % mod order = 16, 2000 symbols

How can I implement the next part of this task, which is to create 8 samples for each symbol?

EDIT 1
I have used CORDIC algorithm and defined amplitude of each samples and phase.
Using them I can define a sample with max amplitude in each symbol.
How can I use this information in implementation the time recovery algorithm?

Preamble: This answer is about timing recovery in a sense of symbol synchronization, i.e. finding the proper sampling phase of a baseband signal.
Based on the stated requirement of only 8 samples per symbol, I will assume that you are employing a fully digital approach to timing recovery. That means that you have no control over the times of sampling by ADC. The fully digital approach looks like this (this nice picture is from slides here):

The "digital processor" is your FPGA-based implementation of the timing recovery algorithm, including the selected Gardner algorithm to find a phase error. "Sampler" - is an ADC with fixed sampling rate, 8 samples per symbol. "Analog processor" - this is an analog front-end, which is implemented in your board outside of the FPGA. And "signal in" is a baseband signal from your target communication channel.
In order to verify an implementation of your algorithm, you will want to model it, e.g. using a Matlab program or a Python script (whatever is more convenient). With the real signal samples from ADC you will be able to see the details of how your algorithm works in a long run, which otherwise is difficult to simulate in RTL and may be hard to debug when programmed on board. With a synthetic signal you can vary the parameters of the signal to see, how the algorithm reacts to it, and also you can plot a sigmoid for the TED, which can be very helpful.
About TED:
Timing error detector (TED) is intended to detect a phase error in sampling. Based on the sign of its output, the timing recovery scheme will incrementally shift the sampling phase by a little step value, until the timing error becomes zero (or close to zero). This step value is usually selected to be 1-2% of a symbol period. Note, that with a fully digital timing recovery approach, when there is no actual control of the sampling phase, an interpolator is used instead, to fill the gaps between the signal samples. As a result the phase shift is made for the recreated version signal and not the real one.
Sampling phase error is defined by a phase difference between the transmitter, which generates symbols, and base sampling frequency of ADC (actual sampling frequency/8). Ideally the algorithm should be able to converge for arbitrary initial value of that error.
The answer to the question:
In order to adequately model the initial phase error, the model must include not only a digital processor, but also a sampler. This means, that the input signal for the model should be continuous or oversampled. The degree of oversampling will define the resolution ability of your tests, i.e. how many different initial values you can provide. The goal can be 50-100 points per symbol period.
If you intended to use a real signal as an input of your model, you would use a logic analyzer to capture the signal at the input of the ADC with high enough frequency.
Alternatively, you may want to use a synthetic signal, that models a real one. Then you have to create its mathematical model. Usually you will start with random symbol values (as you did), then apply the modulation scheme (you will skip this if the symbol is a signal level), then apply a pulse shape (if the target transmitter does), then apply the distortions/pulse shaping relevant to your intended target communicating channel and analog front end (e.g. a low-pass filter can model the distortion of the signal in a twisted pair).
As you can see, the particular way of creating a synthetic signal will strongly depend on the type of the target communication channel.
Edit: Here is an illustration of what the resulting input data $x(k)$ for the TED may look like. It shows two different cases. You can see that they differ in the phase of ADC sampling. Also the examples show that the selection of 2 samples out of 8, which is necessary for the Gardner TED, is random.

The black lines are the exemplary eye diagram of the signal, the green line is the optimal sampling point and the blue lines are the ADC samples.


Autocorrelation of a Shifted Sequence Suppose I have a sequence $x[k]$ with $\mathcal{Z}$-transform

$$ X(z) = x_{0} + x_{1}z^{-1} + x_{2}z^{-2} + \ldots + x_{N-1}z^{N-1}$$

I know that for real-valued $x[k]$ the $\mathcal{Z}$-transform of the autocorrelation of $x[k]$ is given by:

$$ S_{XX}(z) = X(z) X(z^{-1}) $$

Now suppose I shift the signal by $n$:

$$ \tilde{X}(z) = z^{-n}X(z)$$

Would the $\mathcal{Z}$-transform of the autocorrelation simply be:

$$ S_{\tilde{X}\tilde{X}}(z) = z^{-n}X(z) z^{ n}X(z^{-1}) $$
$$ S_{\tilde{X}\tilde{X}}(z) = X(z) X(z^{-1}) $$
$$ S_{\tilde{X}\tilde{X}}(z) = S_{XX}(z) $$

However when I use the sequence $x[k] = \{ 1, 1, -1\}$, shift it by $n=2$ to get $\tilde{x}[k] = \{0, 0, 1, 1, -1\}$, I get:

$$ R_{\tilde{X}\tilde{X}}[k] =\mathcal{Z}^{-1}\{S_{\tilde{X}\tilde{X}}(z)\} = [0, 0, -1, 0, 3, 0, -1, 0, 0]$$

Which is equivalent to $S_{\tilde{X}\tilde{X}}(z) = z^{-2} S_{XX}(z)$ instead of

$$ R_{{X}{X}}[k] = [-1, 0, 3, 0, -1]$$

So how is my math wrong? What is the autocorrelation of a shifted sequence in relation to the original autocorrelation?

Hi: I think there are two different vewpoints to your question so I'll mention both.

First Viewpoint: ( Time Domain )

What you have written for the covariance is not wrong but you need to understand what you're calculating. Those numbers you are calculating are unscaled covariances of an MA(n) process  assuming A) that the variance of the noise term in the  MA(n) model, $\sigma^2$ = $1$ and B) the coefficients of the MA(n) model ( i.e: the $\theta$ ) are all 1.   So, the covariance of the shifted MA(n)series does stay the same because the 5 numbers ( A) the variance and B) the 2 autocovariances at lags 1 and lag 2 and lag -1 and lag -2 which are happily symmetric ) are staying the same. You get the other zeros because, for some reason, when you shifted the sequence in the time domain, some extra zeros arose. I'm not sure how that happened but it shouldn't. atleast not in the time domain.

Second viewppoint ( Frequency domain )

I'm just learning DSP so can't say a lot here but, from a time series frequency domain perspective,  if $y$ has  covariance generating $g_{y}(z)$, then that $X(z)$ filter can be viewed as an MA filter of order(n) so that if the MA(n) filter is applied to $y(t)$, resulting in say, $X(t)$, then the covariance generating function of $X(t) = B(z)B(z^{-1}) g_{y}(z)$ where $g_{y}(z)$ is the autocovariance generating function of $y$.

All of this quite nicely covered in Hamilton's "Time Series Analysis" if you have that book. I highly recommend it for time-series but it's mostly time domain material with frequency domain coverage somewhat scattered. Still, book is really great.


How do I use the step and ramp functions to create a specific signal? I am trying to figure out how to use the step $u(t)$ and ramp $r(t)$ functions to create the function $f(t)$ below:

$f(x) = \begin{cases}0 & \text{for } t < 0 \\ t & \text{for }-1 \lt t \lt 1 \\ 2 & \text{for }t \gt 1 \end{cases}$

The step function $u(t)$ is defined to be 1 when $t \gt 0$ and $0$ otherwise.

The ramp function $r(t)$ is defined to $t$ for $t \gt 0$ and $0$ otherwise.

So here is my thought process:

$f(t)$ has a line with a slope of 1 from [-1,1]. This looks just like the r(t) function shifted to the left by one. Hence, r(t+1). However, the problem arises after $t = 1$.  $r(t+1)$ continues on with a slope of 1 infinitely, but the graph below flattens out at $t=1$. Somehow, I need to cancel out the slope of 1 by adding with a slope of -1. The only way I can get a slope of -1 is by negating $r(t)$ and shifting it to the right by 1.

This yields $r(t+1) - r(t-1)$. This will cancel out the positive slope 1 for t > 1, but this creates a problem on left side because the slope of -1 continues on to negative infinity.

I feel like I'm over thinking this problem.

I would appreciate it if somebody could explain this problem in the simplest possible way.

Also, somebody told me there is an easy way to solve problems like this by looking at where the discontinuities are and where the slopes change.  The discontinuities imply using an amplitude with u(t) and slopes imply using amplitude with $r(t)$.

Would appreciate all / any advise.

Thanks.

You have a confilct in your definition for t < 0.

I'm going to assume you mean f(t) = 0 for t < -1, with the rest of the definition unchanged.

The simplest way I can think of doing this is to use u(t) as a tool to switch functions on and off over the desired ranges of time.  For example, if you want a function to be on over a range of time, say t1 to t2, multply the function by u(t-t1)-u(t-t2).  This turns the function on at t1 and off again at t2.

One representation of your function is:

$$f(t) = [u(t+1)-u(t-1)]t + 2u(t-1)$$

For this representation, we've used two functions: t and the constant function 2.
The function t is switch on at t1=-1 and off at t2=1.  The constant function '2' is switched on at t2=1 and never switched off.

If you need to use the ramp function instead of t, your function is:

$$f(t) = -[u(t+1)-u(t)]r(-t) + [u(t)-u(t-1)]r(t) + 2u(t-1)$$

Here we've used three functions: -r(-t), r(t) and 2.


A few conceptual questions about filter, pole, and bilinear I am working on a school project on converting a 6th order butterworth high pass filter to digital filter using bilinear transformation.

Just got a couple conceptual questions need to be clarified before I continue.

In an analog 6th order butterworth filter, the poles are the same for highpass and lowpass since poles are found in the denominator?
I have found 6 poles, and say my first pole is -0.259+ 0.966i if my cut off frequency is $1 \textrm{ rad/s}$. If my cut off frequency is changed to $10000 \textrm{ rad/sec}$. The radius of the unit circle would be $10000$ and my pole would be -2590 + 9660i?
Is there anyway someone could verify that once applied bilinear transformation without pre-wraping the frequency (professor said do it without pre-wraping first and compare to wrapped after), in the $z$-plane the poles would be:

0.1894 + 0.7499i
0.1894 - 0.7499i
0.1405 + 0.4077i
0.1405 - 0.4077i
0.1222 + 0.1283i
0.1222 - 0.1283i

I'm quite new on DSP, any help is appreciated.

The poles of low pass and high pass Butterworth filters are indeed the same if both filters have the same cut-off frequency. The difference between the two lies in the numerator. A low pass filter has a constant in the numerator, whereas a high pass filter has a constant times $s^N$, where $N$ is the filter order.

The poles of a normalized Butterworth low pass filter (i.e., one with cut-off frequency $\omega_c=1$) lie on a circle with center $s=0$ and radius $1$ in the complex $s$-plane. Of course they are all in the left-half plane:

$$p_k=-e^{j\frac{k\pi}{2N}}\tag{1}\\
k=\pm 1,\pm 3,\ldots,\pm (N-1),\quad N\text{ even}\\
k=0 ,\pm 2,\ldots,\pm (N-1),\quad N\text{ odd}$$

If you have a cut-off frequency $\omega_c\neq 1$, the variable $s$ is transformed as $s\rightarrow \frac{s}{\omega_c}$, so a factor $\frac{1}{s-p_k}$ of the transfer function becomes $\frac{1}{s/\omega_c-p_k}=\frac{1}{\omega_c}\frac{1}{s-\omega_cp_k}$. So the poles are indeed multiplied by the cut-off frequency $\omega_c$, i.e., they lie on a circle with radius $\omega_c$. Their angles of course remain unchanged.

Without pre-warping, the bilinear transformation is usually implemented as

$$s=2f_s\frac{z-1}{z+1}\tag{2}$$

where $f_s$ is the sampling frequency. Expressing $z$ in terms of $s$ gives

$$z=\frac{2f_s+s}{2f_s-s}\tag{3}$$

Equation $(3)$ can be used to determine the pole locations in the complex $z$-plane. If the cut-off frequency of the analog filter is $\omega_c$, its poles are located at $\omega_cp_k$, with $p_k$ given by $(1)$. So the poles in the complex $z$-plane are given by

$$z_{\infty,k}=\frac{2f_s+\omega_cp_k}{2f_s-\omega_cp_k}\tag{4}$$

You can use $(4)$ to check the pole locations you obtained from Matlab (and I don't doubt that they are identical up to numerical inaccuracies).


Complex nature of chirped signal In a book on optical solitons, the author says that a Gaussian signal pulse has the form-
$u(t) = \exp((\eta + i\beta)t^2)$
where $\beta$ is the chirp parameter.
My question is: what does the $i$ signify here? As a chirped signal is still real, what does the complex nature of the signal represent?

Depending on the context, the use of the complex form could be for mathematical convenience or for a no-kidding need for both real and imaginary parts.
When you factor the expression, you get
$$u(t) = e^{{\eta}t^2}e^{j{\beta}t^2}$$
Where the first exponential is a generic magnitude envelope, in this case Gaussian. The second exponential is the chirp itself and is where all the action is.
So to simplify things, lets assume that the envelope is ideal so that the signal is just the chirp
$$u(t) = e^{j{\beta}t^2}$$
To view the chirp, you can take either the real or imaginary part, and it looks something like this

In a system that does not use I/Q, the real part is what you would expect to transmit as your waveform. The type of system will determine whether using a real signal or some type of I/Q is best. I'm going to use a radar example here.
In frequency-modulated continuous wave (FMCW) radars, a real chirp like the one above is transmitted and has the form
$$x(t) = cos({{\beta}t^2})$$
Which is just the real part of the complex form. It is received after a delay and mixed with itself, and without going into the mixing process, produces a single frequency sinusoid that can be used to determine range. Here, using a real part only is practical. Using I/Q in FMCW is also beneficial (SNR improvement), but is not usually necessary and many systems do not use it.
Another type of radar, pulse-Doppler, benefits greatly from using the complex form. The same chirp is considered, except now the imaginary version is used. This is important because pulsed-Doppler radars usually operate on performing pulse compression, which is just correlating the transmitted waveform with the received one.
The autocorrelation of a complex chirp looks like

Using a complex waveform allows us to mix our signals to baseband, which give the classic autocorrelation responses we expect without additional mixing and filtering.


Integral of the Fourier spectrum? The integral

$$\int_{-\infty}^{\infty}|X(f)|^2df$$

of the absolute Fourier spectrum squared is the energy in the signal, but what about the integral of the 'simple' absolute Fourier spectrum?

$$\int_{-\infty}^{\infty}|X(f)|df$$

Does it represent anything?
I couldn't find any information on that point. Any hint would be appreciated!

If $x(t)$ and $X(f)$ are Fourier transform pairs then the value

$$\int_{-\infty}^{\infty}|X(f)|df$$

can serve as an upper bound on the magnitude of $x(t)$:

$$|x(t)|=\left|\int_{-\infty}^{\infty}X(f)e^{j2\pi f t}df\right|\le \int_{-\infty}^{\infty}|X(f)|df$$


The DC bin of an FFT should be the mean of all the other bins in the spectrum, why does the FFT produce 0? I am experimenting a bit with FFT's. From what I have read the first bin of an FFT is called the DC bin -- that is it is the mean of the other components in the FFT. I have also read that it should be close to 0.

When I actually perform the FFT and take the absolute value, what I find is that the first element of this is always 0. This makes sense mathematically if you consider the definition of the DFT itself.

Is there a best practice to reconstructing this first DC bin?

So it seems that because I was making my test FFT with a sine wave as the transient signal, this is what caused the 0 valued first FFT bin. I changed it to a cosine and got a non-zero first element...And the same with any random data. I should have thought a bit first.

Depending on your normalization factor the DC bin is either the sum of the signal values in your frame (1), the average (1/N), a "representation" (1/sqr(N)),

In the DFT formula the exponential value is $e^0$, or 1, for each term.


How to select the center frequency and filter type for 5-band digital audio equalizer? I have a small team project: use Matlab implement a Gui 5-band digital equalizer using FIR and IIR filters (2 cases) - The sample rate of all audio files is standard audio 44.1kHz. Coding is not my concern, however, I have 2 questions

(1) What are 5 center frequencies that I should choose ?
I saw a homework problems select 500Hz, 1.5Hz, 4Kz, 10kHZ, and 16KHz but I am not sure how the author came up with these number.

(2) Should I use 5 bandpass filters? Or should I use 1 lowpass, 3 bandpass, and 1 highpass filter?

The answer may depend on your application. I suggest you look at some commercial audio equipments with EQs to decide for those band frequencies. For example the old analog philips casette-type stereo system in my room has 5 band graph equalizers with the frequencies: [100-300-1000-4000-10000] Hz respectively. I use it as an amplifier for my laptop :)

If you design the first and last bands as Bass/treble shelve types, be careful to limit their cutoffs so that mid bands will have enough spectral intervals.


how to compute the signal passing from the low pass filter? I am currently trying to solve this question.

Let $x[n]=\cos(\frac{\pi}{2}n)$ and $h[n]=\frac{1}{5}\text{sinc}(\frac{n}{5})$.
Compute the convolution $y[n]=x[n]∗h[n],$ and write the value of $y[5].$
Hint - Given for the question is to compute the convolution in the frequency domain first.

when I calculated the Fourier transform of $y[n]$, I got
$$Y(e^{jw}) = X(e^{jw}) \times H(e^{jw})$$
$$X(e^{jw}) = (1/2) * [\delta(\omega-\pi/2)+\delta(\omega+\pi/2)]$$
$$H(e^{jw}) =
     \begin{cases}
       \text{1,} &\quad \text{if } |\omega| <= \pi/5\\
       \text{0,} &\quad \text{otherwise} \\
     \end{cases}$$
But what I don't understand is, since the range of $X$ is beyond the range of $H$, how do we multiply the two FT's? I must be misunderstanding the concept somewhere. Could anyone help explain to me how do we find the output?

Low pass filters get the name because they pass signals at low frequencies and attenuate signals at higher frequencies. The definition of low and high depend on the the cut-off frequency of the filter. In your case, the cut-off frequency is $\frac{\pi}{5}$. Your filter is the ideal low pass filter where signals above the cut-off frequency are completely zeroed out. The input signal is a perfect cosine signal so all of its energy is concentrated at one specific frequency. It happens to be at a frequency above the cut-off frequency, so what should happen?


Confusion understanding phase shift/delay? I am reading proakis, as shown highlighted in attached snap shot,there is '-'sign along with j,but still it is written (underlined red) implies shift in positive n direction,why and how positive?when we have negative sign with j?

Delay is the negative derivative of phase with respect to frequency. Phase that increases linearly negative as frequency increases will result in a positive (later) delay in time. This is given by the Fourier Transform property of time delay:
$$\mathscr{F} \{f(t-\tau)\} = e^{-j\tau \omega}F(\omega)$$
And is easy to prove from the Fourier Transform:
$$\mathscr{F} \{f(t-\tau)\} = \int_{-\infty}^\infty f(t-\tau)e^{-j\omega t}dt$$
$$ = \int_{-\infty}^\infty f(t-\tau)e^{-j\tau \omega}e^{j\tau \omega}e^{-j\omega t}dt$$
$$ = e^{-j\tau \omega}\int_{-\infty}^\infty f(t-\tau)e^{-j\omega (t-\tau)}dt$$
If we substitute $u=t-\tau$, then $du = dt$ and we get:
$$ = e^{-j\tau \omega}\int_{-\infty}^\infty f(u)e^{-j\omega u}du =e^{-j\tau \omega}F(\omega) $$
Also without going into the details of the Fourier Transform, consider how we subtract phase in the time domain to shift a sine wave in time to the right:


Are there preferred patterns for creating texture on objects to aid computer vision? In the comment to this answer using a laser diode and diffraction grating to provide texture on a surface was suggested to aid in height calculations in a computer vision system.

I believe the pattern that I'm familiar with is projecting a chessboard pattern on the object.  I was under the (incomplete) understanding that the deformation of the projected image was somehow used directly.  That is to say the formerly square pattern that is projected becomes a curved polygon and the transformation back to a square gave some information about the 3D structure.  Is this incorrect?

Are there preferred patterns?  What are the dependencies in choosing a pattern?  Does it depend on the shape, material, etc, of the target object or is it more a function of the amount of variability in the features?

There are some generally applicable advices, and some application-specific advices.

Shi and Tomasi's paper, Good features to track explains the criteria for choosing patterns: two-dimensional localizability, or "cornerness".

To put it simply, suppose you are trying to find an object at position (x,y), but instead the object appears in the image at (x + dx, y + dy). It is not very useful if our vision system can only tell us that "no, the position is wrong." Instead, we expect the vision system to be able to estimate the amounts dx and dy provided that it is not too far off.

A sharp point (dot) is the most cornerful, but it is also easily buried in random pixel noise. By following through with the mathematics, we learn that there are other patterns that are just as cornerful as a sharp point. (Think about a 1D "edge" being a 1D delta transformed by integration.)

Some applications would call for localizability in fewer, or higher dimensions.

Added 8/25

Two line-like patterns can also be "intersected" to yield a point during calibration, provided that lens distortion is not significant or has been parameterized.

In deblurring applications, a sharp point is often used for recovering the point spread function (psf). However, in theory any arbitrary-shaped objects could be used, provided that the ground truth is available to the calibration software.

In some applications, we would deliberately make the pattern un-sharp. Depth from defocus uses the bluriness to reason about the position of the focal plane relative to the object, which gives an estimate of the object depth.


Optimal pulse shape for minimal interference in adjacent frequencies I need to find an optimal pulse shape, that convoluted with matched filter on adjacent frequency will give minimal convolution/correlation (Minimal interference) for a given pulse length.

For example Pulse Length of 10 msec, 200 Hz separation (10,000 Hz and 10,200 Hz)

Image below demonstrating the problem.  Currently I have interference with adjacent signal of approx 0.05.

What would be the Optimal pulse shape for given pulse length?

The root-raised-cosine pulse-shape is widely used because

using the root-raised-cosine as both pulse shape and as matched filter yields the raised-cosine pulse-shape, and
the root-raised-cosine and raised-cosine pulse-shapes are both bandlimited.

The frequency response (Fourier transform) of the root-raised-cosine is equal to the square root of the frequency response of the raised-cosine, so these Fourier transforms are nonzero in exactly the same (bounded) intervals.

Suppose that $h_1(t)$ and $h_2(t)$ are two root-raised-cosine pulses with identical rolloff parameter ($\beta$) whose Fourier transforms are nonzero on non-intersecting intervals. Then matched filtering of of $h_m(t)$ with $h_n(t)$ yields
\begin{equation}
g(\tau) = \int_{-\infty}^{\infty}h_m(t)h_n(\tau - t)dt.
\end{equation}
The Fourier transform of $g$ is
\begin{equation}
G(\omega) ~=~ (\textrm{Fourier transform of $h_m\ast h_n$}) ~=~ H_m(\omega)H_n(\omega).
\end{equation}
Since $H_m$ and $H_n$ are nonzero on non-intersecting intervals, their product is zero for all $\omega$. The inverse Fourier transform of an all-zero Fourier transforms is the zero function, so $g(\tau) = 0$ for all $\tau$. Hence, using a root-raised-cosine matched filter will filter out root-raised-cosine pulses from other frequency bands.

See the answers to this question for more information.


Multiple short-range impulse responses aggregated to simulate a long-distance impulse response I'm a biologist and I was thinking of ways that impulse responses could be used to simulate how an animal sounds over various distances in different types of forests.
Ultimately, if a tool was constructed that would allow scientists to simulate animal sounds over various distances, this would be a pretty useful research tool.  Specifically, I'm interested in simulation of the reverberation effects, hence why I'm interested in impulse responses.
Process-based modelling an impulse response is not trivial, however, as an empiracle approximation I was wondering:

If I record an impulse response at 50m from a source in a forest, if I were to convolve a source signal twice, would that result in a reasonable approximation of what the sound would be at 100m distance?

The idea is that if I can accurately measure impulse responses over short distances in different forest types, these could be used to simulate reverberation over a large distance, including sound passing through multiple forest types.  Currently, I can't think of why this would be a bad approach... but I'm a biologist, so I'm hoping others can tell me if this seems like a flawed concept.

if I were to convolve a source signal twice, would that result in a reasonable approximation of what the sound would be at 100m distance?

Yes. For small number of cascades this is a reasonable assumption.
However, you can't cascade the 10 meter response 10 times to get a single 100 meter response. While the general shape and structure of the response from 20m to 30m is the same as that from 30m to 40m, the phase spectra and the fine structures of the reflections are uncorrelatated. If you cascade the same impulse response over and over again, you get spectral thinning and it's starting sound "metallic" and "hollow".
You can try to decorrelate short repsonses from the original or just to start with a parametric model that has tree density, tree reflectivity, leaf/brush absorption and obvioulsy distance as tweakable parameters.


Output of marginally stable systems In' b.p lathi's book of systems and signal it is written that
"marginally stable systems have one important application in the oscillator, which is a system that generates a signal on its own without the application of an external input. Consequently the oscillator output is a zero input response "
By this statement I have two problem:

if let's say oscillator system is a simple L-C circuit then how can it produce output without input and if it does produce then isn't it a violation of linearity?

and if initial conditions are zero then how can L-C  circuits have zero input response? And then in that case how does it generates signal without initial conditions?

An ideal LC circuit has no resistance, and hence no resisitive loss mechanism. Therefore if it starts to oscillate (due to, for example, a non-zero initial voltage or current stored in the L or C), then it will oscillate indefinetely on its own, yielding a situation as described in your question.
Zero input response of a linear system is the response under no input excitation but due to non-zero initial conditions.
Note the distinction between the terms linear time-invariant (LTI) and the broader term linear alone. The discussion in the Lathi's book is about linear systems in the more general sense.
Note also that a linear (but not necessarily LTI) system can have non-zero initial conditions, and hence zero-input response; but a strictly LTI system cannot have non-zero initial conditions, because one property of LTI systems is the initial rest property. Hence zero-input response of LTI systems is zero by definition.
So by definition a system with zero-input response cannot be LTI, but it can be linear nevertheless...


LP to BP transformation steps Considering the following problem, What i know to convert from LP signal to BP signal is to:

Draw the FD of the LP equation
From it I can draw the FD for BP
From this graph I can deduce the equation of the BP signal.

My problem here is that the equation has an Imaginary value, so how can i plot this equation and how can i deduce the graph of BP signal from it؟
because i searched and i found results about Hilbert transform and imaginary odd symmetry but i didn't figure that out.

HINT:

Just multiply the complex baseband signal with a complex carrier $e^{j2\pi f_ct}$, and take the real part to obtain the corresponding real-valued band pass signal. You can sketch either the magnitude or the real and imaginary parts separately.


Terminology: spectrum, spectrogram, spectrograph, sonogram, etc In term of proper or accepted naming conventions of DSP graphics or instrumentation output, what is the difference between the words spectrum, spectrogram, spectrograph, and similar terms, and what type of chart, graph, CRT display, or etc. does each best describe.

ADDED: Also I found the term sonogram used in a couple books for spectrum-vs-time graphics.  When might thus be appropriate in preference to one of the above terms, or vice versa?

It depends on context.

In signal processing, a spectrum (plural is spectra) shows the frequency content of an entire signal.  It's a 1-dimensional function of amplitude (vertical axis) vs frequency (horizontal axis):

Spectra are often shown with a logarithmic amplitude axis (such as dB), but this isn't necessary.

A machine that produces a spectrum is usually called a spectrum analyzer.  In other fields, the machine is called a spectrograph or spectrometer.

A spectrogram shows how the frequency content of a signal changes over time. It's a 2-dimensional function of amplitude (brightness or color) vs frequency (vertical axis) vs time (horizontal axis):

Sometimes this is called a sonogram.  The time and frequency axes are sometimes swapped.  If amplitude is shown as a 3D surface rather than using color, it's called a waterfall plot.

Confusingly, a machine that produces a spectrogram is also called a spectrograph, or spectrograph is used as a synonym for spectrogram.

Also the line is kind of blurred, because if you view the spectrum of a live signal on a spectrum analyzer, it's displaying the spectrum of small chunks of the signal and you're seeing how it changes over time, which is essentially the same thing as a spectrogram.

I think the important distinction is just the way they're displayed:  A spectrum is a 1D plot and a spectrogram is a 2D plot.


DIF based FFT scaling for OFDM Input I'm working on XILINX FFT IP(PG109) in pipeline mode, which uses radix 2 DIF Lattice to perform FFT/IFFT operation. How to perform scaling of each stage in FFT, such that the effect of saturation or underflow will be less?

Input Data: An OFDM Symbol from channel.

My work :
1) Data growth will be due to addition/subtraction and rotation of vector done by twiddle factor.

2) Because of addition/subtraction there will 1 bit growth(when extreme values involves)

FOR DIF ARCHITECTURE BIT GROWTH IN LAST TWO STAGES IS PURELY DUE TO ADDITION/SUBTRACTION AS TWIDDLE VALUES ARE +/-1, +/-j.

3) For rest first N-2 stages, data grows as per note 1 . The max. real/imaginary can go to 2*sqrt(2).

4)Always data growth in every stage can't lead to overflow.(when input values are small).

I have worked with 4096-point Radix-2 based FFT algorithms, and although one might think that we need to adjust the accumulation after every stage by rightshifting the accumulation by 1 bit.

But that leads to over scaling of the magnitude and if the input $x[n]$ has low amplitude, then the accumulation at few stages might become 0 after this division by 2.

I found the scaling 0xAAA as most optimal performing one, which means that we are adjusting the accumulation at alternate stage rather than after very stage.

0xAAA : 101010101010 -> Accumulated register is right-shifted by the bit at bit-position in 0xAAA corresponding to the stage.

You can try with 0xFFF also, and actually see that the FFT result is getting higher error at those FFT Coefficients which are having low amplitude.

Hope this helps.


Heisenberg Uncertainly Principle and wavelet transform I am using CMOR(complex morlet) wavelet in Fourier space in order to reconstruct my signal and also estimate damping and frequencies of the embedded modes. there are two main parameters in cmor. the Fb and Fc. bandwidth and center frequency, with increasing the Fb i get better result for frequency estimation, but i cant get a good damping estimation with high Fb. is it sth to do with Heisenberg Uncertainly Principle? and if yes, would you explain the role of fb and fc in the accuracy of frequency and damping estimations?
 The wavelet transforms provide a unified framework for getting around the Heisenberg Uncertainly Principle.

The HUP follows directly from the properties of the Fourier Transform, because time and frequency are orthogonal bases in which we can expand the co-efficient sequence of our signal.

In fact all pairs of orthonormal bases will have some kind of Uncertainty Principle associated with them.

In traditional Fourier analysis, the either the time axis or the frequency axis will be split into identical sized and spaced regions which permit us to gain information about the signal:

Wavelets allow you to analyse signals by splitting up the time-frequency plane into regions of different sizes.

If you are increasing the bandwidth (i.e. frequency), you will necessarily have worse resolution in time.


How to 'whiten' a time domain signal? I am trying to understand how exactly to implement what is known as a 'pre-whitening' filter or simply a 'whitening' filter.

I understand that the purpose is to make it have a delta as its autocorrelation function, but I am not sure how to do this exactly.

The context here is the following: A signal is received at two different receivers, and their cross-correlation is computed. The cross-correlation can look like a triangle, or some other godforsaken shape. Due to this, it becomes hard to find the peak of the cross-correlation signal. In this case I hear about having to 'whiten' the signals before a cross correlation is performed on them, such that the cross-correlation is now more delta-like.

How is this done?

Thanks!

Suppose you have signals $x(t)$ and $y(t)$ whose cross-correlation function $R_{x,y}(t)$ is not something you like; you want $R_{x,y}$ to be impulse-like.  Note that in the frequency domain,
$$\mathcal{F}[R_{x,y}] = S_{x,y}(f) = X(f)Y^*(f).$$
So you filter the signals through linear filters $g$ and $h$ respectively to get
$\hat{x}(t) = x*g$, $\hat{X}(f) = X(f)G(f)$, and $\hat{y} = y*h$, $\hat{Y}(f) = Y(f)H(f)$, and now their cross-correlation function is
$R_{\hat{x},\hat{y}}$ whose Fourier transform is
$$\begin{align*}
\mathcal{F}[R_{\hat{x},\hat{y}}] = S_{\hat{x},\hat{y}}(f) &= [X(f)G(f)][Y(f)H(f)]^*\\
&= [X(f)Y^*(f)][G(f)H^*(f)]\\ &= [X(f)Y^*(f)][G^*(f)H(f)]^*,
\end{align*}$$
that is, $R_{\hat{x},\hat{y}}$ is the cross-correlation of $R_{x,y}$ with
$R_{h,g}$.  More importantly, you want to choose $g$ and $h$ so that
the cross-spectral density $G(f)H^*(f)$ of $g$ and $h$ is the multiplicative
inverse of the  cross-spectral density $X(f)Y^*(f)$ of $x$ and $y$, or
something close to it.  If you have only one signal and one filter, then
you get the result given by Hilmar (with amendment as given by my comment there).
In either case, the issue of compensating for spectral nulls, or generally,
frequency bands where the signals have little energy still remains.


GNU Radio - PSK Mod block - unexpected constellation diagram I am simulating a QPSK modulation in GNU radio. Here is the flowgraph:

.

The source is a binary stream coming from the Random Source block. There are 4 PSK Mod blocks configured identically (4 constellation points - a QPSK modulation, Grey coding and Differential encoding turned on, excess bandwidth set to 0,35). The modulators differ in the parameter samples/symbol (sps), which is set to 2, 4, 8, and 16. The constellation sinks are configured identically. What is expected to be shown on the constellation diagram is a plot with 4 constellation points that could be the tips of a square or 4 clusters of points that are concentrated around those 4 points as a certain amount of noise or imperfection of the system is added. Instead, we have the following:

Constellation diagram when sps is set to 16:

Constellation diagrams for sps 2, 4 and 8 look the same but with less points on the sides of the rhombus and less points around tips of the rhombus that look like the extensions of the rhombus sides.

The question is why constellation diagrams look like this - why additional points are added and there are more and more of them as we increase the samples per symbol parameter so the diagram with sps 16 looks more like a rhombus than a 4 clusters of points? Why significant number of points appear around 4 tips of the rhombus looking like extensions of rhombus sides?

What you are seeing is the transitions from one constellation point to another.  In order to reduce the signal bandwidth, the baseband signal is low-pass filtered.  This causes the transitions to not be instantaneous (i.e. the I and Q are not square waves), so they take some time.  You are simply seeing those transitions.  The low-pass filtering also causes the overshoot that you see at the constellation points.

Just to see, I tried putting this into R. This is not going to get quite the same result because I am not sure of the parameters GNURadio uses, but it illustrates Jim's point.

As you increase the samples per symbol, the number of "transition" samples increases, making the edges of the constellation clearer.


Filter - Spatial Domain Versus Frequency Domain Note: This is a followup question of

  How can I create a custom band-pass filter?.
  Can you explain the following python source code?.
  Bandpass filter bank.

Suppose, I have an equation of a filter as follows,

$$
h(x, y) = \frac{1}{1 + 0.414 (x + y)^{2n}}
$$

I know that this is an equation in the spatial domain.

So, this is how I obtain a filtered image:

I create a, say, $3\times 3$ kernel using that equation where ($x$: from $0$ to $2$) and ($y$: from $0$ to $2$).
I pad the kernel and the image to the same size.
I apply FFT to kernel and the image.
Multiply them element by element and apply I-FFT to the result.
Convert the result to an image format like Bmp and so on.

.

Now, suppose, I have the following equation,

$$H(u, v) = \frac{1}{1 + 0.414(u + v)^{2n}}$$

I know that this is an equation in the frequency domain.

So, how am I going to create the kernel?

In the previous case, ($x$: $0$ to $2$) and ($y$: $0$ to $2$).

What would be the starting and ending values of $u$ and $v$ in the loop?
What would be the procedure to apply the kernel to an image?
How should I multiply the kernel and the image?

It appears you need to study a bit on convolutional filtering of images, specifically on overlap add/ overlap save methods.

From the links I can see your objective is to apply the filter defined in the second equation i.e, the frequency domain one. The first link's accepted answer seems to be wrong (i didn't go through the code, clearly since the images shown by the paper and answer do not match).

Now the method you are using to apply the filter in the spatial domain is wrong. I'm assuming by padding you mean zero padding. The objective of zero padding before applying fft is to increase the resolution in the frequency domain. So imagine I have a 1D sine wave of 50 Hz. Taking fft of this signal would give a peak at the 50 Hz point with some additional values at different frequencies. This is because the sine is not pure sine since the signal is limited in time. If you had an infinitely long signal you would get a single value at 50 Hz. Padding the signal and then taking fft would still give the similar result but with more resolution on what additional frequencies exists.

In contrast convolution filter is essentially convolving 2 signals.
$$
y(t) = \int_{\tau = 0}^{\inf} h(\tau)x(t-\tau)d\tau
$$
In 2D, if we have a 3x3 filter kernel, we first multiply the first 3x3 block of the input with the kernel and then shift the kernel by one column. Then we multiply the next 3x3 with kernel and add the overlapping pixel values.

for small kernel sizes (say 3x3 or 5x5) this approach maybe faster. Using fft is alluring when working with large kernel sizes. The approach is go block by block and use the overlap save/add method (you can find a lot of info on this just google it for 2D). The reason we need to use overlap method is because fft facilitates circular convolution and not linear convolution (but circular convolution can be made to behave like linear convolution by adjusting the input and kernel sizes).

Once you understand the method then it is just deciding on the kernel size and using the equation for those u,v values. Hope this helps.


Beam width of linear array using MVDR beamformer I did delay and sum beamforming technique and MVDR beamforming technique in MATLAB in time domain, but I could not make a comparison between the beam patterns as both showed same beam width. But when I did it in frequency domain, I could differentiate both the techniques as MVDR beamformer showed better resolution( low beam width) compared to delay and sum bemaformer.
I used 10 array elements, with input signal sine added with noise and the steering angle is 20 degrees. Is there something wrong with my MVDR code? Please help.

    %%
clc;
clear all;
close all;

%% Time specifications:
   Fs = 25e4;                   % samples per second
   dt = 1/Fs;                   % seconds per sample
   StopTime = 1e-3;  %0.952e-3;           % seconds
   t = (0:dt:StopTime-dt)';     % seconds

   %% Sine wave:
   Fc = 8e3;                     % hertz
   signal = sin(2*pi*Fc*t);
   x=signal+randn(size(t));
%%
N = 10;
d = 0.5;
plotangl = -90:90;
tm = plotangl * pi/180;
elementPos = (0:N-1)*d;
thetasteer=20;
c=1150;
k=exp(-1j*2*pi*elementPos'*sin(thetasteer*pi/180)/c);
x1=k*x';%10x250
Sn=(x1*x1');
a=0.01.*eye(10);
Sn1=Sn+a;

%%

as=exp(-1j*2*pi*elementPos*sin(thetasteer*pi/180));%1x10
w1=(inv(Sn1)*as');
w2=(conj(as)*inv(Sn1));
w3=w2*as';
w=w1/w3;%10x1
%w=ones(10,1);%10x1
W=w*as;%10x10
%%
vv = exp (1j *2*pi*d* (0: N-1)'* (sin (tm)));%10x181
B1=W*vv;
B2=B1'*x1;
B3=abs(B2/max(B2))';%normalize
figure(2);subplot(2,1,1);polar (plotangl*pi/180, B3);title('polar plot of mvdr ');

B = 10*log10(B3);
figure(2);
subplot(2,1,2);
plot(plotangl,B);
grid on
xlabel(' Angle (degrees)');
ylabel('Normalized Beam Power[dB]');
legend('20 deg');
title('mvdr Array Pattern, 10 array elements');

I think you're mixing time domain and frequency domain, your signal $x$ is in time and you're multiplying it by a complex exponential, generating a complex signal that is not the Fourier transform of $x$, if you want to delay $x$ in order to simulate the angle of arrival of the signal you have to do $x_{m}(n)=x(n-\tau_{m})$ (in the time domain) or $X_{m}(\omega)=e^{-j \tau_{m}}X(\omega) $ being "m" the microphone number.

Also you have to include the wave number in the steering vector $k=\dfrac{2\pi}{\lambda}$.


Using HMM for speech recognition I'm trying to understand how speech recognition works, when using HMM. I know this post is very long, but that's mainly because I'm a complete beginner, and I since I have no idea what is clear and what is not, I try to write everything as clearly as possible (which takes a lot of space). Help would be really appreciated.

I will be using O for observation sequence, X for state sequence, M for a HMM model with defined parameters.

My understanding - please correct me where I am wrong.

If I understand it correctly, the 3 basic operations with HMMs are:
1, Evaluation - how likely is observation O to happen, given a HMM model M - P(O|M). This is done by using the forward algorithm (dynamic programming).

2, Classification - I'm not sure I understand what this one is - is it finding out which state sequence is the most likely one to cause the observation O, given a model M? If yes, then it it is again done with dynamic programming, very similar to the forward algorithm, except you will be taking maximums instead of sums.

3, Training - You change the parameters of M, to raise the probability of an observation O happening with your given model. The forward-backward algorithm can be used for this, and it is proven that it finds a local maximum.

The actual speech recognition, here in the context of isolated word recognition:

For each word you want to recognize, you have a specific HMM model - here, the example I will use is that if you want to recognize 'one', 'two', 'three' you will have 3 HMM models M_1, M_2, M_3, each with say 3 states and 3 emmissions (that might be enough for 3 words?).
You make training samples of the specific words, and you train the models with these sounds - you do this, by  breaking up each utterance/sound of a certain word into small parts (~10ms), to each part you map a certain vector, with the vector being based on some spectral quality of this 10ms part - this is the observation at a certain time.

Now, if you want to recognize a word, you say the word, and you pick such a model M_i, that it is the maximum of P(O|M_i), i element of L, L being the number of words (and specific HMM models) that you have trained.

What I don't get:

What exactly can the vector/observation be? Can it be a relatively simple function, say approximation of the spectral envelope? Or is it usually a vector of R^n?

How is it, that if you say "ooone" instead of "one", P("ooone"|M_1) will still be the highest of the three trained HMM models? Is it because if the model was trained with 'one', the model already favors the observation that happens when the sound of 'o' is said, and the state transition of x_1 -> x_1 (where x_1 is likely to emit this observation) is also favored? Basically, is there some convincing proof/intuition, that pronouncing a certain word very slowly at certain parts will not cause it to be recognized incorrectly?

In the case of isolated word recognition, does the classification HMM operation happen at all? All I see is training, and evaluation (for finding the maximum of P(O|M_i).

Lastly, could some please explain how continuous speech recognition differs from this? Apparently, unlike the above, classification happens at some point - I don't really understand how it can be used. I can imagine that maybe for each phoneme, a certain HMM model will be used, and then in a long sentence, this sentence will consist of many of these models tied together.
How many HMM models will there be, how many states will there be? How will it be done, that a certain state sequence will define the words in the utterance?

What exactly can the vector/observation be? Can it be a relatively simple function, say approximation of the spectral envelope? Or is it usually a vector of R^n?

You won't go very far by summarizing the entire spectral envelope of a sound by a single real, so yes, observations are vectors - their dimensionality is in the 10 to 100 components bucket...

Typically, the features used are MFCC and their first/second order derivatives. They capture spectral envelope (in short: dimensionality reduction applied to a low-resolution spectrum computed on a perceptual scale) and have several invariance properties that make them relatively robust.

  How is it, that if you say "ooone" instead of "one", P("ooone"|M_1) will still be the highest of the three trained HMM models? Is it because if the model was trained with 'one', the model already favors the observation that happens when the sound of 'o' is said, and the state transition of x_1 -> x_1 (where x_1 is likely to emit this observation) is also favored?

Yes, this is what happens.

The cost of stretching the "o" is high, but still lower than having the emission probabilities totally wrong by matching the "o"s or the "n" that follows with other phonemes.

  Basically, is there some convincing proof/intuition, that pronouncing a certain word very slowly at certain parts will not cause it to be recognized incorrectly?

If the features used for the observation vector are robust, and if the model is trained on the same speaker as for recognition (or if suitable model adaptation measures are used), the emission probabilities can get very low when confusing phonemes.

One thing that could cause confusions would be a HMM trained on utterances of "five" and another trained on utterances of "fiiiiiiiiiiile". In this case, "fiiiiiiiiiiiive" might be recognized as "fiiiiiiiiiiiile" - the cost of repeatedly stretching the "i" might be enough to counterbalance the confusion from "v" and "l" on a handful of frames. So an abnormal utterance could be confused with another abnormal utterance for which a model is available. The thing is, abnormal utterances are "averaged out" during the training process.

  does the classification HMM operation happen at all?

For recognition, it is not relevant to know the most likely sequence of states. However, there are applications for which knowing the sequence of state is needed.

In particular, we sometimes want to exactly synchronize the recognized sequence with the original audio recording. This could be useful for subtitling, or recording a timestamped transcription of a recording. This is also used during Viterbi training, an alternative to the Baum-Welch algorithm that uses state decoding (classification) in the training loop.

  Lastly, could some please explain how continuous speech recognition differs from this?

Left-right models with a handful of states are used to describes diphones or triphones. States have hundreds to thousands mixture components (or nowadays we eschew GMMs and use other emission models). This is already a ridiculously large model space, so there are tied-states - that is to say some states in the model share the same emission probabilities with other states in the model. Anything that pools together the parameters from several parts of the model ensures that more training data will be used.

Then, the triphone models are concatenated together to build word models - using pronunciation dictionaries (lexicons). The lexicon maps a word into a sequence of phones. While there are approaches to automatically learn it, this data primarily comes from linguists.

Then word models are concatenated together to build a language model. The result is an extremely big FST, which is not even fully composed in memory from its component - but only traversed on the fly during recognition.

Of course the whole thing is not trained with Baum-Welch... We train separately each level... Tri-phone HMMs are trained on aligned speech data; lexicons are curated by linguists; language models trained separately on purely textual data...

Typically a small amount of aligned speech data (in which the transcription are aligned with the audio) is used to bootstrap the triphones model. Then this inexact model can be used to align a larger set of unaligned training data, and the process is iterated.


How to Use the Walsh Hadamard Transform to Filter 24 bit RGB Color Values in a Localized Spatial Region? After reading this Signal processing Stack Exchange article tonight, What is the Walsh-Hadamard Transform and what is it good for?, I became interested in applying the Walsh-Hadamard transform to filter one of three classes of 24 bit RGB color values, True Blue or True Green or Tree Red, in a localized spatial region of the color image seen from a commercial airplane cockpit.

In the book Color Image Processing and Applications, the authors Konstantinos Plataniotis, Anastasios N. Venetsanopoulos state that the RGB space is not sn efficient representation for compression because there is a significant correlation between the three color components since the image energy is distributed equally among them both spatially and spectrally.

A solution is to apply an orthogonal decomposition of the RGB color signals in order to compact the color image data into fewer channels. The commonly used YIQ, YUV and YCBCR color spaces are examples of color spaces based on these principles. The resulting luminance chrominance decomposition exhibits unequal  energy distribution from the luminance component in which the vast majority of fine detail high frequencies can be found

I would like to apply color image YUV or YUV (or YCbCrYCbCr which is like HSV) compression to theoretically attain almost real time filtering of RGB color coordinates in a localized spatial region of the color image seen from a commercial airplane cockpit. The reason for the spatial localization is that we do not want to subtract green from everwhere uniformly in the color image projected to the commercial airliner flight team.

Please consider the paper, NTSC component separation via Hadamard transform,
Published in: Vision, Image and Signal Processing, IEE Proceedings -  (Volume:141 ,  Issue: 1 ) Feb. 1994 page 27-32 where the authors, M. Coley and J.Barba of the New York Institute of Technology present a technique for separating the NTSC composite signal into its components (Y, I, Q) by applying an N×N Hadamard transform directly to the composite signal sampled at four times the colour subcarrier.  They show that the Y, I, Q component signals are mapped into specific areas in the Hadamard domain. Component separation is achieved by assigning particular transform coefficients to each component signal. Simulation results of component separation applied to five colour images and one black and white test pattern
are presented. Both signal to noise ratios and subjective results are presented. The results indicate that no noticeable degradation is expected in typical colour images.

  I am interested in restricting RGB color filtering to small localized
  spatial regions of the color image acquired by per pixel processing of
  an array of sensors developed by either by M.I.T Lincoln Labs [27]
  Kelly, M. W. and Blackwell, M. H., “Digital-pixel FPAs enhance
  infrared imaging capabilities,” Laser Focus World 49(1), 90 (2013).
      or Stanford University's DEpartment of Electrical Engineering and Computer Science, [28] Wan, G., Li, X., Agranov, G., Levoy, M., and
  Horowitz, M., “CMOS Image Sensors With Multi-Bucket Pixels for
  Computational Photography,” IEEE J. Solid State Circuits 47(4),
  1031–1042 (2012).

Thank you.

The JPEG-XR Standard is using the Walsh Hadamard transform (Instead of DCT / Wavelets).

It has pretty good results (Unfortunately, it didn't get enough traction, yet it is supported on Windows).

You can follow the standard and have idea how to use it for image compression.
You'll be able to even get some source code.


Differences between the mean filter, trend(mean) removal, and baseline correction filter? I've recorded acceleration signal data using an accelerometer attached to vehicle and I want to know what the differences between the

mean filter,
trend(mean) removal,
baseline correction filter, and
filter to correct displacement for oscillation about zero baseline

are.

A mean filter is one precise type of linear filters that replaces a value by a weighted combination averaging the neighboring values. A constant offset value, or an affine trend, is "quite" invariant by an averaging filter. Therefore, it may estimates a constant or linear baseline well.
Trend removal might be less well-defined. For the mean, there are many alternatives, like the exponential weighted moving average (EWMA), which is an IIR filter.
Baseline is even less  well-defined. It can dwell between simple instrumental offets and complicated "references" that may include noie.
The last one is the fuzziest to me. But there is a notion of oscillation, that could be incorporated into a filtering model with knowledge of the signal. Tracking and informed filtering (eg Kalman) can be useful for that purpose


can $\frac{1}{H(z)}$ be causal and stable? if we have linear phase FIR filter $H(z)$ which is causal and stable
can $\frac{1}{H(z)}$ be causal and stable ?

can it be causal without been stable ?

For 1/H(z) to be causal, it must not have more zeroes than poles. Therefore, H(z) must not have more poles than zeroes.

If H(z) has linear phase, then it must have zeroes either in the unitary circle or simultaneously inside and outside of it. Therefore, 1/H(z) can't be stable if H(z) has linear phase.


Can state covariance matrix in Unscented Kalman filter contain negative values? On the one hand covariance matrix can't contain negative values by definition, but as far as state covariance matrix innovated by subtraction:

$$\mathbf{P} = \bar{\mathbf P} - \mathbf K\mathbf P_{z}\mathbf K^\mathsf T$$

it's possible to get negative value.

Can state covariance matrix in Unscented Kalman filter contain negative values?

actually a covariance matrix can have negative values off the diagonal.  variables can be negatively correlated, but the diagonal terms are variances which must be positive.

The update you show can be a problem which is why updates of $P^{-1}$ are preferred.


Tutorial on 2nd generation wavelets (with lifting)? For some denoising and deconvolution experiments, I'd like to apply a 2nd generation wavelet transform (using lifting steps) to images.

I know that there are several implementations available, but most of them use matlab, while I want to work in C++ with OpenCV.
Since there is no built-in wavelet transform implementation in OpenCV 2.x, I plan to implement it myself (plus, it will make a good exercise for me).
After some research, I've been able to find the original articles about the 2nd generation transform, but I'm still a bit confused about the exact way the algorithm works.

Taking for main reference the paper [1] by Sweldens: The Lifting scheme: a construction of second generation wavelets, I'm still confused by the definition of the index sets $\mathcal{K}(j)$ : what is their size ? how are they built ? ...

Hence my question : does anyone know about some resources about the 2nd generation wavelet transform (papers, tutorials, slides...) that are either in a tutorial-like form, or that provide a more algorithmic view (rather than a mathematical one), which would help me designing my own implementation ?

Thank you in advance.

References

My main reference is:

[1] Sweldens, W. (1998). The lifting scheme: A construction of second generation wavelets. SIAM Journal on Mathematical Analysis, 29(2), 511.

And I am also learning from:

[2] Daubechies, I., & Sweldens, W. (1998). Factoring wavelet transforms into lifting steps. Journal of Fourier analysis and applications, 4(3), 247–269.

[3] Kovacevic, J., & Sweldens, W. (2000). Wavelet families of increasing order in arbitrary dimensions. Image Processing, 9(3), 480–496. doi:10.1109/83.826784

I've bought finally a copy of [Ripples in Mathematics
The Discrete Wavelet Transform][1], and I'm very pleased by this book. The authors explain the DWT with alternating point of views (lifting schemes, filter banks approach, multi-resolution analysis), where each of these viewpoints has its own advantages. Furthermore, the book is implementation oriented, with chapters about boundaries handling and matlab/C implementations.

I'm still looking on a proper way to handle odd-sized signals, but Ripples gave me a good start.

[1]: http://www.control.auc.dk/~alc/ripples.html "Ripples in Mathematics The Discrete Wavelet Transform", by Arne Jensen and Anders la Cour-Harbo


Is PSOLA or WSOLA better for vocal stretching/retuning? I'm looking to implement a time-domain pitch modification algorithm such as PSOLA (Pitch Synchronous Overlap and Add) or WSOLA (Waveform Similarity Overlap and Add).
At the moment I'm using PSOLA, but it seems to me that WSOLA would be more robust to polyphonic signals with complex waveforms, whereas PSOLA works better with monophonic signals such as vocals. I feel it could be better to future-proof by using WSOLA, even though I'm only interested in monophonic signals for now.
I'm wondering if there is any disadvantage to using WSOLA as opposed to PSOLA.
It seems both techniques work in a similar fashion, both identifying repetition within the waveform, although WSOLA does not necessarily lock to the fundamental. Does fundamental locking make PSOLA superior for vocals?
To be more precise, I'm trying to work out the advantages and disadvantages of WSOLA vs PSOLA for the following tasks:

Pitch correction
Time correction
Formant correction

To be more precise, I'm trying to work out the advantages and
disadvantages of WSOLA vs PSOLA for the following tasks:

Pitch correction:

PSOLA - Advantages you can change the pitch directly based in the rate space between  grains, keep formants intacts...  to change the pitch you not need time scale and  resample after.
PSOLA - Disadvantage do you need a super power pitch detector
WSOLA - Disadvantage WSOLA is not one Pitch Shift algorithm, WSOLA is one Time Scale algorithm, to change the pitch do will need time scale and them resample, the formants will be lost, so you just maybe need one kind of autocorrelation to find nice places to slice, so it is not the same as Pitch information, if are you building an Automatic Pitch Correction one basic autocorrelation will not give to you the correct Pitch Information.
WSOLA - one advantages here is that it is not Pitch dependent to slice, just slice in the best match similarity, remember it is just one time scale algorithm, you will need resample to change the pitch

Time correction:
For time scale I think the WSOLA have best advantage here or some TDHS based algortihm, yeahh you can use PSOLA to time scale and apply pitch shift together, just pitch shift or just time scale , but if are you just time scaling, WSOLA is one smart choice, no need build an Ultra Pitch Detector in WSOLA, just one wave form similariry will solve your problem (AutoCorrelation can solve ...)
Formant correction:

PSOLA - by default formants will be untouchable, but you can warp the formants too, it is amazing in PSOLA ...
WSOLA - after Time scale and resample your formants will be lost and to get the original formants you will need get the formants in the spectral envelope from original signal and apply it in the Pitch Shifted signal, you can do it using an low pass filter in the spectral, use a peak picking and interpolate or use a cepstrum to get the envelope ...

I know that PSOLA is infinitely more complex, but I consider this algorithm really great, if you have a decent PitchTrack it will give you great versatility


Why do frequencies of analog signals range from $-\infty$ to $\infty$ while frequencies of digital signals are restricted to $[0,2\pi]$? In Fourier analysis while dealing with discrete-time signals, frequencies range from $0$ to $2\pi$ why? Intuitively how can i understand it?

The digital frequency span of 0 to $2\pi$ is the normalized angular frequency given in units of radians per sample.  For example, if we had a frequency tone that went $0.2\pi$ radians/sample, then it would take 10 samples to complete 1 cycle. So if we sample a sine wave with the result of having 10 samples in one cycle of the sine wave, the normalized radian frequency in that case is $0.2 \pi$.

Consider an analog signal expressed as:

$x(t) = \cos(\Omega t)$

Where $\Omega$ represents the analog angular frequency $\Omega = 2\pi F$, with $F$ being the frequency in Hz.

When we sample this signal at the sampling rate $F_s$, the sampling interval is $T_s= \frac{1}{F_s}$.

Therefore the discrete time signal is given as

$x(nT_s) = \cos(\Omega n T_s) = \cos (\frac{\Omega}{F_s}n)$

Which can be described as

$x(n) = \cos(\omega n)$

with $\omega$ being the "normalized frequency" $\omega = \frac{\Omega}{F_s}= 2\pi \frac{F}{F_s}$

Similarly, we can describe the digital frequency span in cycles per sample, and in this case the unique span would extend from 0 to 1, or if we prefer -0.5 to 0.5. (In comparison to radians per sample with a unique span from 0 to $2\pi$ or $-\pi$ to $\pi$.

Here is a real example I have that will hopefully make this clearer:

The figure below of a different example shows a digital spectrum as an "unwrapped view" extending from plus or minus infinity. This demonstrates the aliases that would exist for an arbitrary spectrum that exists with a spectral density from -3 to +3 Hz when sampled at 10 Hz. This is the meaning on "Unique Span" mentioned above; the portion of the frequency spectrum that is unique (and everything beyond is just replicated).

As @Fat32 will attest, I find it much easier to visualize individual frequencies as complex tones ($e^{j\omega t})$, specifically rotating phasors on an IQ plane, rather than cosines or sines such as $\cos(\omega t)$. (Which we see through Euler's identity consists of two spinning phasors). A single tone $e^{j\omega t}$ is a phasor rotating counter-clockwise versus time on the complex plane at rate $\omega$ radians per second. Consider cases of this as a sampled system in the plot below and see how useful and straight forward the idea of normalized frequency is!

Below we have a 3 Hz analog signal represented by $e^{j2\pi f t} =e^{j6\pi t} $, with different sampling rates resulting in different normalized frequencies.

Woahh!!! What happens here??? Now we see aliasing as it really should be shown wrapping around from the positive to the negative frequency (rather than "folding" back in) but that is another topic!


Why these harmonics always present? I am using Audio Precision 2700 to test some test vectors.
I have connected resistor banks as class - d power amp's load, and tap from this resistor using a XLR cable to the AP.
When my power amp is powered up, but no song is playing(song is paused, not muted) i checked fft spectrum of the input channel , i get the following spectrum

You can see the three peaks in frequency above 4k. (Likely, the 3 peaks in frontier band upto 200 hz). Which causes this harmonics?
The test vector is playing at 48khz.
Note : There's a passive filter after the power amp's output having cut off at 30khz.
Please point me some insight. Does the passive filter aliases the input signal, so that i get undesired signal at folded frequency?[This is my initial doubt, because this is what i believe should be the cause, if it's not then some one please tell me why this couldn't be the problem]
Thanks folks in advance...

If your class-D amplifier's switching frequency is 5 kHz above a multiple of the sample rate of your audio analyzer, then the aliases of the attenuated harmonics of the switching frequency will march upwards in frequency by 5 kHz across (and around) the spectrum displayed by the analyzer, until they disappear below the noise floor of the system and analyzer.


Using Linear Algebra for DSP I am new to DSP in general, but can one use linear algebra by itself to characterize a signal?

My first idea was to transform the signal into a matrix and then use the determinate to characteristic the signal. The big problem being that matrices obey row echelon form which means the rows can be switched and produced the same determinate.

Yes indeed modern signal processing uses matrix methods of linear algebra (or Linear System Theory more correctly), in addition to the classical calculus and harmonic analysis.

Linear algebra (matrix) based DSP considers (1D) signals as vectors and transforms as matrices that map inputs to outputs. As CedronDawg also mentioned, in such a representation the DFT (discrete Fourier transform) is merely a matrix mapping:

$$ Y = \Phi X $$ where $X = [x_1,x_2,...,x_N]^T$ is an $N \times 1$ input vector, $Y = [y_1,y_2,...,y_N]^T$ is the output vector,  and $\Phi$ is the $N \times N$ DFT transform matrix.

For a brief treatment of matrix methods in DSP, have a look at the books:

Multiresolution Signal Decomposition_AKANSU
Statistical Digital Signal Processing_HAYES


Detect playback position in swept sine sweep I am wondering what the most precise method is to detect the playback sample position of a sine sweep? The real playback sample position is unknown and an DSP with constraint computing capabilities shoud be used for the estimation, recording the playback by an external microphone.

In the concrete case the sine sweep is repeatedly played with fixed duration in the ultrasound range and a microphone 50 cm away from the loudspeaker is used to record the chirp and detect the playback position. After each single chirp playback might be a delay, whereby the delay time is known and doesn't change. Real-word noise is given in the background, which might be as loud as the ultrasound playback. One period of the sine sweep looks like this:

As visible, the frequency in the chirp is either linear or exponentially increasing from a low frequency to a high frequency, fade in and fade out is applied and the chirp duration is 100ms with 44kHz sampling rate.

Attached Matlab sample code to apply an matched filter to a sine sweep for finding the playback position (the playback position estimate can be derived as the position of the peak in the matched filter plot):

% params

sampling_rate = 44000;
f_low = 10;
f_high = 100;
shift_percent = 0.0;
shift_samples = 30000;

% init
t = 0:1/sampling_rate:1;
N = length(t);
y = chirp(t,f_low,1,f_high, "linear");
y_inv = flip(y);
if(shift_percent ~= 0.0)
    y = circshift(y, round(shift_percent * N));
else
    y = circshift(y, shift_samples);
end
figure(235)
plot(y)
fade = [linspace(0, 1, (N/2)), linspace(1, 0, (N/2) + 1)];
y = y .* fade;
y_inv = y_inv .* fade;
y_matched = conv(y, y_inv);
figure(236)
plot(y_matched)

% plot
fvtool(y)


Inverse a frequency response I haven't recorded a precise impulse response, but I have carefully written the attenuation/boost of a recording chain for every third of octave:

20 Hz,   25 Hz, 32 Hz, 40 Hz, 50 Hz, 64 Hz, 80 Hz, ..., 640 Hz, 806 Hz, 1016 Hz, 1280 Hz, ...
-1 dB, -1.3 dB, -2 dB, ..., +2 dB, +2.5 dB, +1.5 dB, ...

How to inverse this frequency response? I can do it with an equalizer in any audio editor, but how would you do it with code? (e.g. Python)

How to make use of this list [-1, -1.3, -2, ..., 2, 2.5, 1.5, ...] ?

What you're most probably looking for is a FIR filter designed using a window filter design method.

Essentially, you could argue "Hey, I know how I would like my filter's amplitude vs. frequency diagram to look like" (and use the inverse of your measured gains as exactly that), and then argue "OK, I want to multiply my signal's spectrum with this equalizing spectrum, so let's instead convolve the signal with this spectrum's inverse Fourier transform". And that actually works – but you get into a lot of problems because of effects of the fact that your filter "template" essentially has limited length, and you get spectral leakage. So you add windowing to the method, and tadah, get a usable method of designing a filter.

The article of this name on DSPRelated might be of interest.


How to find Sampling frequency of a video signal i took the Digital Image processing subject in college and i am a little confused that how i can find out sampling frequency of a video signal given its bandwidth. For example what will be the sampling frequency of a video signal of bandwidth 18kHz or 15kHz.

It is not clear from your question if your "video signal" is sampled already or if you are trying to find the sampling frequency for an analog video signal.

In either case, assuming your video signal is real, the sampling frequency must be greater than twice the bandwidth of your signal. To allow for anti-alias filtering we typically sample a little higher than that which is a trade-off with the filter design. A reasonable approach is to go 20% higher, but going more relaxes the filtering requirement at the expense of a higher sampling rate.

If your signal was already sampled, and you are actually trying to determine the sampling rate, then you will need some point of reference (since by implication of the question you just have samples but do not know the time between the samples, which is the sampling rate you are trying to determine). In this case you would need to know the frequency content of your signal, and then you could take an FFT to then determine the sampling rate based on mapping your known signal to FFT bin (given the highest bin of the FFT is one sample less than the sampling rate).


Commonly used test image sets In the image analysis literature there are a variety of ways people show the usefulness of their method. For example, performance on an ideal signal corrupted by noise, an artificial signal encompassing a range of relevant structures, or real world images.

What are some commonly used test image sets for the following methods?

Line and edge detection
Corner and junction detection
Segmentation
Optical flow
English / Chinese character recognition
Handwriting author recognition
Stereo vision / image registration

Segmentation : Berkeley Segmentation Dataset and Benchmark

Stereo : Middlebury Stereo Datasets

Shape from specular flow - Shape from specular flow dataset

Also, check out this link : http://clickdamage.com/sourcecode/cv_datasets.php


Where and how do an Image and a complex number meet? I am trying to implement a Gabor Filter bank. Gabor Filter uses a complex function. How do the filters like Gabor Filter operate on an Image?

Does a complex number represent an entire Image or only a pixel?

Is an image a 2D array of complex numbers?

No it's not a complex number but a real number per pixel for a gray scale image and 3 real numbers per pixel for a standard RGB image which is supposed to be a matrix of N by M pixels, stored in an appropriate storage class in a given technology and programming language, such as an array type in C.

Gray scale digital images are represented as 2D arrays of numbers, each representing the intensity of the image at the given pixel, traditionaly in 8 bit integers between [0,255] and also in [0,1] for floating point implementations.

Color images are treated in a number of ways. Most fundamental approach is to use the RGB convention to represent each pixel with its optical components such as Red, Green and Blue colors. These values are again mostly either 8 bits in [0,255] or floating values in [0,1].

In the old days, color images were also treated via palettes that I don't want to mention here. There are also other represantations for images such as HSV,YUV,YIQ,YPbPr,YCmCr (analog and digital) codecs, commercial and broadcasting applications, but are not much related with your question.

Note that the RGB triad could also be appended with a fourth byte, referred to as the alpha channel which can be utilized by some applications

Also newer imaging standards allow more than 8-bits per color component.


Calculating bit-error-rate of QPSK data I am working on a simulation of QPSK data using the code below (Matlab):

j = sqrt(-1);
N = 10^6; % number of symbols
Es_N0_dB = 12; % ratio of signal energy to noise energy in dB
QP = (2*(rand(1,N)>0.5)-1)+j*(2*(rand(1,N)>0.5)-1); % actual QPSK data
s = (1/sqrt(2))*QP; % normalization of energy to 1
n = 1/sqrt(2)*(randn(1,N)+j*randn(1,N)); % white Gaussian noise, 0 dB variance
y = s + 10^(-Es_N0_dB/20)*n; % recovered noisy data
plot(y, 'b.')

This is my resulting constellation diagram:

To identify the errors (and calculating the bit-error-rate), I have divided the diagram into 4 quadrants and used the following to find out which elements have changed sign (i.e. crossed the decision boundaries):

sr=real(s); si=imag(s);
yr=real(y); yi=imag(y);

for i = 1:N
    if sr(i)*yr(i)<0
        errorI(i)=yr(i);
        if si(i)*yi(i)<0
            errorQ(i)=yi(i);
        end
    end
end

err = errorI+j*errorQ;
BER = numel(err)/numel(y);

However, this does not return any data for the in-phase and quadrature errors. What is the problem with this approach?

Any explanation would be greatly appreciated.

You're not doing anything conceptually wrong, you're just misunderstanding MATLAB. The code as you provided it doesn't execute properly; here's a fixed script:

% don't need this; MATLAB already knows what j is
%j = sqrt(-1);
N = 10^6; % number of symbols
Es_N0_dB = 12; % ratio of signal energy to noise energy in dB
QP = (2*(rand(1,N)>0.5)-1)+j*(2*(rand(1,N)>0.5)-1); % actual QPSK data
s = (1/sqrt(2))*QP; % normalization of energy to 1
n = 1/sqrt(2)*(randn(1,N)+j*randn(1,N)); % white Gaussian noise, 0 dB variance
y = s + 10^(-Es_N0_dB/20)*n; % recovered noisy data
plot(y, 'b.')
sr=real(s); si=imag(s);
yr=real(y); yi=imag(y);

% it seems like you want to predeclare these since you're indexing them below
errorI = zeros(size(s));
errorQ = zeros(size(s));

for i = 1:N
    if sr(i)*yr(i)<0
        errorI(i)=yr(i);
        if si(i)*yi(i)<0
            errorQ(i)=yi(i);
        end
    end
end

err = errorI+j*errorQ;
% BER is the number of nonzero entries in `err` divided by the number of chances
BER = sum(err ~= 0) / length(y);


Why is the continuous time Fourier series of DC signal an impulse? In case of continuous time Fourier transform(CTFT), I can easily calculate the Fourier transform of DC signal by using Fourier duality or inverse CTFT. But I don't know how to calculate the continuous time Fourier series(CTFS) of DC signal. Some people say "One point in the frequency domain which includes the energy of all time should be an impulse." and many textbook says "the CTFS of DC signal is an impulse." But I think the inverse CTFS of an impulse is not DC signal by my calculation. I think I should consider a periodic function that converges to DC, calculate it's fourier series and then take the limit. But all my attempts failed. Who knows correct solution to this problem?

CTFS , $a_k$, of a periodic function $x(t)$ of period $T$ is :
$$ a_k = \frac{1}{T} \int_{<T>} x(t) e^{-j \frac{2\pi}{T} k t } dt \tag{1} $$
and the DC coefficient corresponds to setting $k=0$ above;
$$ a_0 = \frac{1}{T} \int_{<T>} x(t) dt \tag{2} $$
where the integration is along one period of $x(t)$. The DC coefficient cannot be an impulse, (indeed the CTFS coefficients $a_k$ form a discrete sequence, in $k$, of finite amplitude by definition of its existance.)
A signal must have finite energy per period (aka finite avg power) in order to have a valid CTFS (due to Dirichlet conditions) representation, hence the above integral in Eq.2 should be finite; i.e., the integral converges, and cannot result in an impulse.


How to approximate these fuzzy rectangle shapes in binary maps? The following binary maps are the results of a measurement I've carried out. I'm looking for a quantifiable measure of the red areas. I could trace rectangles manually (the green rectangle on the upper-right image), but it's hard to motivate their shapes.

Is there a natural (and well-used) method to detect the shapes of these rectangles?

My initial idea was to sum all columns and all rows, and then find the edges through some thresholding. But I'd like some more ideas, because my plan feels quite ad-hoc.

You might get better answers to this question over at Math Stack Exchange, but let me give it a try.

First of all, you should try to decide what the property that you want your rectangle to maximize actually is.  I'll take a wild guess and say that it might be the number of red pixels minus the number of blue pixels inside the rectangle.

If my guess is right (or at least a good approximation), then the first thing you probably want to do is integrate your data.  That is, let $f(x,y)$ be $1$ if the pixel at $(x,y)$ is red, and $-1$ if it's blue.  Then, for each point $(x,y)$, you want to calculate

$$ F(x,y) = \sum_{i=0}^x \sum_{j=0}^y f(i,j) \;. $$

This is easy to calculate in two passes over the data, e.g. in pseudocode:

// first pass: integrate over x for each y
for y from 0 to y_max:
    sum = 0;
    for x from 0 to x_max:
        sum = sum + data[x,y];
        out[x,y] = sum;

// second pass: integrate over y for each x
for x from 0 to x_max:
    sum = 0;
    for y from 0 to y_max:
        sum = sum + out[x,y];
        out[x,y] = sum;

or, indeed, even in just one pass:

// do first row:
sum = 0;
for x from 0 to x_max:
    sum = sum + data[x,0];
    out[x,0] = sum;

// do other rows:
for y from 1 to y_max:
    sum = 0;
    for x from 0 to x_max:
        sum = sum + data[x,y];
        out[x,y] = sum + out[x,y-1];

(The careful reader will also note that both of these code examples can be made to work in-place just by replacing all instances of out with data.)

Now, if one corner of your rectangle were fixed at $(0,0)$, all you'd need to do to optimally place the other corner would be to find the maximum of $F$.  (This can be done easily and efficiently just by looping over all points and keeping track of the maximum value seen so far and its location.)

However, if both corners can be at arbitrary locations, then things get a bit more complicated.  Essentially, what you want to maximize is

$$\begin{aligned}
G(x_1,y_1; x_2,y_2)
&= \sum_{i=x_1}^{x_2} \sum_{j=y_1}^{y_2} f(i,j) \\
&= F(x_1,y_1) - F(x_1,y_2) - F(x_2,y_1) + F(x_2,y_2)
\end{aligned}$$

where the latter formula can be derived using the inclusion–exclusion principle.

Unfortunately, while the analogous problem in one dimension is relatively easy to solve, in two dimensions it gets trickier.  If you don't have too many points, you could just loop over every possible pair of points and select the one that maximizes $G$, but the time that takes scales with the square of the number of points (and thus with the fourth power of the data resolution).

Instead, you may be better off exploiting the relative smoothness of the integrated data by first sampling $G$ at a lower resolution to find an approximate maximum, and then refining it, either just by sampling the vicinity of the approximate maximum more finely, or perhaps using simple hill cimbing.  Or you could try starting a simple hill climbing process from each low-resolution sample point to find the nearby local maximum, and then take the highest of these maxima.

(You could even consider more advanced optimization techniques, such as simulated annealing, but these may well be overkill for your purposes, and might even be less efficient than the simpler methods described above.)