Periodic Signal Fourier Transform If I have a signal x(t) as follow

x(t)=$\sum\limits_{k=-∞}^∞ R_{T}(t-2kT)$

$R_{T}(t) = \begin{cases}1 ;&[\dfrac{-T}{2},\dfrac{T}{2}] \\ 0; & \text{otherwise} \end{cases}$

$\mathfrak{F}\{R_{T}(t)\}=Tsinc(\pi fT)$

$H(f)= \begin{cases}1 ;&|f| \le \dfrac{3}{4T} \\ 0; & |f| \ge \dfrac{3}{4T} \end{cases} $

$Y(f)=\dfrac{1}{2}\delta(t) + \dfrac{1}{\pi}\delta(f-\dfrac{1}{2T}) +\dfrac{1}{\pi}\delta(f-\dfrac{1}{2T})$

Can Some one Please Help me in how X(f) will Look Like graphically . I am just Poor in the graphical representation of how X(f) multiplied With H(f) to obtain Y(f).
Please this is not a home work question but i really need to know to move forward.i tried a lot but i just don't get graphical representation in frequency domain.
thank You.

The Fourier transform of $x(t)$ can be easily obtained if you realize that $x(t)$ can be written as the convolution of the pulse $R_T(t)$ and a Dirac comb:

$$x(t)=R_T(t)*\sum_{k=-\infty}^{\infty}\delta(t-2kT)\tag{1}$$

Since convolution corresponds to multiplication in the frequency domain, you get for the Fourier transform of $x(t)$

$$X(f)=\mathcal{F}\{R_T(t)\}\cdot\frac{1}{2T}\sum_{k=-\infty}^{\infty}\delta(f-k/2T)\tag{2}$$

because the Fourier transform of a Dirac comb is a Dirac comb (see e.g. this table). And since

$$f(t)\delta(t)=f(0)\delta(t)$$

Eq. (2) can be written as

$$X(f)=\frac{1}{2T}\sum_{k=-\infty}^{\infty}\hat{R_T}\left(\frac{k}{2T}\right)
\delta(f-k/2T)\tag{3}$$

where $\hat{R_T}(f)$ is the Fourier transform of $R_T(t)$. Eq. (3) shows that the Fourier transform of $x(t)$ is a weighted sum of Dirac impulses, where the weights are equidistant samples of the Fourier transform of the pulse $R_T(t)$.

Since the low pass filter $H(f)$ only passes frequencies lower than $3/4T$, only the contributions for $k\in\{-1,0,1\}$ of the sum in (3) are passed. This gives the final result for $Y(f)$. Note that there are two typos in your formulation of $Y(f)$. It should be

$$Y(f)=\frac12\delta(f)+\frac{1}{\pi}\delta(f-1/2T)+\frac{1}{\pi}\delta(f+1/2T)\tag{4}$$


use of complex number in digital signal processing how  complex number used in digital signal processing.
detailed information about on how this method will be useful in radar technology.

Any complex array could be represented as a real array as follows:
$$x_R = \begin{bmatrix}
 \Re(x) \\
\Im(x)
\end{bmatrix}$$
and matrices are represented as
$$A_R = \begin{bmatrix}
 \Re(A) & \Im(A) \\
-\Im(A) & \Re(A)
\end{bmatrix}$$
where $\Re(x),\Im(x)$ are vectors containing the real and imaginary part of complex vector $x$, respectively. Note that addition/subtraction is preserved by simply addition
$$x_R + y_R = z_R$$
You can also verify that matrix multiplications in the real domain also preserves the real/complex information in their respective blocks, namely
$$A_RB_R = C_R$$
where $C_R = \begin{bmatrix}
 \Re(C) & \Im(C) \\
-\Im(C) & \Re(C)
\end{bmatrix}$


Narrowband DOA Beamforming techniques on Wideband signal I'm trying to implement a tracking system (undergrad thesis) using DOA of acoustic signals. I'm testing this by playing a 2kHz sound and receiving the signal with four microphones in a ULA.

So, the band I'm interested in is narrowband but naturally the signal received is wideband.

I have been able to successfully implement a time delay beamformer however when I have tried to implement a phase shift beamformer I have failed.

My suspicion is that you can only perform phase shift beamforming techniques on narrowband signals and will not work even if you are looking for a narrowband within a wideband signal. Am I correct?

If so, will a simple filter be good enough to correct this issue to transform the signal into a narrowband one?

My suspicion is that you can only perform phase shift beamforming techniques on narrowband signals and will not work even if you are looking for a narrowband within a wideband signal. Am I correct?

No; the correlation coefficient of your reference signal to your received signal should, given additive, phase-uniform, white noise, have an expectation value of the phase shift the individual receiver sees, no matter how wide your band is.

Note that this is only true if your transmitter/receiver system can be modeled as linear phase for the bandwidth we're talking about, but without loss of generality, radar/sonar/sound systems are designed to pretty much have exactly that characterestic.


Help me understand this simple code snippet to get the magnitude frequency response I inherited the following FFT code to get the magnitude frequency response of a sample window. I understand what fft() does. What I don't understand are the lines marked with a "???". What do they do? In particular the "real[i]/real.length/2" portion and the multiplication and division by 1000 portion.

double[] real = new double[_FFTSize];
double[] imaginary = new double[_FFTSize];

loadDataIntoReal(real);

fft(real, imaginary); // in-place FFT

for (int i = 0; i < spectrum.length; i++)
{
    double reValue =((int)((real[i]/(real.length/2))*1000)/1000.0); // ???
    double imValue =((int)((imaginary[i]/(real.length/2))*1000)/1000.0); // ???
    double value = Math.sqrt(reValue * reValue + imValue * imValue);
    spectrum[i] = value;
}

Each line seems to have a redundant multiplication and division by 1000. However, if I change the code to be the following, then I get slightly different floating-point results, so I guess they must provide some sort of precision.

    double reValue =(real[i]/(real.length/2));
    double imValue =(imaginary[i]/(real.length/2));

It looks like this mul-floor-div is an incorrect attempt to limit some numeric print-out to 3 or 6 decimal places.  Better to do that during the print formatting.


Phase of the chirp signal I'm having a tough time understanding the phase increments which are quadratic in nature in the chirp signal. When I see consecutive samples, what will I see in terms of phase and if I have to measure the phase difference between the two samples, what will I see ?

Take the chirp signal

$$
x(t) = e^{j\pi\alpha t^2} = e^{j\phi(t)}
$$

Where $\alpha$ is the chirp rate of the signal. You can see that the function describing the phase, $\phi(t)$, is in the form of a quadratic equation. So you can plot

$$
\phi(t) = \pi\alpha t^2
$$

and see the quadratic nature of the phase at any point in time. The actual phase-change between each sample depends on what your sample rate is.


Clarification regarding signal after pulse shaping As can be seen from the figure below, I have plotted a signal before it is pulse shaped and after performing half sine pulse shaping using the formula

$p(t) = \sin(\pi t/2T_c),\ 0<t<2T_c$

$p(t) = 0, {\rm otherwise}$

Now I am confused about how to interpret the signal after pulse shaping. For example, between 0 and 0.002, the digital signal (the one before pulse shaping) has  5 peaks (including the peaks at 0 and 0.002). Shouldn't I expect the analog signal to also follow such pattern with 5 peaks or is this something special to the half sine pulse shaping?

If it matters, I am using OQPSK modulation with DSSS to produce the digital signal and then performing half sine pulse shaping.

I assume in the following that your digital signal is a sequence of width-1 spikes or, equivalently, a sequence of delta pulses.

Your half-sine pulse shaping function gives you a lagged peak for every spike. Consider just your first digital signal spike at $t = 0$. If that were all you have, then the shaped output signal would be
$$g(t) := f(t) \star p(t) = \sin(\pi t / 2T_c),$$
which has a peak at $t = T_c$. That is, the original signal had a peak at $t = 0$ and the shaped signal has a peak at $T_c$. Is this the behavior you intended?

If you only had your first two peaks, then the shaped signal would be the sum
$$ g(t) = \sin(\pi t / 2T_c) + \sin(\pi t / 2T_c - T_c).$$ We can use standard trig identities to rewrite this. Set $\theta = \pi t / 2T_c$ and $\phi = \theta - T_c$, we have $\theta + \phi = 2\theta - T_c$ and $\theta - \phi = T_c$, so we can write
$$\sin(\theta) + \sin(\phi) = 2\sin(\frac{\theta + \phi}{2}) \cos(\frac{\theta - \phi}{2}) = 2\sin(\theta - T_c/2)\cos(T_c/2).$$
This function has only one peak, at $\theta = \pi/2 + T_c/2$. In other words, two peaks in the digital signal correspond to only one peak in the shaped signal. I think this is the behavior which surprised you.

Generally speaking, and without doing the math, I'd expect $N$ (edit) consecutive same-sign peaks in the digital signal to give rise to $N-1$ peaks in the shaped signal.

edit: Without doing the math, I'm not very sure what I expect out of a positive peak followed by a negative peak. I expect you'll get a positive shaped peak followed by a negative shaped peak.


Meaning of channel capacity in the AWGN channel I understand the concept of channel capacity as the maximal rate of the channel code I can apply without making a mistake in the receiver, in that sense the capacity is between 0 and 1.
What I don't understand is the meaning of the capacity in the AWGN channel case where it is calculated by C=(1/2)*log2(1+SNR) where clearly I get a number greater than 1 when the SNR is greater than 3 (linear scale).
How can I pass more than 1 bits/channel use when the SNR is very high?

Thanks

There are two different kinds of channel models that are being confused here.

In the binary symmetric channel, the inputs and the outputs
are constrained to be $0$ or $1$ and the key parameter is the
transition probability: the probability that an input $0$ is changed
to an output $1$ or vice versa.  The channel capacity $C$
is a number between $0$ and $1$ and the meaning is that by using
a suitable code of rate $R < C$ on the channel, we can transmit
at a rate of $R$ information bits per channel use with arbitrarily
high reliability (that is, with arbitrarily small bit error
probability.  A (block) code for such a channel would be a set
of $M$ binary $n$-tuples that will convey $\log_2 M$ bits in $n$ channel
uses at a rate $R = \frac{\log_2 M}{n}$. Usually $M$ is chosen as
a power of $2$, say $2^k$ for some $k < n$.  Note also that the
binary symmetric channel capacity cannot exceed 1 bit per channel use.

The discrete-time Gaussian channel is a different beast. The inputs
can be arbitrary real numbers (it is not necessary that they
be restricted to $0$ or $1$) and if the $i$-th input is $x_i$,
then the $i$-th output is $x_i+n_i$ where the $n_i$ are independent
zero-mean Gaussian random variables with variance $\sigma^2$.
$x_i^2$ is called the power in the $i$-th input and
$\sigma^2$ is called the noise power.  A (block) code for such
a channel is a collection of $M$ $n$-vectors $\mathbf x^{(i)}$
(called a codebook) that can be used
to transmit $\log_2 M$ bits in $n$ channel uses.  Various realistic
constraints are placed on the code.  One
such constraint is that the peak power is bounded by $P$:
the magnitudes of all the $M\times n$ entries $x_j^{(i)}$
in the codebook are no larger than $\sqrt{P}$.
Another is that
that the average power is bounded by $P$:
$$\frac{1}{Mn}\sum_{i} ||\mathbf x^{(i)}||^2
= \frac{1}{Mn} \sum_i \sum_j \left|x_j^{(i)}\right|^2\leq P.$$
The capacity of the channel depends on the signal to noise ratio
$P/\sigma^2$ and it can exceed one bit per channel use.  Think of it
this way: we can use multilevel pulse amplitude modulation to transmit
more than one bit per channel symbol on this channel whereas with
the binary symmetric channel we are restricted to binary modulation
and so cannot exceed one bit per channel use.

For lots of details on this kind of stuff at a level lower than
a typical graduate text, see the following links (lecture notes
of mine).

Lecture 6,
Lecture 7,
Lecture 8,
Lecture 9,
Lecture 10,
Lecture 11,
Lecture 12


comparison between frequency offset estimators I have been working frequency offset estimation in OFDM. The objective was to compare different frequency offset estimation techniques. By using MATLAB, I have simulated three different estimation techniques, namely:

By using Pilot tones
Using phase difference between two repetitive symbols (proposed by Moose)
and by using a blind method (Cyclic Prefix)
The comparison between these methods is done by simulating the mean square error (MSE), and below is the output figure

However, my supervisor highly recommended to do the comparison between the CFO estimation methods using another parameter rather than the MSE(without suggesting another way). Therefore, I would like to request from the respected members to suggest me alternatives that I can do the comparison.

I hope that I'm not violating the rules of this forum by asking this general question.

Well, since you simulate you could calculate the residual frequency offset which is basically the remaining CFO after correction or the difference between the true CFO and the estimated CFO.

Another thing that is often done is to give the resulting Bit Error Rate with different estimators.

But really, please ask your supervisor what he has in mind. Everything else is guessing and probably wasted time.


Transfer function from poles and zeros If we know that a filter (or a system) has two poles at 20 GHz and one zero at 15 GHz, then how do you write the transfer function $H(s)$ for such a system?

I am wondering why sometimes the poles and zeros are given in Hz and other times as complex numbers?

You must remember that $s=j\omega$ and that $\omega=2\pi f$ and the transfer function is:

$H(s)={\prod_n (s-z_n)}/{\prod_m(s-p_m)}$

Where $z$ and $p$ are the zeros and poles respectively. Then you just need to adapt to your case.


Find h[n] using DTFT properties Using the DTFT property, find h[n] of a system where:

Is it an FIR or IIR system?

While this is by your admission homework (and fairly basic), I'll bite. Recall the definition of the DTFT:

$$
X(\omega) = \sum_{n=0}^{\infty}x[n] e^{-j\omega n}
$$

And recall the definition of the frequency response $H(\omega)$:

$$
H(\omega) = \frac{Y(\omega)}{X(\omega)}
$$

where $x[n]$ is the input to the system and $y[n]$ is its output. Combine these two equations:

$$
\begin{eqnarray*}
H(\omega)X(\omega) &=& Y(\omega) \\
\frac{1 - a^4 e^{-j 4 \omega}}{1 - a^4 e^{-j \omega}}X(\omega) &=& Y(\omega) \\
(1 - a^4 e^{-j 4 \omega})X(\omega) &=& (1 - a^4 e^{-j \omega})Y(\omega) \\
X(\omega) - a^4 e^{-j 4 \omega} X(\omega) &=& Y(\omega) - a^4 e^{-j \omega} Y(\omega)
\end{eqnarray*}
$$

Now, perform the inverse DTFT on both sides of the equation. By definition, $X(\omega)$ and $x[n]$ are a transform pair; likewise for $Y(\omega)$ and $y[n]$. For the other two terms, recall the time-shifting property of the DTFT:

$$
x[n-k] \leftrightarrow e^{-jk\omega} X(\omega)
$$

which can be shown easily from the definition of the DFT. Using this property, the equation inverse transforms to the difference equation specification for the system:

$$
x[n] - a^4 x[n-4] = y[n] - a^4 y[n-1]
$$

$$
y[n] = x[n] - a^4 x[n-4] + a^4 y[n-1]
$$

This is the definition of a recursive filter, which are usually IIR; that is the case for this one. Finding the impulse response is easy; let $x[n] = \delta[n]$ and find that the system output is:

$$
y[n] = a^{4n}u[n] - a^{4(n-4)+4}u[n-4]
$$

$$
y[n] = a^{4n}u[n] - a^{4n-12}u[n-4]
$$

The above is plotted for $a=0.99$. It should be noted that the system is only stable for $|a|\le1$.


Why we need fourier transform of periodic signal although we have fourier series for periodic signal? I was going through some of the basics of fourier series and fourier transform. And I came across  one topic "Fourier transform of periodic signal". But I am not able to understnad why we need to go for fourier transform because we already have fourier series to take care of periodic signal.

Any help is appreciated.

If you take the Fourier transform of a periodic signal, you'll arrive at the same result that you would if you had computed the signal's Fourier series. Its spectrum will consist of discrete lines that correspond directly with the terms in the Fourier series. As an example, note the following Fourier transform pair:

$$
\cos(\omega_0 t) \Leftrightarrow \pi \left(\delta(\omega - \omega_0) + \delta(\omega+\omega_0)\right)
$$

If you computed the complex Fourier series of $\cos(\omega_0 t)$, you would find exactly two nonzero coefficients: those that correspond to frequencies $\omega_0$ and $-\omega_0$.


ROC of this LTI system given $x[n]$ and $y[n]$ So I have a system with the following inputs and outputs:

\begin{align}
x[n]&=\left( \frac12 \right)^{n}u[n] + 2^{n}u[-n-1]\\
y[n]&=6\left( \frac12 \right)^{n}u[n] - 6\left( \frac34 \right)^{n}u[n]
\end{align}

The transfer function $H(z)$ can be found to be:

$$\frac{Y(z)}{X(z)}=\frac{1-2z^{-1}}{1-\frac{3}{4}z^{-1}}$$

with a zero at $2$, and a pole at $\frac34$. Now how can we find the ROC of this system given that one of the sequences ($x[n]$) has an anti-causal term?

Let's look at the pole locations and the corresponding regions of convergence (ROCs) of the $\mathcal{Z}$-transforms of the input and output signals:

$X(z)$: poles at $z=\frac12$ and $z=2$ with ROC $\frac12<|z|<2$ (two-sided).
$Y(z)$: poles at $z=\frac12$ and $z=\frac34$ with ROC $|z|>\frac34$ (causal).

Now we know that the ROC of the multiplication of two $\mathcal{Z}$-transforms equals the intersection of the ROCs of the individual $\mathcal{Z}$-transforms, unless there is pole-zero cancellation, in which case the resulting ROC can be bigger than the intersection.

$H(z)$ has a zero at $z=2$ and a pole at $z=\frac34$. Since the zero cancels the pole of $X(z)$, the restriction $|z|<2$ for the ROC of $Y(z)$ is removed. The two theoretically possible ROCs of $H(z)$ are $|z|<\frac34$, corresponding to an anti-causal and unstable system, and $|z|>\frac34$, corresponding to a causal and stable system. In the first case, the ROC of $Y(z)$ would need to be $\frac12<|z|<\frac34$, which is not the case, and which would correspond to a two-sided unstable signal. $^1$ Consequently, the ROC of $H(z)$ must be $|z|>\frac34$, which means that the system described by $H(z)$ is causal and stable.

1. I use the term unstable signal for a signal that, if interpreted as the impulse response of an LTI system, would correspond to an unstable system (i.e., the ROC does not contain the unit circle).


Is there a hypothetical situation where once can use Wiener filter prior Matched filter to get better result? Is there a hypothetical situation where once can use Wiener filter prior Matched filter to get better result? Mostly I am looking for an example or a situation where both the filters are required as ther are doing different thing to the received signal. Under what circumstanctes this would be helpful?

A matched filter is matched to a noiseless known waveform.

A Weiner Filter is based on the correlation structure of a random waveform.

The signal models are different but often a signal has deterministic and random characteristics.

It would be a stretch to say that a Kalman Filter could reduce to a matched filter, but it does work in models without process noise.

The short answer is that one can mix detection and estimation together but the likelihood function is the better way to think about it than terms like “matched” or “Weiner” filter.


Remove Checkered Noise (Pattern) from an Image I'm trying to remove what looks like a checkered noise pattern that has been added to an image. Comparing the 2D FFTs

The one on the left is the original, while the middle one is the corrupted one.
I've been trying to remove the added bright spots on the diagonals. The far right FFT is my attempt at isolating the diagonals, which I did with an X shaped filter.

[f1,f2] = freqspace(300,'meshgrid');
Hd = zeros(300);
d = f1 + f2 == 2 * f1;
f = f1 + f2 == 0;
d = f | d;
d = logical(d);
Hd(d) = 0.75;

h = fsamp2(Hd);

However I seem to be unable to isolate the features. Are there any known algorithms for such a process? Isolating directions or points on an FFT.
The MSE of the middle image is 0.0084 and the lowest I have gotten it to is 0.0031. With a combination of a low pass filter and a median filter. However the image was quite blurry and I don't think it would be possible to get it very close to the original after that much information loss.
This is a HW question, so please only give hints.
Here are the two images, noisy and original in that order

You could model the FFT probabilistically and use your prior knowledge about the FFTs of natural images, such as the one of the left, to model a prior distribution with which to form a MAP estimate.


How to convert magnitude to dB in matlab I exported a diagram from FDTD simulation, and now I want to change the magnitude to dB, here is the exported codes:
ax1 = axes;

plot(lum.x0,lum.y0)

set(ax1, 'XLim', [1.49806e-06 1.60194e-06])

set(ax1, 'YLim', [-0.0232806 1.22328])

set(ax1,'XGrid', 'on')

set(ax1,'YGrid', 'on')

and the ploted spectrum is shown below

I just wondering how to convert the magnitude to dB?
I know there is a function called mag2dB, but i just dont know how to use it.
any help would be very appreciate.


What are good estimates for phonemes length? I am working in a phoneme recognition project, the language that I am using is "English". I have a tool that converts each written phrase in phonemes. Also I have an estimation of the time length of each word. I would like to know also the duration of every single phoneme. Does any one know if there is a good source where I can resort to know the duration of the phonemes in the English language?

NLTK has an interface/downloader for a sample of the TIMIT phoneme corpus (http://nltk.sourceforge.net/doc/api/nltk.corpus.reader.timit-module.html), and it should be straightforward to get the lengths of the phoneme recordings in that corpus.


Envelope behavior dIfference between QPSK, OQPSK and $\pi/4$-QPSK I am trying to analyze whats the major difference between these three variants of 4PSK.
I know that
QPSK -  signal points are at 0,90,180,270
OQPSK - The bit period of Quadrature component is shifted by one bit period
pi/4-QPSK -  Signal points are at 45,135,225,315

Some website says that its to avoid abrupt phase transitions through zero.. I cannot understand and figure it out the actual reason behind it.. Can anyone please point it to me

Let's look at each in turn:
Standard QPSK:
With standard QPSK, each of the signal points are in quadrature (Note that the signal points in the constellation can be at any arbitrary phase really; either 0, 90, 180, 270 OR 45, 135, 225, 315...or any phase offset as long as the four constellation points are always in quadrature-- your interpretation of $\pi/4$ QPSK is actually not correct but I will clarify that later below).
So to be QPSK each symbol is one of four states and they are in quadrature, and any state can transition to any of the other three on each symbol boundary.
The diagram below shows both the I and Q magnitudes versus time along with the same plot on the complex IQ plane as the constellation with state transitions. We can easily see from the constellation the state changes that would cause the signal to go through the origin.

This following plot from this site image link is an actual QPSK constellation with the effect of pulse shaping which controls the transition from one state to another (see more details on pulse shaping at bottom of this post). Observe the similarity to the simplified state transition plot shown above.

OQPSK
In offset QPSK, we completely eliminate the possibility of the signal going through the origin by shifting Q half a symbol relative to I. Since each symbol represents 2 bits, we are shifting 1 bit in the process.

$\pi/4$ QPSK
In $\pi/4$ QPSK we achieve the same objective (avoid the origin!) by rotating every other symbol by 45°. Note in the transitions from any one symbol to the next, as shown in the diagram below, that it is not possible to go through the origin.

Below shows the resulting constellation from any state to any other state and we again see that the origin is never passed through.

Why do we care?
When the signal is band-limited (as is nearly always the case), the envelope of the signal, and notably it's amplitude, will follow the trajectory more slowly in time as it transitions from one state to another. We implement pulse shaping to limit the bandwidth as depicted in the figure below: the upper portion of the plot shows a blue trajectory that is similar to the I and Q versus time plots above in that it is rectangular in shape, while the red shows the same trajectories after pulse shaping. Immediately below that is the corresponding spectrums also in blue and red - notice the significant reduction in occupied bandwidth! (Note that these trajectories shown in this plot are for the I or Q of a QAM waveform, hence the multiple levels versus time, but the point is to show the comparison and effect of the sharp or instantaenous transitions versus the slower transitions with pulse shaping.)

In the transmitter the modulated signal will be passed through a power amplifier, which uses a significant portion of the DC power budget for the device. With handheld mobile devices, power efficiency is a significant concern for long battery life. If we drive a power amplifier into saturation it will run much more efficiently (more transmit power vs DC power needed), but this will eliminate all the careful pulse shaping we did! The effect of the PA saturation is appropriately called "spectral regrowth". At an extreme case, consider the careful processing to implement pulse shaping on QPSK to restrict our transmitted spectrum by converting rectangular pulsed waveforms to slower transitions, then put it through a saturated power amplifier that acts as a hard limiter (in the extreme case) which serves to convert our slow transitions back to a rectangular waveform; and in the process "regrow" the transmitted spectrum back to a Sinc function. For this reason, when power efficiency is a primary concern, we are interested in either using constant-envelope modulations (such as GMSK and other CPM modulations), or waveforms with reduced peak to average power such as OQPSK and $\pi/4$ QPSK. A waveform that transitions (slowly) through the origin would have the worst peak to average power (the envelope and therefore instantaneous power goes to 0!), and OQPSK and $\pi/4$ QPSK significantly reduce this by reducing the total variation in the envelope of the waveform. Further, notice the expanded amplitude range in the continuous time domain waveform due to the pulse shaping in the plots above; the tighter we constrain the bandwidth through pulse shaping, the more overshoot we get which serves to increase the peak-average power in the waveform further restricting our transmit power for a given amplifier (this is a consideration as to how much we can restrict the bandwidth and presents a trade-space). Besides spectral regrowth specifically, power amplifiers when driven toward saturation (and at saturation run more efficiently: more transmit power for a given DC power provided to the amplifier) create AM/AM and AM/PM conversion. This specifically is distortion mechanisms from AM (Amplitude Modulation) in the signal. AM/AM is amplitude distortion due to changes in the signal amplitude (AM), and AM/PM is phase distortion (PM) due to changes in the signal amplitude (AM). By reducing the overall AM content (which is what OQPSK and $\pi/4$ QPSK do!) we can therefore drive power amplifiers further into saturation before the non-linearity effects are not acceptable.
Note too that if the QPSK waveform was not filtered, it IS a constant envelope modulation and none of this would apply as we could simply transmit using a hard saturated (highly efficient) power amplifier and not have any concern for pulse-shaping. The transition diagram shown above for QPSK would transition instantly from one state to another and therefore would go through the origin but in zero time, so there would be no AM modulation to speak of. However the obvious problem is the spectrum would appear as in the blue spectrum in the plot above, which for over the air broadcasts is not acceptable due to regulatory restrictions requiring efficient use of spectrum.


How can i measure soundcard latency? I want to write sonar/radar like sample application. For sake of simplicity (and not to use Arduino and specialized RTOS) - I simply used using audio in/out analog ports of a motherboard integrated card as a sample hardware setup + Windows API with MSVS 2015. So, for my latency tests, I just used jack-jack cable from output(green) to line in (blue) minijacks. And, in my application, I have 2 important callbacks:

void on_playback_finished(short* buffer, int length) {
    for(int i=0;i<length;i++)
        buffer[i]=(i<(length/16))?(short)(sin(i)*32767):0
}

It's quite simple, it produces short "spike" (1/16th of the length of the buffer) of sine signal of some frequency, then makes a long pause of silence for the rest of the buffer(resting 15/16 of it).

And then second callback - which is called just when I recorded buffer from the line in. And this is my question. I must somehow detect this spike and measure latency:

void on_recorded_buffer(short* buffer, int length) {
    // there is our spike, somewhere in the buffer[] samples array
    // how to detect beginning of that spike,
    // considering induced noise, attenuation and possible echo?
}

what to do? need assistance)

UPD. I am not a good drawer, but this could help to understand my question:

Red lines represent software<->hardware transitions.
The sampling rate is 192kHz, Signal length is 256 samples (1.33333ms), pause between signals is at least 768 samples (4ms). I am asking about proper algorithm of doing sample offset calculation!

Here are some approaches to consider:

FMCW:

Transmit a "chirp" signal where the frequency transmitted is slowly ramped from a low frequency to a high frequency (repeatably as a sawtooth or ramp pattern).

Multiply the received chirp with a replica of the transmitted chirp. The frequency of the resulting signal (after low pass filtering) will be proportional to the delay. An FFT is often used at this point to measure the frequency, but other frequency estimation techniques can also be used.

Review more information on "FMCW" radar for more details on this approach. Such as at this link: https://en.wikipedia.org/wiki/Continuous-wave_radar

Ranging with PRN sequences

Transmit a pseudo-random noise sequence that can be generated with linear feedback shift registers as a modulated signal in your audio band (FSK or BPSK are good choices).

Cross correlate the received sequence with the transmitted sequence efficiently with FFT's using xcorr = ifft(fft(rx)*fft(tx)), here assuming both tx and rx are real.

Note if resources are severely constrained but you have time on your side, you can perform the cross correlation over time: To demonstrate this I will give an example:

Given you want to measure 1 cm to 1 m, means you want to measure delays from 300 us to 3ms, this means we would want to use an LFSR that does not repeat from 3ms. Further for an accuracy of 1cm, a chip rate of 4KHz should be adequate. This means an 11 chip barker code such as the one I list below would work but you could also explore longer LFSR sequences.

11 bit barker code: 1 1 -1 1 -1 -1 -1 1 1 1

Note, the unique property of this code like other LFSR maximum length sequences, the cross correlation is maximum when the code is aligned (in this case = 11) and -1 every where else.

Using the above barker code with a 8Khz carrier, When the code is 1 transmit the 8KHz carrier in-phase and when it is -1 transmit it 180 degrees out of phase, thus you are mulitplying your carrier by plus or minus 1. (Note that an XOR gate is such a multiplier if you want to make this "ultra-low" resources). Because of the relationship I used between carrier and code, there will only be two cycles of the carrier for every code chip transmitted.  The sequence is transmitted continuously.

At the receiver multiply in quadrature the received sequence with the same 8KHz carrier you already generated- meaning two multipliers where one of the received sequence is multiplied by the carrier directly and the other is multiplied by the carrier shifted 90 degrees at 8KHz, a quarter cycle.  (again you could do this most simply with square waves and XOR gates by hard limiting the received signal and applying a DC offset to be midrange in the XOR gate, which would really be an "analog" solution as we would be using the XOR gate as an analog multiplier similar to a mixer- given you are using a sound card which will essentially be your A/D converter a digital multiplier will be more likely). The reason for the two multipliers is due to various phase delays, one path will be at a minimum when the other is maximum but with the two you will always have a signal to detect.

The multiplier output is followed with a simple moving average low pass filter, averaging over 1/16KHz = 62.5 us would be a good choice as it will have a null at 16KHz where we expect content of twice the carrier after the multiplication.

After the low pass filter, pass the signal through an FIR filter with the coefficients of the filter equal to your barker code sequence in reverse order. I recommend decimating the received sequence down to twice the chip rate (8 KHz) and then therefore using 2 samples/chip as the coefficients:

1 1 1 1 -1 -1 1 1  -1 -1 -1 -1 -1 -1 1 1 1 1 1 1

The output of the filter will be similar to my plot below that is specific to GPS (1023 length sequence), but to show that you will get spikes that repeat at the symbol rate:

If you compare the spike locations with a reference pulse created by your start of your sequence, you have a measure of delay. To do this, start a counter at the start of your sequence, and have a threshold detector on the filter output which stops the counter. Note as I described earlier that two multipliers and thus two filters are used and one output is "I" and the other filter output is "Q". A single estimate of magnitude can be determined using the following simple magnitude estimator:

$$|I+jQ| \approx Max(I,Q)+ 1/2 Min(I,Q)$$


Audio distortion while implementing a low-pass filter I attempted to implement a real-time audio filtering in a callback function using the simplest moving-average filter with coefficients [0.5 0.5] :

void process(void* in_, void* out_, DWORD length, void* user) {

    //data is float, mono
    int n = length/sizeof(float);
    float* in = (float*)in_;
    float* out = (float*)out_;
    float* buff = (float*)user;

    //buff holds the last input sample from the previous block
    //if this is a first block, it is 0
    out[0] = 0.5*in[0] + 0.5*buff;

    for (int i = 1; i < n; ++i)
    {
         out[i] = 0.5*in[i] + 0.5*in[i - 1];
    }

    buff[0] = in[n-1];
}

however, i get a distortion in the output, like an aliasing of some sort or clipping. But, to my knowledge none of this is possible to occur in this case. I've implemented the same code in Matlab ( although not in real-time ) and everything was fine.

I am using a BASS library.

I would like to know if i am missing something obvious here, cause i am running out of ideas. Thanks in advance!

EDIT: Thanks everyone for the answers. The above function is in the form for simplicity, but in actual from it looks like below (there is no output buffer, filtering must be done in-place) but it does exact same thing, although I will write it for convenience:

void CALLBACK filterDSPMono(HDSP handle, DWORD channel, void* buffer, DWORD length, void* user)
{
    float* x = (float*) buffer;
    int n = length / sizeof(float);
    float* buff = (float*) user;
    float temp;

    int i = 0;
    for (; i < n; ++i)
    {
        temp = x[i];
        x[i] = 0.5*x[i] + 0.5*buff[0];
        buff[0] = temp;
    }
}

EDIT2: Ok, so it turns out the problem was with the buffer size set in the library that my sound card drivers didn't like for some reason. After I changed that, everything is fine. Thank you for all the answers unfortunately I cant mark any as The Answer due to the fault being on my part. This question now can be considered as resolved.

I suspect that most likely the audio coming in is integer, and not floating point. Depending on your device, operating system and type of file input, you can get one or the other. Read documentation for audio input API you're using to see which format the audio is in. Also, if it's integer, it's likely to be 16 bit rather than 32 bit. In any case, again, read the documentation.

EDIT:
BASS library has a flag you can set to get floating point values everywhere:

BASS_SetConfig(
    BASS_CONFIG_FLOATDSP,
    BOOL floatdsp
);


Suggesting method for removing noise and image restoration? Suppose that we have degraded an image with the following:

$$g(n_1,n_2) = f(n_1,n_2)^{v(n_1,n_2)}\text,$$

with $v(n_1,n_2)$ being random noise which is independent from image and we have these relations:

$$f(n_1,n_2) > 1 ,\quad v(n_1,n_2) > 0\text.$$

Suggest the method for removing noise with transforming $g$ and filter this image in the transformed domain.

My suggestion:

I think that I take natural logarithm from equation then I have:

$$\begin{align}
\ln(g(n_1,n_2)) &= \ln\left(f(n_1,n_2)^{v(n_1,n_2)}\right)\\
&= v(n_1,n_2)\ln(f(n_1,n_2))\\
\implies\\
\frac{\ln(g(n_1,n_2))}{v(n_1,n_2)} &= \ln\left(f(n_1,n_2)\right)\\
\implies\\
\mathbb{e}^{\frac{\ln(g(n_1,n_2))}{v(n_1,n_2)}} &= f(n_1,n_2)
\end{align}$$

But I don't have any idea how to estimate noise!

The question is: "Suggest the method for removing noise with transforming $g$ and filter this image in the transformed domain." I believe it is more about "finding a method" than "doing the job practically".

Your first step is quite sound: remove the power complexity. But you are left with a product of a noise and a ($\log$) image. Most well known tools operate in a linear way. So use your hypotheses a little further: since $f>1$, its logarithm is positive. hence $v\log f$ is positive too. So you can turn the noise/signal product into a sum, by applying a logarithm again:

$$ \log \log g = \log v + \log \log f\,.$$

Now  $ \log v$ can be considered as an additive noise to $\log \log f$. With this linearized model, you can try get insights on the nature of the noise, looking at flat regions of the image, checking histograms, etc. Once done,
you can  apply your best filter (adapted to the statistics of $ \log v$, yet unknown), and $\exp \cdot \exp$ the result. The function $v\to \exp  \exp v$ grows very fast (from rechneronline.de):

So it relatively flattens low values with respect to higher values, with can result in  contrast unbalance. This effect could be counterbalanced by intensity bounding, nonlinear mappings, etc., but more importantly  by proper image restoration.

For your information, using linear filters in a $\log$ domain is an instance  of homomorphic filtering, which your solution may belong too:

  Homomorphic filtering is a generalized technique for signal and image
  processing, involving a nonlinear mapping to a different domain in
  which linear filter techniques are applied, followed by mapping back
  to the original domain


Speech Enhancement Libraries Is there any speech enhancement algorithm library?

I need to process a sound to voice recognition, and I want to remove noise from audio source.

Nice a Brazilian Guy lol !

This is a little tricky to tell you without know what language are you working ...

Maybe a simple Filter can solve your problem or you can look at source code from Audacity are a cool code for noise removal.

Another choice is use Sox by command line to process your voice and use noiseprof to try remove noise.


sensing and representation matrices from exhaustive data If I know y (samples) and x (signal), is there a practical matlab algorithm to calculate A in y = Ax?

The solution is probably not unique.

$$    \left[ { \begin{array}{cc} a & b \\ c & d \end{array} } \right] \left[{\begin{array}{c} 1 \\ 0 \end{array}} \right] =    \left[{\begin{array}{c} 1 \\ 1 \end{array}} \right] $$

$b$ and $d$ can be any value provided $a = c = 1$


Relationship between Fs (the Nyquist frequency), and the frequency used with a sine wave Say you are given a vector t =[...] for the time to be used to generate a sine wave of a certain frequency say FHz. I have all the data for t so I can conclude what dt is (my dt = 1d-5). Now I am to do an FFT analysis and plot magnitude vs frequency so I expect to get my original amplitude and also get back the original phase
This what I do:

frequency = ...; %This in Hz

dt = t(2) - t(1); %which is equal to 1d-5

fs = 2*frequency; % should this be 1/dt? or something else?

Nfft = fs;

intensity = sin(2*pi*frequency*t);

FFT_intensity =fft(intensity,Nfft)/Nfft;

magIntensity = abs(FFT_intensity);

f=linspace(-fs/2, fs/2, Nfft);

plot(f,fftshift(magIntensity))

The results I get are sometimes correct and sometimes incorrect, it seems like when I chop the length of the t vector or change fs (say fs = 1/dt instead) I get the wrong amplitude ( and wrong phase, not shown here)... In general how to relate FHz, Fs and t?
I can actually include the numbers I am using if that makes any difference.
Thank you for any input.

The sampling frequency fs is 1/dt. Conceptually, if the time between samples is $T$, then you have $1/T$ samples per second, which is the sampling frequency.

Next, you need to make sure that frequency, the sine wave's frequency, is less than the Nyquist frequency, which is fs/2.

For interpreting the spectrum (FFT_intensity), please see the answer pointed to by Matt L.


Algorithm for object tracking In a video sequence, I would like to track an object. Right now, my detection algorithm can correctly locate the object in approximately 10% of the frames, and incorrectly find something else in 5% of the frames.

The incorrect detections are scattered around in the image, while my object movements are constrained by physics.

I've tried implementing a simple tracker based on distance from last detection (frames and pixels), but since my object moves in a 3d world, the velocity can vary a lot.

I think there is a rather simple approach for tracking this, but I don't have a lot of prior tracking knowledge. Where would I start?

(I appologize if this is not for DSP, please redirect me in that case)

assuming you have a relative high frame rate video (the object is moving relatively slow between two consecutive images) some standard tracking algorithms can be applied.

Mean shift [Wikepedia - Mean shift]
Use some template of the object you are tracking to create a likelihood function for the object location and track the object.

if you need more advanced method that can be more robust, e.g. object occlusions you can use more advanced methods:

kalman filter [Wikipedia - Kalman filter]
In short kalman filter uses a dynamic model that is composed of
prediction step and update step to track an object in an unsertain
enviroment.

More advanced option you can consider, if you have couple of objects you want to track at the same time is particle filter


Scaling Data before Fourier transform I am doing an experiment, where I am trying to get the same results from correlating two matrices in the fourier domain and then scaling the data by some constant.

Method_1: Do FFT, Correlate the data, Do IFFT, Scale data

Method_2: Scale data, Do FFT, Correlate the data, Do IFFT

In other words, I am expecting that if I scale the data before I do correlation I should get the same results if I do correlation first, and then scale the data afterwards(provided that the scaling factor is the same). What's different in my case is that I am doing the correlation in the fourier domain instead of the time domain.

Below illustrates an example of scaling data before or after fft, which will yield the same results

data_1=[1,2,3,6;4,5,9,10;7,8,9,4;10,14,13,12];
%scaling constant
h=magic(4);

%METHOD_1

    %go to fourier domain
    method1_data1_fourier=fft2(data_1,8,8);             %zero padding

    %do ifft
    result1=ifft2(method1_data1_fourier,8,8);

    %scale
    result1=result1(1:4,1:4)./h;

%METHOD_2

    method2_data1 = data_1;

    %scale
    method2_data1 = method2_data1./h;

    %go to fourier domain
    method2_data1_fourier=fft2(method2_data1,8,8);              %zero padding

    %do ifft
    result2=ifft2(method2_data1_fourier,8,8);
    %result2=result2(1:4,1:4);

Below is a matlab example of what I am trying to do in C++ with 32 bit floating point precsion and fft length of 47x47 (padded with zeros to 94x94 before transform)

%correlating data_1 against data_2
data_1=[1,2,3,6;4,5,9,10;7,8,9,4;10,14,13,12];
data_2=[7,2;9,8;3,1];
%scaling constant
h=magic(4);

%METHOD_1

    %go to fourier domain
    method1_data1_fourier=fft2(data_1,8,8);             %zero padding
    method1_data2_fourier=fft2(data_2,8,8);

    %do correlation
    result1 = method1_data1_fourier.*conj(method1_data2_fourier);

    %do ifft
    result1=ifft2(result1,8,8);

    %scale
    result1(1:4,1:4)=result1(1:4,1:4)./h;

%METHOD_2

    method2_data1 = data_1;
    method2_data2 = data_2;

    %scale
    method2_data1 = method2_data1./h;

    %go to fourier domain
    method2_data1_fourier=fft2(method2_data1,8,8);              %zero padding
    method2_data2_fourier=fft2(method2_data2,8,8);

    %do correlation
    result2 = method2_data1_fourier.*conj(method2_data2_fourier);

    %do ifft
    result2=ifft2(result2,8,8);

Unfortunately I am not getting the results I expected to get, as in result1 and result2 do not match at all. Can anyone tell me what is wrong with the above method I am using? I believe it has something to do with the result be being shifted after the correlation in the fourier domain because indice (3,3) of result1 and result2 are the same, and the first matlab example works fine.

I figured out the problem. But there is one constraint

(1) The scaling factor for the entire matrix has to be the same number to do the process efficiently

So the way correlation is done, is multiplication and accumulation of two matrices, take for example the image matrix below correlated with the signal, and the result is from the valid part of the correlation, as in no zero padding is done

image        signal    scaling constant = 3
1  2  3      9 4
4  5  6      3 4

result[1] = (1*9+2*4+4*3+5*4)/3 = 16.3333
result[2] = (2*9+3*4+5*3+6*4)/3 = 23

If I want to scale before the correlation, I would have to scale the whole correlation window, like so

image              signal    scaling constant = 3
1/3  2/3  3/3      9 4
4/3  5/3  6/3      3 4

result[1] = ( (1/3)*9+(2/3)*4+(4/3)*3+(5/3)*4 ) / 3 = 16.3333
result[2] = ( (2/3)*9+(3/3)*4+(5/3)*3+(6/3)*4 ) / 3 = 23

After scaling like above, I could alternatively go to fourier domain with the image and the signal and do correlation there, do ifft, and I still get the same results.

image              signal    scaling constant = 3
1/3  2/3  3/3      9 4
4/3  5/3  6/3      3 4

image_fourier=fft2(image,4,6);   zero padding
signal_fourier=fft2(signal,4,6);
result=ifft2(image_fourier.*conj(signal_fourier));

NOTE: above is matlab syntax

result[1] = 16.3333
result[2] = 23;

But this will not work if the scaling factor is different for result[1] and result[2].


Received Morse Code signal-to-noise ratio? I have some digital recordings of some received over-the-air CW broadcasts of known plain-text (W1AW code practice) at known WPM code speeds, but in moderate to heavy RF noise and fading.

Can I describe a S/N ratio for the Morse Code in these recordings?  If so, how would I calculate this S/N ratio?  What bandwidths should I use for the signal and for the noise?  Should I use recordings with or without the radio receivers narrow-band CW filter enabled?

If possible, I would like to use these S/N numbers to try and plot the decoding error rates versus the S/N ratios for various "demodulators" (including both software decoders, as well as humans at different experience levels).

Personally, I'd use whatever filtering scheme the actual receiver uses. If the narrow-band CW filter is exactly this, or even close, I'd use that. Then try to automatically (or manually) determine the regions where the signal is present. Measure the average power in the "off" regions, and that will be N. Measure the average power in the "on" regions, and that will be S+N. S = S+N - N, and then take 10*log10( S / N ). This will give you S/N ratio for the bandwidth of the narrow-band filter.


What are some really accurate ways to get the value of a peak (local maximum) given some points around it? (To be used for autocorrelation peaks.) I have looked everywhere on the internet for this and, surprisingly, haven't found much useful information.
Given 3 or more points closest to a peak (local maximum) what are some of the most accurate techniques I can use to obtain the value of the peak. In fact, I would like the x axis coordinate of that peak to determine the exact peak locations of the autocorrelation function of a discrete signal.
The best I have been able to do is quadratic interpolation which fits a parabola between three points and finds the maximum as follows:
\begin{equation}
x = \frac{1}{2} \frac{\alpha - \gamma}{\alpha - 2\beta + \gamma}
\end{equation}
Equation obtained here.
I tried cubic interpolation but either I'm doing it wrong or it's just not more accurate than the quadratic interpolation.
Anyway, can I not get better accuracy than quadratic interpolation? I would appreciate even just the names of algorithms or links to resources I can look at. (PS: I have way more than just three points close to the true peak).
Edit:
I'm not sure about the exact model of the peak, but it's one obtained when you do the autocorrelation of a sinusoid. To be more specific, I would like to use autocorrelation to determine, as accurately as possible, the exact fundamental frequency of an electrical grid.

Assuming the autocorrelation of a sufficiently long piece of sinusoidal signal has a good SNR (which it probably will), yes, that quadratic interpolation is quite good, as the sinc (which is actually the interpolating function) is locally quite close to the shape of 1-x². I've done similar things to estimate locations from signal crosscorrelations.
It gets problematic when the noise in your sines is high, or when you try to achieve a resolution much smaller than the sample distance (as the sinc gets quite "flat" when you zoom into it around 0 – and thus, the smallest bit of numerical error or noise leads to a large mistake).
Often, people are pretty happy after sinc-interpolating their signals, and it's not a very computationally hard thing to do:

Take your signal, as is
DFT it (using the FFT of signal length)
add plenty of zeros to the DFT's end (e.g. to pad it to the next larger power of 2, which increases the length by some factor $\alpha$ of your choie)
take the inverse DFT (using the IFFT)
Do your crosscorrelation, without interpolating

Tadah, $\alpha$-increased resolution.
This isn't great in terms of estimation variance (Fisher information didn't change, since the signal energy stays the same).
Since you know you're dealing with a dominant sine, a parametric spectral estimator would perform well. For a line spectrum ESPRIT is an appropriate choice, but autoregressive models lead to other good approaches like Yule-Walker, etc. If you're working in  a framework that can offer ESPRIT, I'd honestly try it.


Reconstruct FIR frequency response from the coefficients There is a lot of tools for designing FIR filters. Using one such tool I have designed a filter that satisfies my needs. It's coefficients are as follows:

 0.000106519966569472
 0.000164942422133671
 0.000154198611135164
-0.000051344150552581
-0.000527514219352601
-0.001228193882673074
 ...
 0.176493072494376790
 ...

In my system, 0.0001 / 0.1764 is about the same as ADC noise (well, not quite, but close). I'd like to know if a shorter filter, with some taps stripped off (zeroed) on either side still satisfies my requirements. This could be done by taking the modified array of the coefficients and reconstructing the frequency response of the modified filter.

Any tools for it (the cheaper the better)?

Edit

The tool works great. I post the result below. Good news is that only higher frequency response is visibly altered. Bad news is that by quite a lot. I am more interested in pass-band and slope, but I will look at more filter design tools to see if any of them has this sort of tradeoff built in. Thanks jojek.

You can try to window the filter with some window function, i.e. Gaussian, to get rid of small coefficients (taper them). Although it won't really work very well and you might really want to think about redesigning your filter - it's the best solution in my opinion.

Here is an example for a bandpass filter, showing that it is not really working well and is affecting the frequency response. Everything is done in Python so presumably  it's meeting your requirements of not spending a single penny.

Code is as follows:

import numpy as np
import matplotlib.pyplot as plt

# Some sampling frequency
fs = 48000.0
# Size of FFT analysis
N = 1024

def fir_freqz(b):
    # Get the frequency response
    X = np.fft.fft(b, N)
    # Take the magnitude
    Xm = np.abs(X)
    # Convert the magnitude to decibel scale
    Xdb = 20*np.log10(Xm/Xm.max())
    # Frequency vector
    f = np.arange(N)*fs/N

    return Xdb, f

if __name__ == "__main__":

    # FIR filter coefficients
    b = np.array([-0.0159147603287189,  0.000745724267485,
                   0.0404251831063147,  0.019042013872129,
                  -0.0644535904607569, -0.054523709591490,
                   0.0787623967281351,  0.105430811100048,
                  -0.0645610841865355, -0.148306808873938,
                   0.0257418551415616,  0.166568575042643,
                   0.0257418551415616, -0.148306808873938,
                  -0.0645610841865355,  0.105430811100048,
                   0.0787623967281351, -0.054523709591490,
                  -0.0644535904607569,  0.019042013872129,
                   0.0404251831063147,  0.000745724267485,
                  -0.0159147603287189])

    # Window to be used
    win = np.kaiser(len(b), 15)
    # Windowed filter coefficients
    b_win = win*b

    # Get frequency response of filter
    Xdb, f = fir_freqz(b)
    # ... and it mirrored version
    Xdb_win, f = fir_freqz(b_win)

    # Plot the impulse response
    plt.subplot(211)
    plt.stem(b, linefmt='b-', markerfmt='bo', basefmt='k-', label='Orig. coeff.')
    plt.grid(True)
    plt.hold(True)
    plt.stem(b_win, linefmt='r-', markerfmt='ro', basefmt='k-', label='Windowed coeff.')
    plt.plot(win*b.max(), '--g', label='Window function')

    plt.title('Impulse reponse')
    plt.xlabel('Sample')
    plt.ylabel('Amplitude')
    plt.legend()

    # Plot the frequency response
    plt.subplot(212)
    plt.plot(f, Xdb, 'b', label='Orig. coeff.')
    plt.grid(True)
    plt.hold(True)
    plt.plot(f, Xdb_win, 'r', label='Windowed coeff.')

    plt.title('Frequency reponse')
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dB]')
    plt.xlim((0, fs/2)) # Set the frequency limit - being lazy
    plt.legend()

    plt.show()


Compute SQNR (Signal to Quantization Noise Ratio) I'm studying the quantization of an audio signal and in particular the SQNR (Signal to Quantization Noise Ratio).

The book on which the study says that:

where:

N is the number of bits in the digital representation
V indicates that the signal voltage varies between -V and +V.

I don't understand how these calculations were carried out.
Is SQNR not simply the "digital" version of the RMS?
I know the RMS is calculated as power signal divided by power error, am I wrong?

The signal voltage is a sinusoid signal variant from $-V$ to $+V$, thus the RMS value

$$V_{rms} = \frac{V}{\sqrt{2}} = \frac{2^N \Delta V / 2}{\sqrt{2}}$$

because $2V = 2^N\Delta V$ where $\Delta V$ is quantization step.

Quantization noise is modeled as uniform random variable in $[-\Delta V /2, +\Delta V /2]$ thus its standard deviation is $V_n = \Delta V / \sqrt{12}$.

SQRN is defined as $\mathrm{SQNR}=\frac{V_{rms}^2}{V_n^2}$

$$\sqrt{\mathrm{SQNR}}=\frac{V_{rms}}{V_n} = 2^N \sqrt{\frac{3}{2}}$$

$$\mathrm{SQNR}_{dB}=20\log_{10}\left(\frac{V_{rms}}{V_n}\right) = 20\log_{10}\left(2^N \sqrt{\frac{3}{2}}\right) = 6.02N + 1.76$$

You must have mistyped your question.


Gardner Timing Recovery for Repeated Symbols I am implementing the Gardner Timing Error Detection algorithm for a QPSK signal. I understand the concepts behind it but I am curious as to what effect repeated signals can have on the phase lock.

If a symbol is repeated (maybe twice, maybe one hundred times) would the error not accumulate out of control? As I understand it the Gardner method looks at peaks and zero crossings to determine the timing error. If for, example, a symbol is repeated many times, would the algorithm not spiral out of control? If the symbol never changes, the algorithm would keep determining that there was symbol error and push the loop filter further and further from the true symbol timing; only fixing itself after a new lock in time when symbol transitions start occurring again.

Am I correct in thinking this will happen? If so, is it mitigated by the assumption that real world data will never (or very infrequently) repeat symbols more than a couple times in a row?

Thank you,

Your Gardner Loop would properly "Flywheel" in the absence of symbol transitions. So it would slowly drift off of synchronization, but would not spiral out of control as you describe. This is ideal for a control system with "missing updates"; to have the contributions for those updates with no additional information to contribute 0 to the accumulated error, resulting in a "flywheel" behavior. There are other detectors used to discriminate the control term in loops (such as a digital phase-frequency detector used in PLL implementations (such as this device: Maxim Phase Frequency Detector) that produce maximum error in the absence of transitions, and these systems would move off of lock at the maximum speed possible - this would be bad for timing recovery where we expect to have missing bits!

See my picture below demonstrating the operation of the Gardner Loop. In the absence of transitions, the Timing Error Detector would be zero on average as the Late and Early samples would be the same, and therefore $Late - Early = 0$.

Here is another post I made earlier with more information on the Gardner Loop:

Location of Matched Filter


Time Period of one Cycle is way too large I am calculating Time period of Sine/Cos Wave with different frequencies while keeping the Sampling Frequency at 1MHz and plotting it, but in some cases the value of time period is too big,

have a look at some of the examples.

n=[0:0.1:10] %because Time period is 10

F1=100*10^3; %100KHz  and t=1*10^-5
Fs=1*10^6;  %1MHz => Sampling Frequency
T=1/Fs;    %1*10^-6

%Time period = t/T => 10

xn=cos(2*pi*(F1)*n*T)   %T=1/Fs
stem(n,xn);

ylim([-1.5 1.5]);
xlim ([0 10]); %because Time period is 10
grid on;

Here Time period = 10 which means one oscillation would take 10seconds to complete.
and here is the output as i expected,

Similarly in second case when i changed the frequency to 10Hz while the Sampling frequency remains 1MHz the value of Time period is way too big,

n=[0:1000:100*10^3] %because Time period is 100*10^3

F1=10; %10Hz  and t=0.1
Fs=1*10^6  %1MHz
T=1/Fs;    %1*10^-6

%Time period = t/T => 100*10^3

xn=sin(2*pi*(F1)*n*T)

stem(n,xn)

ylim([-1.5 1.5])
xlim ([0 100*10^3])  %because Time period is 100*10^3
ylim ([-1 1])
grid on

and the output wave is,

is it the correct value or am i doing some sort mistakes calculating the time period ?

The equation xn=cos(2*pi*(F1)*n*T) expects n to be a sample index vector (n=0:N) where N is the total number of samples.  n*T or n/Fs produces the time vector.

Here is a way to plot a single cycle:

F1 = 100e3;
Fs = 1e6;
t = 0:(1/Fs):1/F1;  %Creates time vector for a single cycle
x = cos(2*pi*F1*t);
stem(t, x);

Note that 1/F1 gives the time period of a single cycle.


Algorithm for detecting whether a signal stored in rolling buffer is Increasing, Decreasing or staying Constant I have a data signal which is generated at constant time interval. It is a real world measurement and thus exhibits some noise in its value (it fluctuates approx + or - 20 units). The noise is fairly consistent.

I obtain the measurement at a regular time interval and stash it inside a rolling buffer of sample data (done via a computer). This rolling buffer is a list of values which are located at positions (0,1,2,3,4,5... -> n) in the list. Once the nth position in this list is reached, my buffer rolls over and overwrites the value at the 0th list position followed by the value at the 1st list position and so on until the nth value is overwritten (then rollover again). The most current data point value is always tracked (so I always know the position in the list of values at which I am writing my data to). This data buffering process repeats over and over. I use the rolling sample data to calculate a moving average and a std deviation.

Now for my problem. The signal has the potential to exist in 3 states, increase at a constant rate or to decrease at a constant rate at any time t1 or to remain constant. An increase or decrease can suddenly happen at any point in my list and will continue to some time t2 after t1.

I want to detect whether the signal is Increasing, Decreasing, or is remaining Constant (within the typical noise profile) with reasonable statistical certainty. What sort of algorithm would effectively notice this?

What sort of algorithm would effectively notice this?

What you are describing, is an online algorithm that uses a small set of "past" measurements to make an inference, in this case, about the slope of the signal.

Since you are keeping track of the current beginning and end of the rolling buffer, it is essentially like having access to a linear buffer of N past samples.

So, provided that the phenomenon you are trying to study (your signal) evolves slower than the length of the buffer you have chosen, you can apply a plain simple slope estimation on the values of the buffer and then use the slope estimate value (and its residuals) to decide if the signal is mostly heading downwards or upwards or remains constant. The most simple "classifier" here, would be a plain simple threshold one, possibly with the addition of hysteresis as it has been suggested already. More complex "classifiers" could be monitoring both the first and second derivative and returning additional information on whether the signal also changed direction within a buffer.

If you want to be getting one new estimate every time a new sample is received, you would have to look for iterative (or "recursive") methods, (for example, iterative least squares).

Hope this helps.


Reconstruct DWT for each cD1,cD2,cD3 and cA3 signals This question must be basic for this forum, but I'm only start working with DWT recently, and I was working with CWT before.
I think I can decompose my signal to L level and then get the scalogram with the following scripts.
 my question is how can I reconstruct the signal for each details{cD1,cD2,cD3} and approximation CA3?

Here is what I have tried so far!

t=0:0.1:300;
f=[0.4 0.6 1];
y=sin(2*pi*f(2)*t)+sin(2*pi*f(2)*t)+sin(2*pi*f(3)*t);
len=length(y);level=10;
[c,l]=wavedec(y,level,'coif5');
cfd=zeros(level,len);
for k=1:level
    d=detcoef(c,l,k);
    d=d(ones(1,2^k),:);
    cfd(k,:)=wkeep(d(:)',len);
end
image(flipud(wcodemat(cfd,255,'row')));

P.S: Actually I want to estimate the frequency components in a signal, so I thought I must decompose my signal to few levels then take FFT from the details to estimate the components that involve is the main signal. Is that right? or a better idea?
for the CWT case I was estimating the frequency components directly.

Alright, I'll try to take a shot at giving you an explanation of what's going on with the DWT.

So in the CWT, basically what you are doing is you are generating the wavelet coefficients by convolving the signal with each scale and shift of the mother wavelet function. However, the problem with doing this is that when you analyze the low frequency components of the signal, you are oversampling, or performing the wavelet analysis too often. Also, the relationship between time and frequency localization is much like the relationship between velocity and position localization in quantum mechanics, and when you perform the CWT this inequality is violated for some time/frequency scales. Essentially, with regard to the CWT, the transform is redundant at low frequency scales.

However the DWT doesn't have this problem because the DWT effectively downsamples the signal at each level. The DWT is basically just two filters, a highpass and a lowpass filter. The highpass filter produces the "detail" coefficients, and the lowpass filter produces the "average" coefficients. However, keep in mind that calculating the moving average of a signal is a form of lowpass filtering. Thus when we perform the lowpass filter and get the "average" DWT coefficients, those coefficients are basically just the low frequency components of the signal. And, since we know they are lower in frequency, we know that in the digital case we need less data points to represent a signal lower frequencies. The lowpass and highpass filtering is then repeated for each level of the DWT that can be computed. The level is dependent on the length of the signal. If a signal has $$k$$ data points, where $$k = 2^N$$
then the DWT can be computed to N levels.

I'll try to give a real life example. Lets say you have a signal of 1024 samples, sampled at a rate of $$4\pi \text{ Hz.}$$ Based on the Nyquist criterion, the highest frequency that can be resolved by the signal is $$2\pi \text{ Hz.}$$Computing the first level of the DWT will result in two vectors, each of length 512. One is from the highpass filter, the other from the lowpass (these vectors are sometimes called high/lowpass filter banks). The highpass coefficients are the cD1 coefficients, representing frequencies in the range of one to two pi, and the lowpass coefficients are the cA1 coefficients representing frequencies in the range of zero to one pi. The DWT is applied again to the lowpass coefficients resulting in two vectors each of length 256, and this time the lowpass filter bank represents frequencies in the range 0 to 1/2 pi. This process is repeated 10 times (for this case, 2^10 = 1024), as the last time will result in lowpass and highpass filter banks each of length 1.

Basically at higher scales of the DWT, we get better frequency resolution, but poorer time resolution. Lower scales have better time resolution but poorer frequency resolution. I'll also give you the link to a good tutorial on this stuff by Polikar : http://person.hst.aau.dk/enk/ST8/wavelet_tutotial.pdf

To reconstruct the signal, simply apply the inverse DWT for the number of levels which you want to reconstruct.

EDIT: Forgot to mention that the FFT is not needed, you are already filtering based on frequencies.


Energy of Sinc function How to find energy of sinc function $x(t) = \frac{\sin \pi t}{\pi t}$ with out the help of Fourier transform or Parseval's theorem?

The energy of $x(t)$ is given by

$$E_x=\int_{-\infty}^{\infty}x^2(t)dt=\int_{-\infty}^{\infty}\frac{\sin^2(\pi t)}{(\pi t)^2}dt\tag{1}$$

If we may assume that we know that $x(t)$ is the impulse response of an ideal low pass filter, the integral $(1)$ can be computed without using the Fourier transform and Parseval's theorem by noticing that it can be represented as a convolution evaluated at $t=0$:

$$E_x=(x\star x)(0)=\int_{-\infty}^{\infty}\frac{\sin(\pi\tau)}{\pi\tau}\frac{\sin(\pi(t-\tau))}{\pi(t-\tau)}d\tau{\huge|}_{t=0}\tag{2}$$

The convolution in $(2)$ can be seen as filtering an ideal low-pass signal with an ideal low-pass filter with the same cut-off frequency as the signal. Clearly, the filtering does not affect the signal, and the result of the convolution equals $x(t)$:

$$(x\star x)(t)=x(t)\tag{3}$$

And, consequently,

$$E_x=x(0)=1\tag{4}$$

(Note that the value of $x(t)$ at $t=0$ is defined as a limit.)


Computation of Only Even or Odd Frequency Bins of DFT I have an algorithm where I am computing the FFT of a large signal. However, I desire only the even or odd terms of the DFT of the signal, but not both. Currently, I discard these undesired terms. Is there a way by which I may appreciably reduce the runtime complexity when only every other DFT term is desired?

How might I begin to find such an optimization? If it is unclear, I will look for how the larger, encompassing algorithm could be improved, but I'm starting by examining the smallest, most expensive portion before thinking about more fundamental changes to the system.

Given a signal $x[n]$ of length N, assume $X[k]$ is its N-point DFT $$X[k] = \sum_{n=0}^{N-1}{x[n]e^{-j\frac{2\pi}{N}nk}}$$ where $k = 0, 1,..., N-1$  computed via N-point FFT. And if we wish to compute separateley those even and odd indexed samples of the DFT $X[k]$ , we can proceed with the following:

Consider those even indexed samples of $X[k]$ which are denoted as $X_e[k]$ where:
\begin{align}
X_e[k] = X[2k] &=  \sum_{n=0}^{N-1}{x[n]e^{-j\frac{2\pi}{N}n(2k)}} &\scriptstyle{\text{term 2k computes even samples of X[k]}}\\
X_e[k] &=  \sum_{n=0}^{N-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} &\scriptstyle{\text{factor 2 is moved into N/2 term}}\\
X_e[k] &=  \sum_{n=0}^{\frac{N}{2}-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} + \sum_{n=\frac{N}{2}}^{N-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} &\scriptstyle{\text{Divide the sum into 2 halves}}\\
X_e[k] &=  \sum_{n=0}^{\frac{N}{2}-1}{x[n]e^{-j\frac{2\pi}{N/2}nk}} + \sum_{n=0}^{\frac{N}{2}-1}{x[n+\frac{N}{2}]e^{-j\frac{2\pi}{N/2}nk}} &\scriptstyle{\text{adjust the 2nd sum range}}\\
X_e[k] &=  \sum_{n=0}^{\frac{N}{2}-1}{ (x[n]+x[n+\frac{N}{2}])e^{-j\frac{2\pi}{N/2}nk}}  &\scriptstyle{\text{merge the 2 sums}}\\
\end{align}

The final form suggest that the required even samples of N-point DFT of signal $x[n]$ can be computed from an N/2 point DFT of a new signal $x_{half}[n]= x[n]+x[n+ N/2]$ whose length is half that of $x[n]$

The signal $x_{half}[n]$ is simply computed by adding the first half of the signal $x[n]$ into its second half. (if the length N of signal x[n] is not even padd a zero to its tail to make it even)

In short the following matlab code gives you the desired even samples of X[k] without computing a full N point DFT (via FFT) of x[n]:

Xe = fft ( x(1:N/2) + x(N/2 + 1 : end), N/2);

Note that because of the addition of those half length signals efficiency will degrade from a simple N/2 point FFT.

The case for the odd indexed samples proceeds similary but results in a more complex form which I would like to summarize in the following Matlab script:

% Let x[n] be the signal of length N
xc = x .* exp(-j*2*pi*[0:N-1]/N) ; % Multiply x[n] with a complex phase term
Xo = fft (xc(1:N/2) + xc(N/2 + 1, N), N/2); % odd indexed samples of X[k]

Due to a complex multiplication before the FFT, in this case of computing odd indices performance is further reduced but still preferable over a direct N point FFT.

The following Matlab excerpt demonstrates both even and odd samples (with odds requiring only half the complex multiplication wrt the above line)

N = 1024;
x = randn(1,N);
X = fft(x,N);       % original DFT of X

Xe = fft( x(1:N/2) + x(N/2 + 1: N) , N/2);                          % Even   samples only from N/2 point FFT
Xo = fft( ( x(1:N/2) - x(N/2+1:N) ).*exp(-j*pi*[0:2:N-2]/N) , N/2); % Odd samples only from N/2 point FFT

figure,stem( abs ( Xe - X(1:2:N)) ); title('Xe - X(1:2:N)');
figure,stem( abs ( Xo - X(2:2:N)) ); title('Xo - X(2:2:N)');

Note the absolute value of the difference. It is neither zero (down to roundoff error), nor too big. But for the purpose of particular computation this "practical" error shall be considered carefully.


Local Image Features Extraction for Road Signs I'm working on a project to detect a road sign from the Image. Consider the sample image below:

I'm using a heuristic search to detect the road sign within an image. Basically, I want some weak features (colour, shape, zero symmetry within road signs) which can be calculated very fast and their weighted average can be used recognize the road sign.

Currently, I'm searching the entire image to find the road sign which could be of different scale. Since the domain of search is a lot, I tried to reduce it bu by using some heuristics:

Location of sign detect in previous frame (its mostly likely to be near)
Road signs are usually on right sign
Road sign contains enough white pixels

But eventually search needs to stop when it finds the region that contains a road sign. For this purpose, I'm relying on some fast and robust local features that predicts likelihood of containing road sign within given region of an image. Currently, I'm relying mainly on colour and some region props that Matlab provides.

I tried looking for lines and corners within edge map of the region even tried thresholding and looking for rectangle using morphology. But seems like they are too binary either its YES or its NO. In fact in most cases it fails, so my questions what are some nice local features that I could rely on to detect road sign?

I would really appreciate if somebody you guide in the right direction or provide some links or papers that could be helpful for this problem.

Do you know what signs you are looking for? If yes, maybe you could do template matching (e.g. a normalized cross-correlation, available in matlab). It won't work great when signals are getting closer, since the perspective projection will change their appearance, but it should work for mid-range detection.

You can limit your template-matching search to the right side of the road.

There are also template matching algorithms for affine transformations that might be a next step.

Christos


Image getting completely white on downsampling I performed an image downsampling with averaging operation(The image is divided into 2 by 2 blocks,then it is downsampled by 4 times.New pixel values od downsampled image is replaced with the mean of each block) as follows.

img = imread('cameraman.gif');
avgfilter = fspecial('average', [2 2]);
avgimg = filter2(avgfilter, img);
down1 = avgimg(1:2:end,1:2:end);
imshow(down1);

I used a graysacle image as input.

but I'm getting a complete white image as output on downsampling.
I tried with jpg image also. Why this happened?
Please help.

When you open an image using imread, you get the image as a uint8 matrix with grayscale (or color) values of 0-255. However, filtering is no possible for this type. Therefore, filter2 converts your image to double.
The function imshow accepts both uint8 and double type. When the input is uint8 it also assumes color values of 0-255. But, when the input is double, it assumes the values are normalised to the range of 0-1. Therefore you should tell imshow that although your input is of type double it has color values of 0-255 by using

imshow(down1,[0 255]);


System characterization given pole-zero mapping The figure shows the pole–zero plots for three different causal linear time-invariant systems with real impulse responses.
Indicate which of the following properties apply to each of the systems pictured: stable, IIR, FIR, minimum phase, all pass,linear phase, positive group delay at all $w$:

I am having problems specifically finding out if the system is linear phase and has a positive group delay, and all pass or not.

In the first plot: the system is stable since we have a region of convergence containing the unit circle, and FIR because it is stable, and minimum phase because we have a pair of complex conjugate pairs.

And so on with the other figures...But I still cannot quite figure out how to determine the other properties given the pole-zero mapping.

Ok, so I've gotten to this point:
Plot 1 is not stable because the ROC doesn't contain the unit circle and since all zeros are in the unit circle, we have a minimum phase system. Additionally, we have a IIR system because of the 3 poles so we have feedback. Its also not an allpass filter because poles and zeros are not reciprocals of each other. I hope those are right and I understand the concept so that I can continue with the other plots. But, I'm really having difficulty trying to find out if the system is a linear phase system or not and if it has a positive group delay. I know that linear phase translates in delay in time for all frequencies by the same amount and that group delay $\tau_g = -d\theta(w)/dw$.

Thanks in advance guys.

I will not do the homework for you, but I'll make a few comments on your observations to clear up some of your misconceptions. Concerning the first plot (top left I suppose) you wrote:

"the system is stable since we have a region of convergence containing the unit circle"

This would be correct if the system were non-causal (i.e., with a two-sided impulse response), but since - according to the first sentence in the question - the system is causal, it can't be stable if it has poles outside the unit circle.

" ... and FIR because it is stable"

This is wrong. What you're saying is that any stable system is FIR, which is definitely not the case. Of course there are stable IIR filters. Here you'll have to do some more studying to learn what stability means and to understand that it is not necessary (only sufficient) for a stable system to be FIR.

"... and minimum phase because we have a pair of complex conjugate pairs"

That's also wrong. The fact that a system has complex conjugate pairs of poles and zeros has nothing to do with the minimum-phase property. You need to read up on minimum-phase systems. They can be inverted by a causal and stable system. What does this tell you about its zeros?
Any system with a real-valued impulse response has either real-valued or complex conjugate pairs of poles and zeros. That's because $H(z)$ is a real-valued function of a complex argument.


Why 8-PSK modulation has higher Eb/N0 than QPSK or BPSK? I run a transmitter/receiver simulation using QPSK, BPSK, and 8-PSK modulations. As far as I have noticed, the best SNR per bit results are obtained when simulating with 8-PSK modulation (and also at longer distances). For example - at a distance of 10 m (transmitter to receiver) Eb/N0 for QPSK was 22.9 and for 8-PSK it was 26. Does it depend on transmitting more bits symbols per cycle (8-PSK transmits 3 while QPSK/BPSK only 2)?

example - at a distance of 10 m (transmitter to receiver) SNR for QPSK was 22.9 and for 8-PSK it was 26. Does it depend on transmitting more bits symbols per cycle (8-PSK transmits 3 while QPSK/BPSK only 2)?

This makes no sense, physically: The attenuation through the channel is the same for all these signals, since you'd put them all through the same pulse shape to have a fair comparison.
The fact that they don't transport the same number of bits per symbol doesn't matter to physics.
Maybe you're confusing $E_b/N_0$ with SNR?

Ah yes, you were confusing $E_b/N_0$ with SNR. Now, you didn't compare fairly. In BPSK, a bit gets the energy of a full symbol. That's the transmit power times the symbol period.
Since you need to compare fairly, you need to keep both the transmit power constant (obvious) and either the bit rate constant, or the symbol rate constant:

constant power, constant bit rate: Now your 8-PSK transmit signal has only 1/3 of the symbol rate of your BPSK signal. That means narrower bandwidth, 1/3 of BSPK. Which means 1/3 of the noise. You cheated here by reducing the bandwidth.
constant power, constant bandwidth: your system now divides the transmit power on three times as many bits as BPSK, so $E_B/N_0$ actually is reduced, since within the same bandwidth, the same original noise occurs.

Note that in (almost) any sensible comparison of transceiver systems, you keep the $E_b/N_0$ constant and vary the bandwidth or necessary transmit power – your receiver cares less about the SNR than it cares about how likely it is that a bit gets mis-decided, which is a function of the energy per bit. This becomes slightly more obvious when one considers what interleavers and channel decoders do: they distribute noise energy that were concentrated on individual bits or symbol-bit-groups over a larger duration, and then correct these less-than-perfect bits using the knowledge from all the other bits.


How to Extract Nonorthogonal PCA Principal Components My problem is essentially a 'blind source separation' problem. I have 3 non-orthogonal sources (or basis functions) and $N$ random linear combinations (mixes) of said sources. My problem is to obtain the sources from the mixes.

Figure A shows the sources, B shows the mixes.

Approaches taken:

$\textbf{PCA}$ - I tried PCA on the mixes (fig.C), but the issue is that PCA will only give orthogonal bases, while my sources/bases are non-orthogonal. This issue with PCA is shown in figure D, where the data is clearly described by 2 non-orthogonal basis, but PCA (the solid lines) cant reconstruct them!
$\textbf{Factor Rotation}$ - I tried applying some solutions form Factor Analysis (not my forte). The promax rotation (matlab: nw = rotatefactors(cov,'method','promax') ) is shown in fig. D with the dashed lines. As far as I can tell, factor rotation works with the principal components, not the original matrix and thus I have no idea how it could reconstruct the right basis. I think this only works with factor matrices, not with generic ones...
$\textbf{ICA}$ - I tried overdetermined ICA with the fastICA algorithm (also not my forte, but I think I understand it better). I was hoping this would work since independent components (ICs) are non orthogonal. The solution is shown in fig.E. While the ICs are indeed nonorthogonal, they are NOT my original sources :-(

Any other potential tips or leads or solutions would be greatly appreciated.

More mathematical Detail

The 3 source signals are themselves linear combos of the 5 laguerre basis functions. Thus, the problem can be formulated as such:

I take l=5 laguerre basis functions of length m=30 to get B, an [l,m] matrix
I randomly mix then is my x=3 sources, S: S=B*M1, where M1 is a random [l,x]=[5,3] mixing matrix. S is displayed in fig.A
I mix S n=100 times to get my mixed signals, K (fig.B), where K=S*M2, a [x,n]=[3,100] randomly generated mixing matrix.
I have the original basis, L, & the mixed matrix, K. K=S*M2=B*M1*M2. My goal is to recover S.

NOTE: I understand that this may seem impossible since there are many possible S's to choose from (ie the factorizatio S*M2=B*M1*M2 is not unique), however, the true S, will dramatically decrease model complexity by having only x=3 rather than l=5 basis functions which can exactly reproduce all the kernels.

Thanks!

Let's say x,y and z are your basis vectors and you mix them to get what you are calling your mixed sources. If your basis vectors were (x-y), y, and (z+y) then you could mix them in a different way and end up with the same mixed sources.

Basically, there are many ways to come up with a set of basis vectors for your signals. You need to say what kind of basis you want. It is not enough to say that the basis vectors are non-orthogonal.


Least Square Error Estimation: Conditions for $ (A^TA)^{ -1} = A^{-1}(A^T)^{-1} $? I am new to linear algebra and have this simple question...
in least sqaure estimation...the best estimation of the equation $Ax = b$ is $x_{Estimated} = A (A^t A )^{-1} A^t b$...the projection of $b$ on the column space of $A$ is $p = A (A^t A)^{-1} At b$ which is written as $p = P b...P$ is the projection matrix = $A (A^t A)^{-1} A^t$
now using inverse rule $(AB)^{-1} = B^{-1} A^{-1}$  .......$P$ reduces to $I$ ....so is the projection matrix $P$ just the identity matrix $I$ ..?? am i missing out some basic or is my concept not clear? can we write $(A^t A)^{-1}  = A^{-1} (A^t)^{-1}$ always??

First of all, the minimum norm least square solution is $A^+b$, where $A^+$ is the pseudoinverse. Only when the left inverse $A_L^{-1}$ (or right inverse $A_R^{-1}$) exists, you have $A^+=A_L^{-1}=(A^TA)^{-1}A^T$ (or $A^+=A_R^{-1}=A^T(AA^T)^{-1}$).

Likewise, the projection matrix onto the column space is $P=AA^+$, or $P=AA_R^{-1}=AA^T(AA^T)^{-1}=I$ if $A_R^{-1}$ exists. In other words, if A has full row rank, then the projector $P$ is indeed an identity matrix. If only left inverse exists, then $P=AA_L^{-1}=A(A^TA)^{-1}A^T$ is not an identity matrix.

Finally, $(AB)^{-1}=A^{-1}B^{-1}$ is valid only when $A$ and $B$ are square matrices.


Implementation of a filter in a software application I have to embed a bandpass filter into a software application.

All the signal processing will be done in 32 or 64-bit double precision.

The filter coefficients for an IIR Butterworth filter are computed offline in Matlab.

What is the best way to implement the filter in the application?

Right now my plan is to implement the filter as cascaded second order sections in transposed direct form II, but I have no practical experience (and literature on the practical part of implementation is also kind of sparse), so I'm just kind of guessing.

For example: Second order sections are used to avoid numerical problems, but do I really need to take care of that when the processing is done in floating point precision, also which form should i chose for the realization of the filter?

I would be thankful for any tips or literature hints.

If your filter is not of insane order, than double precision (64bit on most machines) floating point will totally suffice, and you really shouldn't worry too much. I agree with Arnfinn here.

I don't 100% agree on the optimization recommendation – in fact, if you want my two cents on general purpose CPU-based software implementations of filters,
then simply don't use IIRs; they might be shorter than FIRs with the same flatness/transition band steepness/suppression, but the fact that they are recursive makes them a whole lot harder to execute fast, not to mention that shortcuts like polyphased decomposition don't generally work, and thus, you can't save as much CPU when e.g. building a decimating anti-aliasing filter.

My honest advise is:

Don't worry too much about CPU performance. If you're not trying to do this in JavaScript, and you got but a couple of kilosamples to second to worry about, a simple filter implementation will do.
If you're worried about numerical accuracy, don't worry more than at the level that's sensible: unless your filter is multiple thousand taps long or your taps are extremely bad-behaved, a FIR won't really hurt much. If you're worried about stability, simply don't use an IIR.
Don't reinvent the wheel if you don't have to. There's excellent optimized DSP libraries out there, and as long as the programming language you use supports native libraries with a C calling convention ABI (it probably does), you can use most of them. Just let whatever tool you chose provide the taps for your IIR or FIR, and that library will do its job. I personally like libVOLK a lot, but it concentrates on FIRs, for the reasons given above.

Generally: You're asking about numerical accuracy without indicating which boundaries you need to work in and what filter orders you're talking about. This might indicate that you're not really having any problem at all, because otherwise, you might have noticed already that these are the critical aspects here :)


Is the basis of the sparse signal assumed known in compressed sensing? I'm new to compressed sensing, and am a little confused about the assumption of the basis matrix $\Psi$.  Is $\Psi$ assumed known in compressed sensing?

Specifically, suppose that a signal $x$ is sparse in some basis, say $\Psi$, i.e. $x=\Psi\alpha$, where $\alpha$ is $k$-sparse, i.e. $\|\alpha\|_0\le k$.  My understanding is that in compressed sensing, we store $y=Ax$, where $A$ is the $m\times n$ sensing matrix with RIP.  Later on, we can recover $x$ theoretically by

$$\hat\alpha=\underset{\alpha\::\:y=A\Psi\alpha}{\arg\min}\|\alpha\|_0,$$

(or use $l$-1 norm) and then $\hat x=\Psi\hat\alpha.$

In doing so, we need to know $\Psi$ to recover $x$, don't we?  But if we know $\Psi$, why don't we just measure $z=\Psi^{-1}x=\Psi^{-1}\Psi\alpha=\alpha$, and store the value and index of the non-zero components of $z$ ($=\alpha$)?  Can the dimensions ($m$) of $y$ be less than $2k$?  Or is it that $z$ won't be perfectly equal to $\alpha$ due to noise, and we don't want to do thresholding on the fly?  Or some other reason?

I'm confused about the rationale of compressed sensing, and would appreciate any pointers, comments, and clarifications.  Thanks a lot!

Compressive Sensing is an approach to reconstruct sparse signals from incomplete set of measurements.

  In doing so, we need to know $Ψ$ to recover $x$, don't we?

yes, we do.

  But if we know $Ψ$, why don't we just measure
  $z=Ψ^{−1}x=Ψ^{−1}Ψα=a$, and store the value and index of the
  non-zero components of $z (=α)$?

Compressive Sensing  as the name indicates is a method to sense the incoming signal not to compress it. Which means the signal is not present yet and we are going to sample it, this sampling is efficient and requires minimum number of samples for signal reconstruction. There is a misunderstanding here, we might know $Ψ$ in advance but we usually do not have the signal itself $x$, if we do, clearly compression techniques are preferred. The other reason is that sampling domain could be different from the domain in which the signal has a sparse representation. For example, we might sample a fourier sparse signal in time (e.g. radiowaves) or sample a fourier sparse image in space (e.g. Cameras, MRI,..). And then reconstruct it it in fourier domain.

  Can the dimensions ($m$) of $y$ be less than $2k$?

Observing this bands guarantees perfect reconstruction (with overwhelming high probability) for any $k$-sparse signal, but reconstruction is still probable for $m < 2k $; So, as the number of measurements ($m$) is reduced probability of reconstruction reduces as well. But if you mean maintaining the perfect reconstruction probability while minimizing $m$, I think it would depend on your matrix design.


Fourier coefficients of |A cos (x)| I have to find Fourier coefficients of $x(t) = |A \cos (2 \pi f t )|$. The problem also give me $y(t) = |x(t)|$. To find Fourier coefficients I wrote

$$ \frac{1}{T_0} \int_{0}^{T_0} |A \cos (2 \pi f t )| e^{2 \pi k f t}$$

Solving this I found that $k$ should be equal to  $1$ or $-1$. First I assume that $=1$ and I found the result $A/2 $ but on the book should be $2A/3\pi$. $$ T=\frac{T_0}{2} $$ Can someone help me?

Sketching the signal always helps in such cases. Let's consider a cosine of period $T_0$, that is, of fundamental frequency equal to $f_0=1/T_0$. That is the signal on the left.

As you can see on the right, $|A\cos(2\pi f_0 t)|$ does not share the same period with $A\cos(2\pi f_0 t)$. The period of the signal you are looking for is $T_0^{\prime}=T_0/2$ and its fundamental frequency is $f_0^{\prime} = 2f_0 = \frac{2}{T_0}$.

Assuming that you're looking for the exponential Fourier Series, the coefficients should be $$X_k = \frac{1}{T_0^{\prime}}\int_{T_0^{\prime}}x(t)e^{-j2\pi kf_0^{\prime}t}dt$$
and for $k=0$, which is usually calculated seperately, it is $$X_0 = \frac{1}{T_0^{\prime}}\int_{T_0^{\prime}}x(t)dt$$

Notice that you can integrate over any time interval of duration equal to $T_0^{\prime}=T_0/2$. I will pick $[-T_0/4, T_0/4]$ since the cosine is always positive in this interval and thus removing the $|\cdot|$ is trivial.

I will help with $X_0$ and you can do the $X_k$s. It is:

$$\begin{align}
X_0 &= \frac{1}{T_0^{\prime}}\int_{T_0^{\prime}}x(t)dt = \frac{2}{T_0}\int_{-T_0/4}^{T_0/4}A\cos(2\pi f_0t) dt \\
&= \frac{2}{T_0}\frac{A}{2\pi f_0}\sin(2\pi f_0 t)\Big]_{-T_0/4}^{T_0/4} = \frac{2A}{2\pi f_0T_0}\Big(\sin\Big(\pi f_0\frac{T_0}{2}\Big) -\sin\Big(-\pi f_0\frac{T_0}{2}\Big)\Big)
\end{align}$$

We know that $T_0 = 1/f_0 \Longrightarrow f_0T_0=1$ and that $\sin(-x)=-\sin(x)$. Simplifying:

$$\begin{align}
X_0 &= \frac{A}{\pi}\Big(\sin\Big(\frac{\pi}{2}\Big) + \sin\Big(\frac{\pi}{2}\Big)\Big) = \frac{2A}{\pi}
\end{align}$$

Now for the $X_k$s, I've already mentioned that you have to use
$$X_k = \frac{1}{T_0^{\prime}}\int_{T_0^{\prime}}x(t)e^{-j2\pi kf_0^{\prime}t}dt$$
You can use the same interval for your integration, keep in mind that $f_0T_0=1$ to simplify your expressions, and treat $k$ as a constant until the end. Judging from your question, you should get something that includes $\frac{2A}{(2k+1)\pi}$.


variable length no-refrence filter algorithm I'm searching for a an algorithm or an idea of implementation. I need to filter a signal in real time, and want to change the length of it (number of taps) in order to change the delay.

In low signal environment I want a high order filter and care less
about delay
In high signal environment I want a low order filter in order to
minimize the delay.

I've tried to search "variable length filter..." on google but all of the results were versions of the LMS algorithm, but I don't have a "clean noise" reference signal which I can use.

Have anyone heard of something like this?
Thanks :)

One possible implementation would be to run two (or more) different order filters in parallel, and use a mixer to blend them by ramping the gains up and down depending on your selection criteria.


Difference between a wavelet transform and a wavelet decomposition I'm confused about the difference between a wavelet transform and a wavelet decomposition is. For example

load woman
[cA1,cH1,cV1,cD1] = dwt2(X,'db1');
[c,s] = wavedec2(X,2,'db1');

What's the difference between these two matlab commands, and when would you want to do one over the other?

I don't think there is any difference.  The documentation for dwt2 says

  Single-level discrete 2-D wavelet transform

  The dwt2 command performs a single-level two-dimensional wavelet decomposition...

While the documentation for wavedec2 says

  Multilevel 2-D wavelet decomposition

The difference is that dwt2 is single-level (produces a single A, H, V, D output):

and wavedec2 is multilevel (produces array C output, which contains multiple A, H, V, D inside it):


How do I obtain the frequencies of each value in a 2D FFT and what do they mean? My goal is to analyze the frequencies in an image representing water simulation data before and after, for example, a Gaussian filter. The direction of these frequencies is not important to me. Ideally, I think I want to plot these frequencies in a 2D graph, where x represents the frequency, ranging from 0 to 0.5, and y represents the amplitude.

I understand how to obtain the 1D frequencies, as is explained here:
https://stackoverflow.com/questions/4364823/how-do-i-obtain-the-frequencies-of-each-value-in-a-fft. From this, I also understand how to obtain the frequencies in x- and y- direction.

What I do not fully understand is what the direction-independent frequencies of my output values are. Intuitively, I would expect to calculate them like so:

// from output to frequencies to periods
fx = kx / N; fy = ky / M
px = 1/fx; py = 1/fy

// direction independent frequency of output value
period = sqrt(px*px + py*py)
frequency = 1/p

However, when I read about this particular subject, I find something different:

From what I understand, u and v are the x- and y- frequencies respectively, in cycles per unit distance. This would mean that, if I have a signal with x- and y- frequencies of both 0.5, the resulting frequency is 0.707... ? This is over the Nyquist limit, even more so diagonally! I have a feeling I am missing something, but I can't find what.

Additionally, what do the frequencies with an x- or y- frequency of 0 or 0.5 mean? In 1D, they would be the DC and Nyquist component. How does this translate to 2D, where the other may have a perfectly fine value? More practically, which frequencies can I use for my graph?

It's not over your Nyquist limit, because you're doing sampling at your nyquist rate in two independent dimensions and look at the combined result.

Your square-root formula is easily explained: It's pythagoras, matching your combined frequencies in x and y to a single frequency.


Exact indices of $\tt fft$ Let's consider the cyclic frequency transform and let's try to make it so that the transform of the function is itself.

t=linspace(-32,32,4096);
dt=t(2)-t(1);
fq=linspace( -1/(2*dt) , 1/(2*dt) , length(t) + 1 );
fq=fq(1:  length(t)  );
x=exp(-pi*t.*t);
X=fftshift(fft(x)).*exp(-1i*2*pi*fq*t(1))*dt;
plot(t,real(X));hold on;plot(t,imag(X));hold off;

Now this actually works, and believe me it was not easy for me to get these right. I basically resorted to trial and error, tried adding $\pm 1$ everywhere until I got something that looked reasonable at all. I would like to learn the method behind this madness.

Is is a good idea to use something like t=linspace(-32,32,4096);
Considering how the frequency is asymmetrical, perhaps time should be as well?
What frequencies does the result of X=fftshift(fft(x)) actually correspond to?
Any other suggestions, etc.?

The FFT calculates the spectrum at N points in the range from $[0, 2\pi)$ where due to the cyclic property the relative frequency 0 must be equal to the frequency $2\pi$. This is best shown via an example. Note that $2\pi$ corresponds to your sampling frequency.

Example:
For a length 3 impulse response, say $h=[0.5, 1, 0.5]$, the FFT calculates the spectrum at the relative frequencies $[0, \frac{1}{3}2\pi, \frac{2}{3}2\pi]$, instead of the frequencies $[0, \frac{1}{2}2\pi, 2\pi]$. One way to calculate these frequencies in Matlab for an $N$ points FFT is as follows:

N = 3;
h = [0.5 1 0.5];
fq = linspace(0, 2*pi, N+1);
fq = fq(1:end-1);
H = fft(h);
figure;
plot(fq, abs(H));

The function fftshift takes the right part of the spectrum and places it on the left such that you obtain the frequencies $[\frac{2}{3}2\pi, 0, \frac{1}{3}2\pi]$.

When you are not using the fftshift function and just plot the range $[0, 2\pi]$, then you can use positive time indices for t.


Why should the number of subcarriers of OFDM system be a power of 2? I notice that almost all papers use the number of subcarriers of 64, 128,256.. so on, which is a number of $2^k$ and $k$ is an integer.

Why is the number of subcarriers always like that? Can we use any number or it must be only using one of those numbers?

Implementing an IFFT and FFT engine which is a power of $2$ is easier in hardware and hence if an OFDM system is talked about, it is talked about in $2^k$ length FFT. All practical Communication system based on OFDM use $2^k$ length FFT-IFFT engines.

However, for OFDM in principle, it is not required to have FFT-IFFT in power of $2$.

There are different factors which may affect the choice of number of sub-carriers based on multi-path channel response in the frequency range of operation and maximum doppler spread the channel expect to suffer.