Applying Nyquist theorem to digital vs. analog audio quality I had a recent conversation with a friend about analog vs. digital recordings.  He shared his opinion that music recorded by analog techniques is of superior sound quality to music recorded by digital techniques.  His reasoning was that whenever music is turned into digital, not all of the information is captured, and thus, some of the music is lost.

I was not too sure about his statements as I recently took a DSP course and learned a bit about the Nyquist-Shannon theorem.  As human hearing is typically in the range of 20 Hz to 20 kHz, I was under the impression that the common sampling rate of 44.1 kHz should be enough to capture all information that is needed to reproduce the sounds that we can hear.  Therefore, I thought that differences between music reproduced from digital and analog should unnoticeable to human ears.

Who is correct?  If anyone, could shed some light on this and clear up any misunderstandings, I would greatly appreciate it.  Thank you very much.

By the way, I am new to Stack Exchange, and I have already tried to find similar questions but to no avail.  Hopefully, this is not a repost.

The Nyquist-Shannon reconstruction theorem applies to perfectly bandlimited signals and perfectly flat reconstruction filters.  Unfortunately, neither exists in the real world.  Nor do digital system that use an infinite number of bits per sample to avoid quantization noise.

The headroom allowed by using higher sample rates than 44.1k may allow a closer approximation (for the bandlimit and reconstruction filter response) to the required perfection.  So does using more bits than 16 to reduce quantization noise.

But analog recording is also not perfect, as noiseless analog recording media and head amplifiers do not exist (because of the existence of atoms and non-zero K temperatures).

So there can be a trade-off.  Choose your preferred sounding reproduction artifacts(s).


Normalized cross-correlation in detail I'm trying to implement a normalized cross-correlation algorithm but I don't get what in fact is this measure. What confuses is the wikipedia definition:

$\frac{1}{n} \sum \frac{(f(x,y)- \overline{f})(t(x,y)- \overline{t}) }{\sigma_{f}\sigma_{t}}$

Which result is an scalar (AFAIK)

But then adds other way to measure it:

$\left \langle \frac{F}{\left \| F \right \|},\frac{T}{\left \| T \right \|} \right \rangle $

Where $F$ and $T$ are normalized vectors and $\left \langle . , . \right \rangle$ is the inner product. But the output will be a vector, isn't? Isn't supposed to give me a scalar as well? Am I getting something wrong?

The idea is to implement this formula and use it with matrices with same dimensions.

Thank you

An individual inner product does produce a scalar, but often when a cross correlation is calculated multiple individual cross correlations (i.e. dot products) are calculated at different time offsets.  These individual scalar results form a vector that is indexed by the relative time offset.


Chosing polyphase interpolation coefficients I'm writing code for a 3x polyphase interpolator using a total of 9 coefficients. This is organized as 3 parallel FIR branches with 3-taps each that are summed at the 3x sampling rate.

The problem I have is the impulse response looks correct (I get my 9 filter coefficients as output) but when I use a real signal that is a smooth Gaussian area near the peak has 'digital noise'. I have confirmed by debug printouts the code is calculating what I want it to. It seems the combination of the input signal & the coefficients really does generate a non-monotonic output.

as reference I used the info at this site: [polyphase filter faq][1]

What am I missing? How do I avoid this? My coefficients are:
[4066   6261    8133    9391    9834    9391    8133    6261    4066]

The inputs sample rate is 100Mh and the output is 300Mh. I choose the filter coefficients to have a bandwidth < 50Mh.

Adding followup information. I looked at the step response and it oscillated every 3 output samples. The output to a step should be the sum of the coefficients in a branch. I noticed the 3 coefficients of the 3 branches did not all add to the same value. By deduction, I changed the coefficient 6261 to 5878. Now the coefficients are:
coeff = [4066   5878    8133    9391    9834    9391    8133    5878    4066] and the response is much better.

See new plot:

Thanks,
Brian

At Peter K's suggestion I am submitting my solution as the answer. In the poly phase filter the 9 coefficients were all symmetric about the center but gave poor ripple performance. When arranged in the format of a poly-phase filter the 3 coefficients in the 3 branches did not add up to the same value for all three branches. By making a slight edit to the 1 coefficient at +/-3 positons from the center all the branches added up to the same value and the results are as shown.

To get better performance then that shown I need to increase the # of coefficients to N=15 (wanting to keep N a multiple of L=3 and odd).

Brian


Can an Delta Dirac function (not response) ever be sampled? This is purely theoretical question but can an actual delta-dirac function ever be sampled? To make it clear again, I'm not talking about an impulse response from the impulse function. Lets say we have a truly instantaneous impulse as defined with delta-Dirac that the impulse with infinite(or finite) height and no width, an instantaneous pulse. This pulse is to be detected by sampling, but to find this instant signal, wouldn't this require an infinitely fast sampling or infinite sampling rate? Thinking this in terms of frequency domain, I believe delta function is 1 or has infinite spectrum. So would delta function ever be sampled unless we have infinite sampling rate?

A Dirac delta impulse cannot be sampled because it is a distribution and as such it does not make sense to talk about the value of $\delta(t)$ for e.g. $t=0$ because that value isn't defined. However, sampling extracts values of a function for certain values of $t$. The Dirac impulse is no function in the conventional sense. It's only under an integral that $\delta(t)$ makes sense.

Furthermore, the mathematical description of sampling is a multiplication with a Dirac comb. However, the expression $\delta^2(t)$ is undefined, so sampling a Dirac impulse can't be described mathematically in any meaningful way.

Note that the discrete-time unit impulse $\delta[n]$ can be sampled without any problem, because it has well-defined values for any $n$, unlike its continuous-time counterpart.


How can I calculate the amplitude of a signal from its power in dBm? $$\begin{align}
P_\text{out} &= 51.787\text{dBm}\\
P_\text{out} &= 10\log_{10}\frac{A_\text{out}^2}2\\
A_\text{out} &= \sqrt 2 \cdot 10^{\frac{P_\text{out}}{20}}\\
&= 549.3701
\end{align}$$

Is this calculation correct to yield the amplitude of my output signal when the output power in dBm is given? Or do I have to calculate power in dBW first? I assume an average output power of $P_{out} = A_{out}^2/2$

The Decibel Miliwatts Scale $dBm$ is the power ratio in Decibels, considering a reference of $P_0=1mW$.

$$P[dB]=10\text{log}_{10}(P/P_0)=10\text{log}_{10}(P[W]/0.001)$$

The general context is considering the signal amplitude $V$ in $volts$ and the power as an impedance load system $R$ in $\Omega$:

$$P=\frac{1}{R}V^2$$

Moreover, sinusoidal excitations are considered (for a given frequency also!), so we also can have:
$$P=\frac{1}{R}V_{rms}^2=\frac{1}{2}\frac{1}{R}V_{pk}^2$$

Thus:

$$
P[dBm] = 10\text{log}_{10}\left(\frac{V[V]^2}{R[\Omega]}\frac{1}{P_0}\right) = 20\text{log}_{10}\left(\frac{1}{\sqrt{R[\Omega]}}\frac{V[V]}{\sqrt{P_0}}\right)
$$

So, as reference:

$1W$ equals $30dBm$,
for a $50\Omega$ impedance system, this also equals $7.07V$, and
under the sinusoidal case, this also equals $10V_{pk}$ and $7.07V_{rms}$.


Filters - Hilbert Transform and Absolute Value of Magnitude Response Following relationship between magnitude response and phase response for minimum phase:

I have implemented in C++ code using the Hilbert Transform C code from file found online ht.c. However, I am not getting the phase response as expected. My magnitude response follows raised cosine formula:

Essentially I provide $|H(j\omega)| = \text{Gain}(f)$ to be all real values as per this cosine formula. The phase response is not as expected. Also, how am I suppose to handle $\log(|H(i\omega)|)$ when $|H(i\omega)| = 0$? Am I missing any steps to obtain phase response in this process?

Here is how phase and magnitude response should look like:

Neb, i dunno where you got that frequency response expression from, but it's flawed.  from the description (of parameters like $W$ in "octaves"), it should be like this:

$$ Gain(f) = \begin{cases}
10^{Boost/20} \times \frac{1}{2}\left\{1 + \cos\left(\pi \frac{\log_2(f/F_c)}{W} \right) \right\}, & \text{if }|\log_2(f/F_c)|<W \\
0, & \text{otherwise } \end{cases} $$


The inverse of an orthogonal matrix is its transpose In the following statement I don't understand the case for $\ i = j$:

Let $\mathbf A$ be an $\ m \times \ n$ orthogonal matrix where $\ a_i$ is the $\ i^{th}$ column vector. The $\ ij^{th} $ element of $\mathbf A^{T}\mathbf A$ is

$$
\left(\mathbf A^T \mathbf A\right)_{ij}  = a_i^{T} a_j =
\begin{cases}
1 &\text{if} &i = j\\
0 &\text{if} &\text{otherwise}\\
\end{cases}
$$

Therefore, because $\mathbf A^{T} \mathbf A = \mathbf I $, it follows that $\mathbf A^{-1} = \mathbf A^{T}$.

An orthogonal matrix has orthogal columns, i.e. the scalar product of two different columns is zero (the case $i\neq j$). For the case $i=j$ you have $a_i^Ta_i=\|a_i\|^2>0$.

So, all you can say about an orthogonal matrix with colums $a_i$ is:

$$
a_i^Ta_j=\begin{cases}c_i>0 & i=j\\ 0 & i\neq j\end{cases}$$

where $c_i>0$ is some constant, depending on the matrix.

Then, there are orthonormal matrices. These fulfill, in addition to the previous statement, $c_i=1$ and hence $A^TA=I$.

Finally, only if your $A$ is a square orthonormal matrix, it holds $A^{-1}=A^T$.

The task you describe has two flaws:

it only talks about ortogonal matrices, where we can just say $A^TA=\text{diag}(\vec{c})$ with $\vec{c}$ contains the norms of the columns.
It states $A$ is an $n\times m$ matrix (i.e. not square). Rectangular matrices dont have an inverse. So, the last statement is not correct.

Edit: Normally in literature, orthogonal matrices are considered to be square with unit-norm columns. Then $A^TA=I$ and $A^{-1}=A^T$. However, as the question explicitely mentions rectangular matrices, I want to point out that it really depends on your definition of what an orthogonal matrix is. And, if you just define it as orthogonal columns, the arguments above hold.


Neural Network: Spectrogram Dimension I like to work with a convolutional neural network in combination with a spectrogram as input.

Assuming the spectrogram has the dimensions $T\times F$ (time$\times$frequency). Is it more natural to see the spectrogram as a 1D image with $F$ channels and the size $T\times 1$ or as a 2D image with the size $T\times F$?

Depending on this I work with 1D convolutions or with 2D convolutions. But I'm not sure whats the more natural way.

Thank you very much

Depends on on what kind of feature(s) you want your CNN trained.  With 2D convolution, convolution in the F dimension of a T-F spectrogram will produce something like a real cepstrum, which has proved useful in speech recognition and musical pitch detection.


Auto/cross correlation of a cosine in the context of an adaptive beamformer I am currently reading chapter 13 of 'Adaptive Signal Processing' by Widrow and Stearns, so if anyone happens to have a copy of the book to hand it could be helpful.

I am reading the chapter on "Introduction to Adaptive Arrays and Adaptive Beamforming" and am struggling understanding the auto-correlation/cross correlation between the primary and reference signals at a delay of one (as shown in the diagram below) and in the more general case (where the signal is arriving at an angle), where the general cross correlation vector is described by:

$$
\mathbf{P} = \begin{bmatrix}
\phi_{dx}(0) \\ \phi_{dx}(1)
\end{bmatrix} = \begin{bmatrix}
E [C \cos k\omega_0 \cdot C \cos (k + \delta_0)\omega_0] \\ E [C \cos k\omega_0 \cdot C \sin (k + \delta_0)\omega_0]
\end{bmatrix}
$$

where $k$ is the sample number, $C$ is a constant amplitude, $\delta_0$ is the time delay of arrival between the primary and reference sensor and $\omega_0$ is the angular frequency.

The book assumes that the refence signal to the array is a cosine of the form:
$$
\text{reference signal} = C \cos [(k + \delta_0) \omega_0].
$$
(The above equation is 13.3 in the book).

My question is, I don't understand why, for $\phi_{dx} (1) $ the 'lag' part of the equation becomes a sine rather than a cosine with a shift 1 (sample)? It's as if there is a shift of $-\pi/2$, hence the sine, but I can't see where this is from? Why for a lag of 1 does this cosine seemingly change to a sine?

My understanding of cross (or auto) correlation is that for $\phi_{dx} (1)$ you compare one signal with another with a shift of 1 sample (hence a lag of 1).

I understand that a sine is the quadrature of a cosine and understand why these equal zero (as the book describes; "$\phi_{dx} (1) $, is zero because it represents the correlation of (13.3) with its quadrature (sine) component"). My question is why is a sine introduced for a lag of 1?

Diagram for completeness:

Many thanks

The $90$ degree shift introduces the sine component to the cross-correlation vector $P$ at lag 1 ($\phi_{dx}(1)$).

The block should be considered to be a 'high-level' description, and does not mean that a shift of one sample is equivalent to  a $90$ degree phase shift, or anything like that. Rather it should be thought of as "shift by as many samples as necessary for a phase shift of $90$ degrees".


Limiting bandwidth by adjusting sample rate I have an SDR that is capabple of 2.4Ms/s.  I am using it to measure the power of signals.  Right now I am setting the sample rate to 2M, collecting a burst of samples, and calculating the average power using

10 * log10(sum(abs(sample)**2 for sample in samples) / len(samples))

This gives me the average power over 2MHz, I think.  I'd like to limit that.  One approach is to use a filter.  But it would be simpler just to sample less frequently (sample at my target bandwidth).  But will the output measurements be as accurate?

Depends on your signal. Typically we try to sample at the Nyquist rate, which is equal to two times the maximum frequency of your signal. If you sample less often than this, you will lose information and your answer will not be correct, I think.


How to test if error has zero mean? I know how to calculate the mean of a data set. In DSP, I think zero-mean should means that the error fluctuates around zero, and the mean of errors at any interval is zero. So what about the error in the picture that I show below? The noise in the picture doesn't look like a zero-mean error at all. But the overall mean of all the error values is zero (near zero). How do I mathematically say that the error in the picture is not zero-mean? Usually, which statistical tool is used for testing if a error has zero-mean?

(I think my problem related to the interval at which the mean is calculated. With the error I show in the figure, the overall mean is zero, but the mean at smaller interval is not zero. But I don't know if the official definition about zero-mean care about such thing.)

Edit: to be more correct, I'm asking about zero-mean error, not noise. The figure I show below is the error of my measurement (only error, not include the true value). The error is not stationary, not ergodic. The mean and variance of error change over time but following a repetitive patterns. However, the mean of error over a long period is near zero.

Once you remove the systematic error (or trend) you want to test whether the process is zero mean or not. In hypotheses testing this would assume the following formulation:

$H_1: \quad x[n] = \theta + w[n], \quad n=0, \dots, N-1$,

$H_0: \quad x[n] = w[n], \quad n=0, \dots, N-1$,

where $w[n] \overset{i.i.d}{\sim} N\left ( 0, \sigma^2 \right )$, $\sigma^2$ is an unknown nuisance deterministic parameter and $\theta$ is another deterministic parameter that bares information on determining which hypothesis we choose.

We will develop a generalized Likelihood Ratio Test (GLRT) for this problem. In the GLRT framework the unknown deterministic parameters are estimated using maximum likelihood criterion, i.e.:

$\quad \widehat{\theta}_{ML} = \frac{1}{N}\sum\limits_{n=0}^{N-1} x[n]$, is the sample mean.
$\widehat{\sigma}^2_{ML} = \frac{1}{N} \sum\limits_{n=0}^{N-1} \left ( x[n] - \widehat{\theta}_{ML}  \right )^2$ is the sample variance under $H_1$, and,
$\widehat{\sigma}^2_{ML} = \frac{1}{N} \sum\limits_{n=0}^{N-1}  x^2[n]$, under $H_0$.

After some algebra dividing the two likelihood functions we get the following test:

\begin{equation}
Z = \frac{\widehat{\theta}_{ML}^2}{\hat{\sigma}^2_{ML}|H_0} = \frac{\left (  \frac{1}{N}\sum\limits_{n=0}^{N-1} x[n]\right )^2}{\frac{1}{N} \sum\limits_{n=0}^{N-1}  x^2[n]} \overset{\underset{>}{H_1}}{\underset{H_0}{<}} \gamma,
\end{equation}

where $\gamma$ is determined by the Neyman-Pearson's theorem given a desired probability of false alarm $P_{FA}$.


How does voice affect a sound wave or its corresponding waveform? Someone in the physics stack exchange told me to come here for this question. My question is how does a person's voice affect how a soundwave looks? For example, if I say the letter "a" it would sound different than if you said the letter "a", but my question is if there is a pattern to this, or maybe if we are able to derive a waveform based on voice (however it would be affected by multiple other factors due to when I say "a" it is different than when I say "b").

I would view it as the following:
The letter A is the input to a system, one is your vocal chords the other one mine. Since the system response of your and my vocal chords are different, we would get different outputs.
Now, what would characterize the system response, it would be factors like, rise time/bandwidth, decay response etc. All of this would impact the tone and pitch terminologies usually used in audio/speech. These variables will vary for each user. With different input the response will then also be different . So "b" will not be same as "a"


How can I approximate a signal I have a signal as the below figure. I am looking for some schemes to approximate the signal. I have found some ways for that task:

Using basis function (polynomial) to approximate the signal
Divide the signal to many subrange, and assume that the signal in the range is constant.

Could you suggest to me the other way to approximate the signal? Thanks

Assuming you want to remove the DC bias from your signal, a few thoughts come to mind. Your signal sequence looks like a noisy periodic wave riding on a DC bias. And the bias is trending upward in value from the beginning to the end of your sequence. Computing the average value of your signal sequence, a single value, and subtracting that value from your sequence will produce a sequence whose average value is very close to zero. But the upward bias trend will remain in your new sequence.

A process called "DC blocking" will remove both the average bias and the upward bias trend from your signal. If you have Matlab software, then use Matlab's detrend() command to implement DC blocking. If you don't have Matlab available then have a look at the following DC blocking networks (copied from Chapter 13 of my DSP book).

Other than different DC gains those above networks have the same frequency response. (Variable $\alpha$ is in the range of just greater than zero -to- one.) Note, those networks have very nonlinear phase response near zero Hz. If linear phase DC blocking is what you want then have a look at:
http://www.dsprelated.com/showarticle/58.php


What techniques/algorithms can I use to reduce noise when recording audio? I'm making a handmade guitar chords dataset by recording 20480 samples of a 44.1kHz microphone to later on make a chord classification using CNN.
But there's a problem... I'm making this dataset on my room. Quiet.
When other noise enviroments comes to the case the classifier will give poor results. I want to apply some technique to reduce the noise, even if little.
Do window functions reduce noise? Can someone help me by some advice of techniques/algorithms? Using Python, if possible...

It's probably better to train your classifier with recordings made in a very large variety of typical noisy situations, rather than to try to modify the input to match that with which the classifier was trained.


Bandwidth of Information Signal I have trouble finding the bandwidth of a signal. Say I have an info bearing signal m(t)=sinc(2t/pi). I found the fourier transform of the sinc function and found that the angular frequency was 1/pi. I am confused whether the the bandwidth is w or 2w. Since it is band-limited, does that mean only the positive w is counted and the bandwidth is just 1/pi?

When you look on a signal's DFT transform what you see is given in normalized units in the range $ [ -\pi, \pi ] $ where $ \pi $ equals the sampling rate.

Since the DFT is periodic and the way FFT work when one calculate the DFT in modern software (MATLAB / Python, etc...) the DFT is given in the $ [0, 2 \pi ] $ range.

In order to view it in the more intuitive way one could use fftshift() (MATLAB).

When you display it then, probably you won't see it gets to zero anywhere.
What you'll probably see (If the signal is band limited with bandwidth lower than the sampling rate) is a big drop in the magnitude.


RMS calculation in frequency domain after windowing I can calculate RMS in frequency domain as derived from Parseval's Theorem.

But what if I have applied a windowing function before doing the FFT (in my case a Hann window)?. Now the RMS values are not correct any more. How do I calculate the correct RMS values?

For a reasonably stationary time domain signal you can assume

rms(signal*windows) = rms(signal)*rms(window)

and can estimate the RMS of the un-windowed time signal this way. For a hanning window that means calculate the RMS using Perceval and then multiply with sqrt(8/3). It's not exact but will be fairly close if the signal is stationary with respect to the window length.
If it isn't stationary the interpretation of RMS isn't particularly well defined and it would helpful to learn more about your application and what exactly you are trying to do.


Compare the quality of IQ data from different Software Defined Radios (SDRs) I want to know what methods are best to scientifically evaluate the quality of different Software Defined Radios (SDRs) that spit out IQ samples.

For this question, I would like to set the boundary conditions as follows:

Their ADCs are from different vendors but contain the exact bit resolution.
They may have very different RF front ends or even be based on different architectures such as Direct conversion, Direct sampling or even Super-Heterodyne. We do not care about this as long as they tune to the band we specify and provide IQ.
During the test, they have same exact type of antenna plugged in and all other environment parameters are maintained same.
SDRs are tuned to the exact band to RX at an officially supported band in the specification.

After gathering, lets say, 2-3 IQ recordings (one from each SDR) through GNU Radio, how can we compare the quality of these IQ files?

Doubts: Can we tune the receivers to the same FM channel and ask some block in GNU-Radio to report the SNR so we can chose the one with highest SNR as the SDR that offered the most quality IQ?

Appreciate your time and knowledge.

First and foremost, I would recommend against over the air testing for this given the significant challenge in really being able to provide the same signal to each radio (since you have both temporal and spatial constraints that you cannot simultaneously meet).

I would instead use one GNU radio as a transmitter (or any other repeatable high quality source) to provide a consistent and repeatable waveform that can be used for each of the radios under test with a cabled test where the attenuation can be carefully and repeatably controlled over a wide range.

With this, there are several tests that immediately come to mind:

An ideal test waveform would be band limited complex Gaussian noise (pseudo-random noise as a repeatable pattern) since it will best occupy all possible positions on the complex IQ plane as well as all frequencies in the band of interest. This would then test the radio for a wide range of modulations, but if you have specific modulations in mind then those such waveforms for test purposes would also be ideal. (OFDM waveforms are well represented as complex Gaussian distributed bandlimited waveforms by the way, but it would be sufficient enough to simply create a complex random waveform and band-limit it for this generic case.)

An ideal test metric would be the correlation coefficient, $\rho$, as it will show errors from all sources including noise figure, quadrature error, LO phase noise, ADC quantization noise, IP3, IP2, etc.  For more information on that see Noise detection and How can I find SNR, PEAQ, and ODG values by comparing two audios?

The test should be done over a wide range of input power levels to show the related sensitivity and maximum power handling for comparative receivers. Thus the test metric for all receivers would be $\rho$ vs $P_{in}$.

The $\rho$ may be limited by the transmitter if that used is not of sufficient quality, so care must be taken to confirm that the $\rho$ for the transmit waveform is better than that which would result for all the tests performed (or represent the ceiling of the measurement). For the same reason start with the highest test power with known confirmed $\rho$ and adjust the power level in the cabled test receiver using an inline passive variable attenuator (with a rated much higher input P1dB) since any active device or variable amplifier is more likely to introduce additional distortion.

Additional tests involve how the receiver rejects interference such as intentional jammers or co-site transmitters. This is usually done with one tone testing in an adjacent band to help distinguish those receivers with better LO phase noise, and two tone interference testing in such a way to test for IP3 (and IP2 for zero-IF receivers in particular). The same test as above with the correlation coefficient vs input power is repeated, but in the presence of the jammers which can be done with unmodulated tones. Different recievers can have different sensitivities based on their architecture (or mis-managed spurs from a poor design), so a robust test would sweep the test jammers over a wide frequency range. If the radios are intended for specific applications, then the interference tests can be limited to cases appropriate for their known environment.

Instantaneous dynamic range in band would also be of interest.  As a simple test for this an in band single tone could be used at maximum input and the noise floor and spurious free dynamic range could be assessed in an FFT of the captured IQ signal. (You could also do specific linearity and single-tone and two-tone in band testing if you wanted to isolate these characteristics from what was otherwise in the catch-all $rho$ test described above).


Fourier Transform of the Hilbert Transform of cos(t) (using Fourier time-shifting property) If $x(t)=cos(t)=\frac{1}{2}e^{jt}+\frac{1}{2}e^{-jt}$, then $X(\omega)=\pi \delta(\omega-1)+\pi \delta(\omega+1)$. If $y(t)=cos(t-\frac{\pi}{2})=\frac{1}{2}e^{j(t-\frac{\pi}{2})}+\frac{1}{2}e^{-j(t-\frac{\pi}{2})}=\frac{1}{2}e^{jt-j\frac{\pi}{2}}+\frac{1}{2}e^{-jt+j\frac{\pi}{2}}$, then $Y(\omega)=\frac{e^{-j\frac{\pi}{2}}}{2}F\{e^{jt}\}+\frac{e^{j\frac{\pi}{2}}}{2}F\{e^{-jt}\}=-j\pi \delta(\omega-1)+j\pi \delta(\omega+1)$. $y(t)$ is the Hilbert transform of $x(t)$. This makes sense so far.

The trouble is when I try to evaluate $Y(\omega)$ using the time-shift property of the Fourier transform. $Y(\omega)=F\{x(t-\frac{\pi}{2})\}=X(\omega)e^{-j\omega\frac{\pi}{2}}=(\pi \delta(\omega-1)+\pi \delta(\omega+1))e^{-j\omega\frac{\pi}{2}}$. If I were to plot the phase of $Y(\omega)$, then we would see a line with slope $-\pi/2$. This does not seem to match the Hilbert Transform.

I forgot about the property of a dirac delta function: $f(t)\delta (t)=f(t)|_{t=0} \delta (t)=f(0)\delta (t)$

$Y(\omega)=F\{x(t-\frac{\pi}{2})\}=X(\omega)e^{-j\omega\frac{\pi}{2}}=(\pi \delta(\omega-1)+\pi \delta(\omega+1))e^{-j\omega\frac{\pi}{2}}=\pi \delta(\omega-1)e^{-j(1)\frac{\pi}{2}}+\pi \delta(\omega+1)e^{-j(-1)\frac{\pi}{2}}=-j\pi \delta(\omega-1)+j\pi \delta(\omega+1)$


CDMA and Pulse Shapping I am a little confused on CDMA and its interaction with the typical pulse shaping filter ex, raised root cosine.
Is the RRC filter bandwidth selected to be consistent with the CDMA spreading/update sequence? Or is it kept the same as prior to CDMA implementation so the number of pulses is consistent with the number of transmitted bits?

Pulse shape filtering is used to constrain the signal bandwidth (with what would otherwise be a Sinc function in frequency given the rectangular pulse) and as a Nyquist filter done in such a way as to not introduce inter-symbol interference.
For CDMA and specifically direct-sequence spread spectrum in many cases we have no concern for spectral containment (GPS for example since by the time the signal reaches any receivers it is well below the noise floor and therefore not an interference to other receivers in adjacent bands). But IF we did, then it would definitely be the individual chips that should be pulse shaped, since it is their base pattern that would define the spectral occupancy.


Conjugate symmetry of real-coefficient filters in Oppenheime's Discrete Time Signal Processing Is someone able to tell me where and how I could look up the yellow part in Oppenheime's Discrete Time Signal Processing, or explain it?

The Fourier transform of a real function $x[n]$ is conjugate symmetric. That is,
$$X^*(\omega)=X(-\omega)$$
This is easy to observe. From the definition of FT: $$X(\omega)=\sum_{n=-\infty}^{+\infty}x[n]e^{-j\omega n}$$
Conjugating both sides:

$$\begin{align}
X^*(\omega)=&\sum_{n=-\infty}^{+\infty}x^*[n]\left(e^{-j\omega n}\right)^*\\
&=\sum_{n=-\infty}^{+\infty}x[n]e^{j\omega n}\\
&=X(-\omega)
\end{align}$$

A direct consequence of this property (that you can easily verify) is that

  1- The real part and the magnitude of the Fourier transform of a real valued function are even functions of $\omega$. That is, $\Re\{X(\omega)\}=\Re\{X(-\omega)\}$ and $|X(\omega)|=|X(-\omega)|$.

  2- The imaginary part and the phase of the Fourier transform of a real valued function are odd functions of $\omega$. That is, $\Im\{X(\omega)\}=-\Im\{X(-\omega)\}$ and $\angle X(\omega)=-\angle X(-\omega)$.

Based on the above explanation, for a real-coefficient filter $h$, the magnitude of the frequency response is an even function.

$$\begin{align}
|H(e^{j(\omega-\pi)})|&=|H(e^{j(-\omega-\pi)})|\\
&=|H(e^{j(-\omega-\pi+2\pi)})|=|H(e^{j(\pi-\omega)})|
\end{align}$$
and the second line is because DTFT is $2\pi$-periodic. Hence we can add $2\pi$ and the result does not change.


What is the interpretation of the discrete-time spectrum? The CTFT of an analog signal is a representation of that analog signal in terms of the frequency parameter of sinusoidal (cosine specifically) functions whose weighted sum make up that signal. The cosines are continuous functions of time and the amplitude of a particular cosine at a given frequency is given by the amplitude of the spectrum at that frequency.

Is the interpretation for the DTFT of a sampled signal analogous where the spectrum indicates the frequency of (discrete-time/sampled) cosines that sum to the sampled signal, the spectrum amplitude indicates the magnitude of that cosine, etc.?

For instance, if $F=100\textrm{ Hz}$ and $|X(F=100\textrm{ Hz})| = 50$, then we know we have the component cosine signal: $50\cdot\cos(100n)$. Assuming $X(F)$ is the DTFT of $x(n)$.

Let's first consider the continuous-time case. From the formula for the inverse Fourier transform

$$x(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(\Omega)e^{j\Omega t}d\Omega\tag{1}$$

you can see that at a given frequency $\Omega$ (in rad/s), the value of the spectrum $X(\Omega)$ can be interpreted as the "strength" with which the complex exponential $e^{j\Omega t}$ contributes to the signal $x(t)$. You cannot conclude (as you seem to suggest) that the signal contains a sinusoid with amplitude $|X(\Omega)|$. This would only be the case if the signal $x(t)$ were periodic and, consequently, the integral in $(1)$ could be written as a sum, namely the Fourier series. Also note that even in that case the sinusoid is not just a cosine with zero phase, because the Fourier coefficients are generally complex and, hence, determine the amplitude as well as the phase of the respective sinusoid.

With this in mind, we can interpret the formula for the inverse discrete-time Fourier transform (IDTFT):

$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{jn\omega}d\omega\tag{2}$$

Note that I use a different variable here because $\omega$ in $(2)$ is not the same as $\Omega$ in $(1)$ (see below). Also note that the limits of the integral in $(2)$ are $\pm \pi$ instead of $\pm \infty$. The reason for the latter is that the sampling frequency sets an upper limit for the frequencies that can be uniquely represented by the DTFT (check out Aliasing and the Sampling Theorem if this is new to you). Apart from that, the structures of $(1)$ and $(2)$ are the same. So you can interpret the DTFT of $x[n]$ at a certain "frequency" $\omega$ as the "strength" with which the complex exponential $e^{jn\omega}$ is represented in the sequence $x[n]$. Again, it's generally not the amplitude of some sinusoid contained in $x[n]$, for the same reasons as outlined above.

The last thing that remains is the explanation of the difference between the frequency variables $\Omega$ and $\omega$. $\Omega$ is simply the frequency in radians per second, so we have $\Omega=2\pi f$, where $f$ is the frequency in Hertz. Since in discrete-time all frequencies must be seen as relative to the sampling frequency $f_s$, the variable $\omega$ is a normalized frequency in radians: $\omega=2\pi f/f_s$. So the integration limits in $(2)$ correspond to $f=\pm f_s/2$. The frequency $f_s/2$ is called the Nyquist frequency, and frequencies above the Nyquist frequency cannot be represented uniquely for a given sampling frequency.


Make a Movie from Multiple Images in MATLAB I have 10 images in the folder say "D:\images\" named as frame1,frame2,...,frame10. I want to make a movie from the color images with frame rate as per choice.Also i want to save created video file in .mp4 format at the location "D:\video\".

Is it possible to do so in MATLAB?how?

It is easy to do using the VideoWriter class.

You easily add frame (Image) using the proper method which is writeVideo.

After creating a video object, videoObject and given an image by mImage001, adding frame is done by:

writeVideo(videoObject, mImage001);

Enjoy...

Code Sample

Here is a code sample based on Mathworks documentation with some remarks of mine:

% Creating Video Writer Object
writerObj = VideoWriter('peaks.avi');
% Using the 'Open` method to open the file
open(writerObj);

% Creating a figure.
% Each frame will be a figure data
Z = peaks; surf(Z);
axis tight
set(gca,'nextplot','replacechildren');
set(gcf,'Renderer','zbuffer');

for k = 1:20
   surf(sin(2*pi*k/20)*Z,Z)
   % Frame includes image data
   frame = getframe;
   % Adding the frame to the video object using the 'writeVideo' method
   writeVideo(writerObj,frame);
end

% Closing the file and the object using the 'Close' method
close(writerObj);

This should be straight forward.


Designing butterworth filter with pole placement I am trying to design a Nth order continous time Butterworth filter in state space, with pole placement technique.

n = varargin{1}
Fc = this.fc
wc = 2*pi*Fc;
T = 1/(2*pi*Fc);
%----------------Finding Poles--------------------------
 p = exp(1i*(pi*(1:2:n-1)/(2*n) + pi/2)); %n is the order of the filter
                 p = [p; conj(p)];
                 p = p(:);
                 if rem(n,2)==1   % n is odd
                     p = [p; -1];
                 end
%----------------Finding Filter co efficents---------------
                 if length(p)<= 2
                    a = conv([1,-p(1)],[1,-p(2)]);
                end
                if length(p) > 2
                    a = conv([1,-p(1)],[1,-p(2)]);
                    for i = 3: length(p)
                        a = conv(a,[1,-p(i)]);
                    end
                end
                a = real(a)
                a = fliplr(a);
%----------------- calculating the matricies%-----------------

            %--------------A matrix---------------
            N = length(a)-1;% order of the matrix
            last_row = (a(1:end-1))*(-1/a(end));% creating the last row

             last_row = (wc.^(N:-1:1)).*last_row;% creating the last row

            k = ones(1,N-1);
            this.x = diag(k,1) ; % diag(vector, k) produce a matrix filled with zero's and k'th off diagonal as vector, in this  case, 1's.
            this.x(end,:) = last_row;% adding the last_row

            %-------------desiging B matrix--------------
            this.y = zeros(N,1); % creates a coloumn matrix with zeros
            this.y(N,:) = 1/(a(end)*wc^n); % writes last elements as a

            %--------------- writing C matrix--------------
            this.z = zeros(1,N); %creates a row matrix
            this.z(1,1) = 1;

            %------------ writing the D matrix----------------

            this.u = [0];

I know the code is bit lengthy, but i have tried everything, but can't figure out where i am going wrong. Because I am not getting any error. The resultant output is just 0. The code is pretty simple, I find the poles, then the filter coefficients. with the filter coefficients. i design A,B;C,D Matrices.
This is a continuous butterworth filter, Can anyone point the mistake i am doing?

edit1:
A low pass PT2 filter:

A low pass 2nd order butterworth filter:

I think I found the problem. Everything was fine except for the denormalization step for arbitrary cut-off frequencies wc. I changed the following lines to obtain the correct state-space matrices (and to handle the case n=1):

...
if (n == 1)
    a = [1,1];
else
    p = ...
    ...
end

a = real(a);
a = a./(wc.^(n:-1:0))   % new!
a = fliplr(a);
...

% last_row = (wc.^(N:-1:1)).*last_row;    % removed!

if( n == 1)
    this.x = last_row;
else
    k = ones(1,N-1);
    this.x = diag(k,1) ;
    this.x(end,:) = last_row;
end

...

this.y(N,:) = 1/a(end);    % changed!

...

I used the function ss2tf to compute the polynomial coefficients from the state-space matrices computed above, and the coefficients were correct for different values of wc. This means that the matrices are correct. If the output signal is still different from what you expect, the problem must be in the filtering routine.


Filtering out unwanted square wave (Radio: DCS/DPL signal) I've need to remove a DPL/DCS signal from a recording. Roughly speaking, a DCS signal is a low frequency (67.15Hz) square wave sent at the same time, over the same carrier, as speech. Because it is a square wave, it has many strong harmonics that overlap with speech. Obviously, the speech must be preserved.

On the surface, this seems like a problem that might be solved the same way as removing 60Hz power-line noise: lots of notch filters. However, power-line noise tends to be weaker and comes from a source that is much closer to a sine wave.

So, my question is: is there a better way to do this?

This link contains more info about DCS: http://onfreq.com/syntorx/dcs.html
It mentions "Since DCS creates audio harmonics well above 300 Hz (i.e. into the audible portion of the band from 300 to 3000 Hz), radios must have good filters to remove the unwanted DCS noise." Ha!

EDIT: DCS is not a "square wave". I think the correct way to explain it is to say it's a NRZ binary code clocked at 67.15 Hz.

If you know that it is a square wave, and the frequency is roughly constant, I would think about synthesizing an inverted copy of it and adding it to the signal.  The remainder would be your speech, plus some residual harmonics of the square wave frequency, which might be dealt with by a 1-channel adaptive filter (see Sec. IIC of this link: http://www.cs.cmu.edu/~aarti/pubs/ANC.pdf).  It's possible that the adaptive filter will do it all by itself--I have used them with good results on periodic interfering signals with attention to the delay setting involved.
In order to get the parameters of the inverted signal, you would do one of two things:
1. Fit a two-parameter square wave to the signal using least squares minimization (two parameters being the initial time delay and the amplitude, or add other parameters as you see fit).
2. Use an FFT to estimate phase and amplitude of the fundamental sine wave, then use the theoretical ratio of the amplitude of a square wave to its fundamental (isn't it 2/pi?) to reconstruct the square wave from that.

Which of these two approaches you pick depends on the bandwidth of the speech and its amplitude relative to your DCS signal.  Good Luck!


qam demapper why are soft-outputs used? I'm starting to learn about OFDM QAM demapping, and a lot of documentation points to using "soft-output" algorithms that calculate Log-Likelihood Ratios (LLRs).

Why is this method advantageous?

My naive approach would be to just place a grid over the IQ map, and find which grid-square my symbol lands in.

My naive approach would be to just place a grid over the IQ map, and find which grid-square my symbol lands in.

And that's (more or less) exactly how it's done: it's called a hard decision, and your  grid is what we'd call the decision boundaries. It's a result of applying the Maximum Likelihood criterion to the decision process.

Of course, in that decision, a sample that hits a constellation point (nearly) exactly is not distinguishable from a sample that was somewhere "in the middle" between two constellation points (and hence close to one of the decision boundaries).

You can't tell "absolutely certain" from "I'm about 51% sure" decisions that way! So, this hard deciding throws away valuable information.

If you find a way to encode the decision information and the certainty of that decsion into numbers, then you get a *soft-output decider". The most common way of finding a single number that describes that is through likelihood ratios, i.e. divide how likely one symbol is by the likelihood of a different symbol. Because ratios are unelegant to deal with, apply a logarithm to it, and get a LLR.

Now, you want bits, not something like "I'm almost certain this was a 0 or 1", but don't forget:

In modern systems, you use channel coding. If you have an algorithm that doesn't only consider single symbols, but sequences of symbols, then you can use these soft-outputs to infer the likeliest sequence that was sent more reliably.

As said in my comment to your question, you really shouldn't start learning about QAM mapping by reading about soft-output decision on QAM. Start with soft-output BPSK, and then learn about the hard- and soft-output Viterbi sequence estimator. (For Viterbi, you'll need working knowledge of simple convolutional channel codes, but that should be comparatively simple when considering you're doing LLRs on QAM...)


Understanding the Matched Filter I have a question about matched filtering.  Does the matched filter maximise the SNR at the moment of decision only?  As far as I understand, if you put, say, NRZ through a matched filter, the SNR will be maximised at the decision point only and that is the advantage of the matched filter.  Does it maximise the SNR anywhere else in the output function, or just at the point of decision?

According to Wikipedia

  The matched filter is the optimal linear filter for maximizing the signal to noise ratio (SNR) in the presence of additive stochastic noise

This to me implies that it maximises it everywhere, but I don't see how that is possible.  I've looked at the maths in my communications engineering textbooks, and from what I can tell, it's just at the decision point.

Another question I have is, why not make a filter that makes a really tall skinny spike at the point of decision.  Wouldn't that make the SNR even better?

Thanks.

Edit:
I guess what I'm also thinking is, say you have a some NRZ data and you use a matched filter, the matched filter could be implemented with an I&D (integrate and dump).  The I&D will basically ramp up until it gets to the sampling time and the idea is that one samples at the peak of the I&D because at that point, the SNR is a maximum.  What I don't get is, why not create a filter that double integrates it or something like that, that way, you'd have a squared increase (rather than a ramp) and the point at which you sample would be even higher up and from what I can tell, more likely to be interpreted correctly by the decision circuit (and give a lower Pe (probability of error))?

Since this question has multiple sub-questions in edits, comments on answers, etc., and these have not been addressed, here goes.

  Matched filters

Consider a finite-energy signal $s(t)$ that is the input to a (linear
time-invariant BIBO-stable) filter with impulse response $h(t)$, transfer function $H(f)$,
and produces the output
signal
$$y(\tau) = \int_{-\infty}^\infty s(\tau-t)h(t)\,\mathrm dt.\tag{1}$$
What choice of $h(t)$ will produce a maximum response at a given time
$t_0$? That is, we are looking for a filter such that the global maximum
of $y(\tau)$ occurs at $t_0$. This really is a very loosely phrased
(and really unanswerable) question because clearly the filter
with impulse response $2h(t)$ will have larger response than
the filter with impulse response $h(t)$, and so there is
no such thing as the filter that maximizes the response.
So, rather than compare apples and oranges, let us include the
constraint that we seek the filter that maximizes $y(t_0)$ subject
to the impulse response having a fixed energy, for example, subject to
$$\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt = \mathbb E
= \int_{-\infty}^\infty |s(t)|^2 \,\mathrm dt.\tag{2}$$

  Here onwards, "filter" shall mean a linear time-invariant filter
  whose impulse response satisfies (2).

The Cauchy-Schwarz inequality provides an answer to this question. We have
$$y(t_0) = \int_{-\infty}^\infty s(t_0-t)h(t)\,\mathrm dt
\leq \sqrt{\int_{-\infty}^\infty |s(t_0-t)|^2 \,\mathrm dt}
\sqrt{\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}
= \mathbb E$$
with equality occurring if $h(t) = \lambda s(t_0-t)$ with $\lambda > 0$
where from (2) we get that $\lambda = 1$, that
is, the filter with impulse response $h(t) = s(t_0-t)$ produces
the maximal response $y(t_0) = \mathbb E$ at the specified time $t_0$.
In the (non-stochastic) sense described above, this filter is
said to be

  the filter matched to $s(t)$ at time $t_0$ or
  the matched filter for $s(t)$ at time $t_0.$

There are several points worth noting about this result.

The output of the matched filter has a
unique global maximum value of $\mathbb E$ at $t_0$; for any other
$t$, we have $y(t) < y(t_0) = \mathbb E$.
The impulse response $s(t_0-t) = s(-(t-t_0))$
of the matched filter for time $t_0$ is just $s(t)$ "reversed in time"
and moved to the right by $t_0$.

a. If $s(t)$ has finite support, say, $[0,T]$, then the matched filter is
noncausal if $t_0 < T$.

b. The filter matched to $s(t)$ at time $t_1 > t_0$ is just the filter
matched at time $t_0$ with an additional delay of $t_1-t_0$. For this
reason, some people call the filter with impulse response $s(-t)$,
(that is, the filter matched to $s(t)$ at $t=0$) the matched filter for $s(t)$ with the
understanding that the exact time of match can be incorporated into
the discussion as and when needed. If $s(t) = 0$ for $t < 0$, then
the matched filter is noncausal. With this, we can rephrase 1. as
The matched filter for $s(t)$ produces a unique global maximum
value $y(0) = \mathbb E$ at time $t=0$. Furthermore,
$$y(t) = \int_{-\infty}^\infty s(t-\tau)s(-\tau)\,\mathrm d\tau
= \int_{-\infty}^\infty s(\tau-t)s(\tau)\,\mathrm d\tau = R_s(t)$$
is the autocorrelation function of the signal $s(t)$. It is
well-known, of course, that $R_s(t)$ is an even function of $t$
with a unique peak at the origin.  Note that the output of the
filter matched at time $t_0$ is just $R_s(t-t_0)$, the autocorrelation
function delayed to peak at time $t_0$.
No filter other than the
matched filter for time $t_0$ can produce an output as large
as $\mathbb E$ at $t_0$. However, for any $t_0$,
it is possible to find filters that
have outputs that exceed $R_s(t_0)$ at $t_0$. Note that $R_s(t_0) < \mathbb E$.
The transfer function of the matched filter is $H(f)=S^*(f)$, the
complex conjugate of the spectrum of $S(f)$.
Thus, $Y(f) = \mathfrak F[y(t)]= |S(f)|^2$.
Think of this result as follows. Since $x^2 > x$ for $x > 1$ and $x^2< x$ for
$0 < x < 1$, the matched filter has low gain at those frequencies where
$S(f)$ is small, and high gain at those frequencies where $S(f)$ is large.
Thus, the matched filter is reducing the weak spectral components
and enhancing the strong spectral components in $S(f)$. (It is also
doing phase compensation to adjust all the "sinusoids" so that
they all peak at $t=0$).

-------

But what about noise and SNR and stuff like that which is what the OP
was asking about?

If the signal $s(t)$ plus additive white Gaussian noise with
two-sided power spectral density $\frac{N_0}{2}$ is processed
through a filter with impulse response $h(t)$, then the output
noise process is a  zero-mean stationary Gaussian process with
autocorrelation function $\frac{N_0}{2}R_s(t)$. Thus, the
variance is
$$\sigma^2 = \frac{N_0}{2} R_s(0) = \frac{N_0}{2}\int_{-\infty}^{\infty} |h(t)|^2\,\mathrm dt.$$
It is important to note that the variance is the same regardless
of when we sample the filter output.  So, what choice of $h(t)$
will maximize the SNR $y(t_0)/\sigma$ at time $t_0$?  Well, from the
Cauchy-Schwarz inequality, we have
$$\text{SNR} = \frac{y(t_0)}{\sigma}
= \frac{\int_{-\infty}^\infty s(t_0-t)h(t)\,\mathrm dt}{\sqrt{\frac{N_0}{2}\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}}
\leq \frac{\sqrt{\int_{-\infty}^\infty |s(t_0-t)|^2 \,\mathrm dt}
\sqrt{\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}}{\sqrt{\frac{N_0}{2}\int_{-\infty}^\infty |h(t)|^2\,\mathrm dt}} = \sqrt{\frac{2\mathbb E}{N_0}}$$
with equality exactly when $h(t) = s(t_0-t)$, the filter that is matched
to $s(t)$ at time $t_0$!! Note that $\sigma^2 = \mathbb EN_0/2$.
If we use this matched filter for our desired sample time, then at other
times $t_1$, the SNR will be
$y(t_1)/\sigma < y(t_0)/\sigma = \sqrt{\frac{2\mathbb E}{N_0}}$.  Could
another filter give a larger SNR at time $t_1$? Sure, because $\sigma$
is the same for all filters under consideration, and we have noted above that
it is possible to have a signal output larger than $y(t_1)$ at time
$t_1$ by use of a different non-matched filter.

In short,

"does the matched filter maximize the SNR only at the sampling
instant, or everywhere?" has the answer that the SNR is maximized only
at the sampling instant $t_0$. At other times, other filters could give
a larger SNR than what the matched filter is providing at time $t_1$,
but this still smaller than the SNR $\sqrt{\frac{2\mathbb E}{N_0}}$
that the matched filter is giving you at $t_0$, and if desired,
the matched filter could be redesigned to produce its peak at time
$t_1$ instead of $t_0$.
"why not make a filter that makes a really tall skinny spike at the point of decision. Wouldn't that make the SNR even better?"
The matched filter does produce a spike of sorts at the sampling time
but it is constrained by the shape of the autocorrelation function. Any
other filter that you can devise to produce a tall skinny (time-domain)
spike is not a matched filter and so will not give you the largest possible
SNR. Note that increasing the amplitude of the filter impulse response
(or using a time-varying filter that boosts the gain at the time
of sampling) does not change the SNR since both the signal and the noise standard deviation increase proportionately.
"The I&D will basically ramp up until it gets to the sampling time and the idea is that one samples at the peak of the I&D because at that point, the SNR is a maximum."
For NRZ data and rectangular pulses, the matched filter impulse response is
also a rectangular pulse. The integrate-and-dump circuit is a correlator
whose output equals the matched filter output only at the sampling instants, and not in-between.  See the figure below.

If you sample the correlator output at other times,
you get noise with smaller variance but you can't simply add up the samples
of I&D output taken at different times because the noise variables are highly correlated, and
the net variance works out to be much larger.  Nor should you expect to be able
to take multiple samples from the matched filter output and combine them
in any way to get a better SNR.  It doesn't work. What you have in effect
is a different filter, and you cannot do better than the (linear)
matched filter in Gaussian noise; no nonlinear processing  will give
a smaller error probability than the matched fiter.


Matlab: Different 1-D FFT results with and without zero-padded Recently I am implementing FFT with C code. So I use Matlab for algorithm verification. As Radix-2 FFT was chosen, I've learned that zero-padding technique could help to reach 2^N limitation without messing up the result. However, hesitation has shown.

Let's say a set of acceleration data with the size of 150 (x) comes across. If FFT function is called directly like:

X = fft(x);
abs_X = abs(X);
figure
plot(abs_X(2:75)); % discard first data (DC)

The result looks reasonable:

With this result I can easily indicate the highest amplitude and get the main frequency with frequency step. However, after zero-padding the data to the size of 256 and then apply fft, a confusing result is obtained.

x = [x; zeros(106,1)]; % pad 0 to 256
X = fft(x);
abs_X = abs(X);
figure
plot(abs_X(2:128)); % discard first data (DC)

Now I cannot identify the main frequency since the decreasing pattern is shown.

From here I have learned that the fft function in Matlab will automatically adjust the data size with zero-padding in order to perform FFT. So theoretically two operations demonstrated above should perform same results. I know that zero-padding will increase the visual resolution (not data resolution). But the difference between the beginning of the figures, where main frequency can't be determined, is confusing me.

Can anyone help to clarify the situation? Or do I misunderstand something important about the statement "Zero padding is going to perform a same result."? Any comment is appreciated. Thank you for browsing my question!

I don't know how MATLAB handles this internally. What will help in your situation however is multiplying your unpadded signal by a suitable window function, like a hann window.

What you see in your diagrams is a convolution with a $\text{si}$ function, which is the fourier transform of a rectangle function. If you just take the signal as it is and transform it, you essentially multiply it by a rectangle function.

If you try this, for both of your cases, you should get matching results.


Minimization of Essential matrix A problem in computer vision and 3d reconstruction is getting the camera's intrinsics  parameters. A common solution is to use an object in which one knows the measurements of the shape before hand, such as a checker board. The issue with this method is that it must be done every time the camera parameter's change such as focal length and magnification.

I'm trying to implement the camera self-calibration discussed in A Simple Technique for Self-Calibration. The essential matrix is constrained by its two singular values. This can be used to recover the camera's intrinsics without doing a manual calibration (i.e., with a checkerboard). I am a little confused by how the cost function can be minimized. Here's what I understand so far:

essential matrix
$$E=K_2^TFK_1$$

intrinsic matrix
$$K=\begin{bmatrix}\alpha_x & s & u_0 \\ 0 & \epsilon\alpha_x & v_0 \\ 0 & 0 & 1\end{bmatrix}$$

$\alpha_x$ product of focal length and magnification factor [solve]
$\epsilon$ aspect ratio [assume provided, I guess from camera or EXIF data?]
$u_0 v_0$ are the coordinates of the principal point [assume 0, 0]
$s$ skew [assume 0]

cost function
$$C(K_i,i=1..n)=\sum_{ij}^n(\sigma1_{ij}-\sigma2_{ij})/\sigma1_{ij}$$
the $\sigma$s are the singular values of $K_j^TF_{ij}K_j$

Question: How is this cost function being minimized?

I guess this is a straightforward non-linear optimization problem (to be solved with Newton variations, such as Trust-Region methods), where you don't even need to compute the Jacobian analytically. It appears to me that the optimization problem is written over $K_i$, and thus is the input to the cost function. To compute the cost, at each call to this function, you basically compute the singular values of $K_i$ and compute the cost according to equation $4$ (in the paper). As your input parameter is $K$, the derivatives are computed over the elements of $K$. That makes your optimization $4$ dimensional (or $5$ if you consider the skew) per camera. The derivatives are computed automatically, and you don't need to care about that. If you are using MATLAB,lsqnonlin would work for you.

The weight computation is explained in a detailed fashion in the paper, so I skip this part.

Having checked the paper once again, I have noticed that the authors are actually using the numerical differentiation scheme I have mentioned. If you would like to understand more deeply how to differentiate an SVD, you might like to check out this one or this.


Why is the size of the output of the FFT of a signal is same as the size of the input signal? I am kind of new to the DSP domain. I was trying to get the frequencies associated with a signal by performing FFT over it. I used numpy.fft.fft for this.

Input signal frequency = 1000 Hz

Noise signal frequency = 200 Hz

The input signal, noise signal and combined_signal all have 2400 samples. When I apply FFT on the combined_signal, i get 2400 size array as output. According to my understanding the FFT operation returns the frequencies present in the signal (in this case, 1000 Hz and 200 Hz). Also, when I plot the output of the FFT, I get some weird plot. Instead of giving two peaks at 1000 Hz and 200 Hz it is giving some random frequencies. Also, why is it giving are two peeks at the beginning and 2 at the end ?

I have used the following jupyter notebook code.

import numpy as np
from matplotlib import pyplot as plt
import matplotlib
%matplotlib inline
matplotlib.rcParams['figure.figsize'] = (16, 9)

freq_noise = 200
freq_signal = 1000

sampling_rate = 48000.0
seconds = 0.05
num_samples = int(sampling_rate * seconds)

signal = [np.sin(2 * np.pi * freq_signal * x1/sampling_rate) for x1 in range(num_samples)]
noise = [np.sin(2 * np.pi * freq_noise * x1/sampling_rate) for x1 in range(num_samples)]

# convert to numpy array
signal = np.array(signal) * 1000      # Increase the amplitude of the signal
noise = np.array(noise) * 500         # Lower amplitude for the noise

combined_signal = signal + noise

print "The length of the signal is : ", len(signal)
print "The length of the noise is : ", len(noise)
print "The length of the combined signal  is : ", len(combined_signal)

plt.subplot(4, 1, 1)
plt.plot(signal)
plt.xlabel('Sample number')
plt.ylabel('Amplitude')
plt.title('Original Signal')

plt.subplot(4, 1, 2)
plt.plot(noise)
plt.xlabel('Sample number')
plt.ylabel('Amplitude')
plt.title('Noise Signal')

plt.subplot(4, 1, 3)
plt.plot(combined_signal)
plt.xlabel('Sample number')
plt.ylabel('Amplitude')
plt.title('Combined Signal')

combined_fft = np.fft.fft(combined_signal)
combined_freq = np.abs(combined_fft)

plt.subplot(4, 1, 4)
plt.plot(combined_freq)
plt.xlabel('Frequency')
plt.ylabel('Amplitude')
plt.title('FFT of Combined signal')

plt.tight_layout()
plt.show()

The output of the above code is:

The length of the signal is :  2400
The length of the noise is :  2400
The length of the combined signal  is :  2400

To answer the question in your title:

Because it is defined to do exactly that.

Think about it: the DFT is a revertible linear operation on all vectors of the $\mathbb C^N$. In other words, for every vector that you apply the DFT to, there is exactly one vector that is the result, and for every result vector, there is exactly one vector that maps to it.

The operation can only exist on sets of the same dimensionality. And thus, the input and output vectors need to be of the same size.

The fact that you can interpret the DFT of a time sample vector as a frequency sample vector doesn't change that.

Your "why are there two peaks" question can be answered by: because your signal is real valued; real-valued signals always have symmetry in Fourier domain.

The "odd" shape is probably Gibb's/leakage phenomena. That always happens when your signal is not perfectly periodic within a circular repetition of your vector.


How to set borders for image translation What i am trying to do is to build a function similar to how imtranslate works in matlab.What I have until now thanks to FAT32's answer, is a function for translation which have the following effects:

Since the limits of the image are reached when it is translated from left to right, we can see that the translated pixels are in the left side.
What should be done to prevent such phenomenon?

What you observe in the image is the effect of the circular shift implied by the multiplication of the DFT of the image by the complex exponentials; i.e., the last (right most) sample becomes the first (left most) sample, hence circularly shift around the image, after a single pixel shift to the image is performed.

Theoretically you cannot prevent the shift phenomenon if you use the same method to translate the image; i.e., by using the inverse DFT of the exponentially multiplied forward image DFT.

However you can practically make it unobservable if you pad the original image by enough number of zeros at the ends so that the shifts of zeros will be unnoticed and only the linear translation of the image will be apparent.

To do so requires the amount of shift to be known in advance. Then given the shift amount and direction, you must pad that number of zeros in the correct end of the image so that when those zeros are circulated to the opposite end, they will not create a noticeable effect and the image will seem to be linearly shifted all.


Why square-root of Hanning Window I just studied an implementation (Matlab) of a windowed FFT and saw the following line of code:

w = sqrt(hann(fft_length, 'periodic'))

and later:

spectogram(:, column) = fft(ringbuffer.*w, fft_length) %why .*w ?

The ringbuffer obviously contains the current frame of the signal and the overlap so its's obvious to me that this needs to be multiplied with the Hanning Window  but why is the square-root being taken first?

The purists will tell you it should be called a Hann window, or VonHann, but not Hanning.

From a math perspective, the Hann window is also equivalent to a sine squared window at half the frequency, so by taking the square root you are getting a sine window.

Why use a window at all?  The seminal paper on the topic is:

harris, fred. (Jan 1978). "On the use of Windows for Harmonic Analysis with the Discrete Fourier Transform"

A link to the paper can be found in reference [10] of the Wikipedia article on windows.

Hope this helps.

Ced


Applying Nyquist's sampling theorem to a real signal I'm struggling to fully understand the Nyquist-Shannon sampling theorem.

For some message input signal $m(t)$ that is infinite in time (i.e. is not identically $0$ for any interval $t_1<t<\infty$), and that is hence necessarily band-limited (by the properties of the fourier transform), the sampling theorem states that if $f_{sampling-frequency}>2f_{maximum-fourier-component}$, then the original signal $m(t)$ is fully determined by the samples.

However, for any real signal that is transient and aperiodic, we will always take a set of samples in some range $t_1<t<t_2$, and hence the message is time-limited, excluding the possibility of it being band-limited.

For example, here it states:

  The Shannon’s sampling theorem was derived
  using the assumption that the signals must exist over
  infinite time interval. But all of our applications are
  based on finite time intervals. The objective of this
  research is to correct this inconsistency. In this paper
  we show where and how this infinite time assumption
  was used in the derivation of the original sampling
  theorem and then we extend the results to finite time
  case. Our research shows that higher sample rate is
  necessary to recover finite duration signals. This
  paper validates, with detailed theory, the common
  industrial practice of higher sample rate. We use the
  infinite dimensionality property of function space as
  the basis of our theories. A graphical example
  illustrates the problem and the solution.

However, I don't understand why this should be the case. If we have some $m(t)$ defined in the range $t_1<t<t_2$ then we can postulate an oscillatory function that extends to infinity that is made up of direct repetitions of the waveform with period $t_2-t_1$. Now, the fourier transform of the function defined over infinity is a set of discrete values that range over some bound (band-limited). For this infinite function, the sampling theorem should hold in its current form, and the Discrete Fourier transform (DFT) translates between the set of discrete time-domain samples ${m_n}$ and the discrete set of complex frequency domain samples $M_b$. As long as we specify the time range over which our original signal was taken, since the DFT does not specify whether the original samples was actually periodic or not), it would seem that Nyquist's theorem holds in its current form.

To summarise, any signal we take is not a time-limited function if we hypothetically extend it in a periodic manner. However, I know there is something wrong with my understanding - it would imply that for any real signal or use of the theorem, we only need consider the DFT, and if so why is it routinely derived for infinite functions at all? What am I missing i.e. how is what I think incorrect, or how is it in fact consistent with the quote above?

I do not see specifically what you think is incorrect. I follow your thinking and do not see a flaw but let's review the key characteristics below in case that fills in a hole in the thought process:

For periodic systems we can consider a time span that extends over one period (0 to T) or from $-\infty$ to $+\infty$. Both are equivalent as long as we agree that the $\infty$ case must be periodic.

Similarly for sampled systems we can consider a frequency span that extends over one frequency range, 0 to $F_s$, where $F_s$ represents the sampling frequency, or from $-\infty$ to $+\infty$. Both are equivalent as long as we agree that the $\infty$ case must be similarly periodic in frequency.  (Note the base frequency that repeats can also be given as $-F_s/2$ to $+F_s/2$ but the range 0 to $F_s$ is most common when referring to the DFT.)

So the DFT itself is both periodic in time and in frequency, representing a sampled system that is periodic in time.

For any other systems; aperiodic or not sampled, we cannot use the DFT and call that equivalent, we can only approximate those systems with the DFT.

Below is a summary of the various Fourier Transforms that we can use in those other cases:

CASE 1: Continuous and Aperiodic in time: The Fourier Transform

First case we have is the Fourier Transform, where the time function is by definition continuous and extending from $-\infty$ to $+\infty$. Since it is continuous in time (not sampled) the frequency result will be aperiodic (not repeating). Since it is aperiodic in time, the frequency result will be continuous.

CASE 2: Continuous and periodic in time: The Fourier Series Expansion

Note that x(t) can still be periodic in the case above; there is nothing that excludes taking the Fourier Transform as written in the graphic above on a periodic function. This will result in the continuous Fourier Transform result as being 0 except for the frequencies at multiples of 1/T, including zero, which may have non-zero values at only these locations.  Knowing this, we can solve for just those frequencies in this case, resulting in the Fourier Series Expansion:

CASE 3: Discrete and aperiodic in time: The Discrete Time Fourier Transform

The third case is to now consider the sampled version of our first case, where we have an aperiodic time domain function again, but this time it is sampled in time. This is the Discrete Time Fourier Transform. We have a similar result to the first case, the Fourier Transform, with the only change that it is now repeating in frequency. As previously described the repetition boundaries occur at multiples of $F_s =  1/T$, so if our frequency spectrum in our first case extended beyond $F_s$, we would get frequency folding (aliasing) as a result and therefore would not be able to represent our time domain waveform exactly (the sampling rate was not high enough!).

(Side note: Zero padding prior to taking the DFT approximates the DTFT, the more zeros you add in, the more the frequency result approximates a continuous function, and that function is samples of the DTFT!).

CASE 4: Discrete and periodic in time: The Discrete Fourier Transform

Finally, we come to the DFT (the FFT is an algorithm to compute the DFT exactly). In this case, when considering the same time axis as above that extends from $-\infty$ to $+\infty$, the time domain waveform must be repeating (periodic).  The periodicity in time results in discrete frequencies that can only exist at integer multiples of 1/T Hz where T represents one period in time in seconds. Note that like the Fourier Series Expansion, the periodicity is implied; the actual transform considers N samples in time and provides N samples in frequency, but the periodicity described is evident when both axis are extended to infinity.

I hope this answered your question. If not please comment on the remaining points of confusion and I will try and update accordingly.


Why does downsampling stretch a signals frequency response and upsampling shrink and create images of a signals frequency response? I am learning some basic DSP and I have a pretty good intuition as to why sampling creates spectral images of the frequency response at intervals of the sampling frequency (convolution with pulse train).

I am reading about sample rate conversions (downsampling/decimation and upsampling/interpolation) and I can't seem to get intuition as to why downsampling/decimation stretches the frequency response by M (the downsample factor) or why upsampling/interpolation compresses it by I (the upsample factor) and creates more images of the originally sampled datas frequency response

I try to think of it this way: In the time domain, when downsampling occurs, the signal gets compressed; while on upsampling, the signal gets stretched.

Then, from the Fourier transform we know that time stretching means frequency compression, and vice-versa.

It may not be a rigorous answer, but hope this helps!


Variance of Integral of a real white Gaussian Noise Process In this question, is the answer not equal to infinity ? Answer is mentioned as 6. But my doubt is cant we think of it like a linear combination of many independent random variables each having infinite variance, so the resulting random variable also has infinite variance. This was a question asked in GATE  exam conducted in India for Electronics and Communication stream

Since $W(t)$ is assumed to be zero-mean, also the RV $Y$ is zero-mean. Hence, the variance of $Y$ is given by
$$\begin{align}\sigma_Y^2&=E\left\{Y^2\right\}\\&=E\left\{\int_{-\infty}^{\infty}W(t_1)\phi(t_1)dt_1\int_{-\infty}^{\infty}W(t_2)\phi(t_2)dt_2\right\}\\&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\phi(t_1)\phi(t_2)E\big\{W(t_1)W(t_2)\big\}dt_1dt_2\tag{1}\end{align}$$
where $E\big\{W(t_1)W(t_2)\big\}$ is the auto-correlation function $R_W(t_2-t_1)$ of $W(t)$. Now you just have to figure out the expression for $R_W(\tau)$ and solve the integral $(1)$.


Rollercoaster: determining subject position from inertial data? I have inertial data (accelerometer, gyroscope, magnetometer) from phones of human subjects riding a single rollercoaster repeatedly.

The ride has temporal variance (ex. sometimes getting from A to B takes 300ms longer)
Phones may have different initial orientations
Phones stay in the subjects' pockets, but jostle a bit

Using this data, I want to

Learn the "mean" rollercoaster ride

as a time series
as a sequence of significant rollercoaster events

On a new ride, in real-time, robustly determine where the subject is along the path so we can align music to their experience \m/
Bonus: 3D model the rollercoaster path

How do I do this?

  My guess:

  Low-pass-filter data to remove jostle/noise
  Something to normalize the different initial orientations in subjects' pockets
  Mark significant rollercoaster events by peaks in some kind of flux
  Some kind of Dynamic-Time-Warping?

  Some kind of Expectation Maximization?
  Test accuracy by splitting rides into training/test sets

  ?

The ride has temporal variance (ex. sometimes getting from A to B takes 300ms longer)

The carts are driven by a motor on certain parts of the ride (e.g. up a slope). The acceleration should be fairly stable during those parts. Any temporal variation, if not attributed to different acquisition timings, is probably because of aerodynamic factors. For a given ride, that could be temperature of the air (which affects density), number of passengers in the carts (body parts disturbing the flow of air) and direction of wind. With an average speed of 82 mph (~132kmh), the roller coaster advances $\frac{132000m}{3600000ms} \approx 0.03666m/ms$. Therefore, a 300ms variation accounts for a "position error" of $300ms \times 0.03666m/ms = 11m$. Depending on the actual length of the ride, this might be negligible.
In addition to this, you have some jitter on your sampling frequency because the phone's OS doesn't guarantee that it will sample at an absolutely strict frequency.

Phones may have different initial orientations

It doesn't matter, if you subtract the first measurement from subsequent measurements then you have a relative indication rather than an absolute. For example, assuming that the first few measurements of your bearing were [270, 265, 260, 255, . . .] then the series [0, -5, -10, -15, . . .] indicates relative turning.

Phones stay in the subjects' pockets, but jostle a bit

Similarly, compared to the G forces that characterise the turning points in the ride, jostling of a phone in someone's pocket might be negligible. Having your phone tucked snugly in a front pocket is one thing, securing it in a zip pocket of a jacket that might flail about during the ride is another.

Learn the "mean" rollercoaster ride

as a time series
as a sequence of significant rollercoaster events

The "mean" roller coaster ride is the average of all available recordings for a given ride, per sensor. There is nothing to "learn" there, the more recordings you have, the closer your average approximates the "average" experience on the ride because you get a number of different accelleration / rotation / bearing samples for each time instance on the ride.
A sequence of "significant roller coaster events" is easy to extract once you have a clear definition of what a "significant roller coaster event" is.
For example, a "significant roller coaster event" might be a sharp drop, or a tight turn. Both of these are characterised by rapid rates of change in the total accelleration of a body on the cart. Therefore, you can derive a total accelleration from your 3D accelerometer as $A_{total}(n) = \sqrt{x(n)^2 + y(n)^2 +z(n)^2}$ and then apply some "event theshold $T_e$ on your $A_{total}(n)$. The indicator function that results from this (some $Z_n = aT_n \ge T_e$) points to some "significant roller coaster event". In the above, $x,y,z$ are the individual components of the accelerometer's measurements. Obviously, it is also possible to combine more than one sensors to extract a "significant roller coaster event".

On a new ride, in real-time, robustly determine where the subject is along the path so we can align music to their experience \m/
Bonus: 3D model the rollercoaster path

You will not be able to extract an accurate 3D model of the roller coaster path from accelerations because it is impossible to do accurate integration of the phone's sensors to a proper 6 DOF inertial nav system. If you want an accurate model of the roller coaster path, simply record GPS coordinates and you have the shape of the ride in 3D in zero time and with minimal data processing.
Tracking where a user is on a given ride is trivial once you have that "average track" signal generated. With that available, tracking a user on the ride is translated to "given a user's "past history" of a few points, where are they likely to be along this specific one dimensional curve?".
So, given that a "random" phone's signals are expected to vary in approximately a given way (the ride's shape), all you have to do is check where they are. You can do this by establishing a Kalman Filter on "position" which monitors current total accelleration and rate of turn. The Kalman will constantly be making a guess of where it should be along the curve and also be correcting this guess with the actual signals it receives from the sensors. You can achieve the same with a particle filter too, arguably with less processing.
Obviously, you can then set markers on the one dimensional curve to trigger your music replay.
Hope this helps.


What is the difference between ISI and MAI in MIMO system? ISI ： inter-symbol interference

MAI ： Multiple Access Interference

I know that multi-path would caused delay spread, and when signal symbol period is less than delay spreaad, ISI would happened.
In MIMO system, can I say that owning to the differenct location of antenna element that caused different path length, which make different delay and caused MAI ?
I also want to ask :
When would MAI happend?

MAI is when there is Interference due to multiple users accessing a common wireless media. So if you consider an OFDM system then as long as users are trasmiting in their allocated sub carriers there is no multi access interference irrespective of the MIMO technique used, delay spread etc.

If each user is using the same frequencies to transmit simultaneously, then there is Interference due to multiple access, this interference will happen irrespective of multipath, placement of antennas.

In case of Space division multiple access the it is theoretically possible to allocate precoding matrices to users such that they don't interfere with each others data streams while utilizing the same frequency resources. But this is only possible after suitable post processing at the receiver to separate out the data.

So anytime non orthogonal resources (time, frequency, codes, space) are used by users simultaneously there will be multi user/access interference, irrespective of multipath, delay spread, antenna placement.


What plotting / visualization tools are useful in generating good plots? It has been said a picture is worth a thousand words!
There have been a large number of softwares to write dsp codes and analyze signals and algorithms. It will be good to share on data plotting tools, that is - given a data set (1-d, 2-d, 3-d, ..higher dimensions (creative ways to visualize this)) what are the different tools available and advantages of each for plotting such data sets.
Advantages can be with respect to:
(a) quick/easy plotting
(b) formats of fugures (.eps, .svg) file which can be generated
(c) latex font embedding easiness
(d) post script figures
Tools like: matlab, octave, gnuplot, tikz, pgf etc are widely used by many of us. It would be useful to share your experiences (advantages/disadvantage) and best plot you have drwan using the tool you prefer!

I guess Python has great tools for this, especially IPython Notebook: http://ipython.org/notebook.html . It has good visualization opportunities and at the same time you could document your code and make it re-configurable and online. Of course it best for visualizing data low dimensions and creative visualization algorithms can be prepared either in OpenGL or high level languages such as MATLAB.

In 3D, if appropriate, time varying surfaces can be embedded in 2D and shown as a 3D volume. For generic 4D visualizations you could either use videos, or accumulated volume maps. For visualization of high dimensional data, always make sure you check methods such as Multi Dimensional Scaling.


How to do edge detection using Curvelet transform? I wanted to use Curvelet thresholding for edge detection, that is first take discrete curvelet transform then choose coefficient which are for curve and then reconstruct back to get an edge map.

Please help me on how to do Curvelet thresholding.

If you decide to stick to curvelets, you can find a precise description of an edge dectection algorithm is  Edge detection in microscopy images using curvelets, Tobias Gebäck and Petros Koumoutsakos, BMC Bioinformatics, 2009. The pdf is there. And you are very lucky: you can get, for free, an  implementation, CurveletUtils. It provides a Matlab source code for a GUI implementing the edge detection method. I hope you can get started with that.

If you care of other multiscale transforms providing different oriented  edges at different scales, you can take a look at a review paper  on 2D directional wavelets (including curvelets, shearlets, contourlets, etc.). It includes morphological and non-linear wavelets.

By the way, I do not understand why you call Canny and Sobel "morphological edge detectors".


How to normalize the power of a QAM signal? If I have, for example a 16-QAM signal, and the constellation points are (3,3),(3,1),(1,3),(1,1),...(-3,-3).

Now I have a simulated in coming signal with arbitrary power, how can I normalize the signal to this constellation so I can make correct decisions?

The first thing I could think of is to calculate the average power of the in coming signal and normalize it to the average energy of the constellation. But then when I have a very noisy signal, the noise power will effect the scaling factor considerably.

So is there some kind of least mean square method that I can implement so I can really get the scaling correct?

BTW I'm using Python/NumPy/SciPy for my simulation.

In $16$-QAM, four constellation points
$\left(\pm\frac{d}{2},\pm\frac{d}{2}\right)$ are "interior" points while
the remaining $12$ are "edge" points or "corner" points for which
at least one of the coordinates is $\pm\frac{3d}{2}$.
Thus, in
$n$ measurements $(X_i,Y_i)$ of successive transmitted symbols,
some of the $X_i$ will equal $\pm\frac{d}{2} + NI_i$ while
others will be $\pm\frac{3d}{2} + NI_i$.  Similarly,
some of the $Y_i$ will equal $\pm\frac{d}{2} + NQ_i$ while
others will be $\pm\frac{3d}{2} + NQ_i$.  Here, $NI_i$ and $NQ_i$
are independent zero-mean Gaussian random variables with variance
$\sigma^2$ representing the noise. Note that we can expect
about $\frac{3}{4}$-th of the
samples to have at least one large coordinate. Now, consider
that if we want to estimate the value of $d$ from the measurements,
and want to use the average energy calculations, then
$\sum_i (X_i^2+Y_i^2)$ is one possibility, but we are weighting
equally the symbols with large means and the symbols with small
means, and the noise component is a more significant part of
the symbols with small means than it is of symbols with large
means.  It might be better to use a weighted sum to estimate
$d$ (or $d^2$), and one extreme possibility is to give zero weight to
interior points ($(X_i$ and $Y_i$ are both small) and to use
only the large values of $X_i$ and $Y_i$ (edge and corner
points) in the estimation.  If the number of samples is large,
this could be as simple as taking those (approximately) $3n/4$
samples in which either $X_i$ or $Y_i$ (possibly both) is
large, and averaging the absolute values of the largest $X_i$ and
$Y_i$ to get an estimate of $3d/2$.

Once upon a time, when
analog circuits were in use, these kinds of calculations would
have been difficult to implement in hardware, and this effect
carries over into some present-day implementations that faithfully
follow the methods of the past that were originally devised
because of hardware limitations.
But these days, with digital computation being ubiquitous,
perhaps something new can be attempted.


Multi-input Single-output communication system In case of communication in MISO system, if we have a system of 4 x 1,  which means 4 transmitter with 1 receiver. Will the received signal be a vector of dimension 4 x 1 ? Or it will be 1 x 1 since the transmitter are sending to the same receiver?
Thank you

Multi-input single-output (MISO) is describing the inputs and outputs of the system. In communication, this is the channel, so if you have a single-output it means that the received vector is $1 \times 1$.
The general model for a MISO system is:
\begin{split}
y &= \mathbf{h^Tx}+\text{noise} \\
  &= (h_1x_1+h_2x_2+...+h_Nx_N)+\text{noise}
\end{split}
The single output in a MISO system is a combination of the multiple inputs plus the noise.


Determine minimum distance between two symbols in 16-QAM using the average symbol energy If the probability of bit error for a square M-ary QAM is

$P_M = (1-(1-P_\sqrt{M})^2)$

where

$P_\sqrt{M} = 2(1-\frac{1}{\sqrt{M}})Q(\sqrt{\frac{3E_s}{(M-1)N_0}})$

and $E_s$ is the average symbol energy, can I assume that $E_s=10A^2$?. $2A$ is the minimum distance between two adjacent symbols.

A symbol with coordinates $(x,y)$ has energy $x^2+y^2$. In 16-QAM, minimum distance $2A$ implies that the values of $x$ and $y$ are restricted to the set $\lbrace\pm A,\pm3A\rbrace$, and in consequence the possible symbol energies are $2A^2$, $10A^2$ and $18A^2$.

Furthermore, there are 4 symbols with energy $2A^2$, 8 symbols with energy $10A^2$, and 4 symbols with energy $18A^2$. Calculating the averagey symbol energy, we conclude that $E_s=160A^2/16=10A^2$.

Note that a similar procedure can be used to calculate the average symbol energy of any quadrature modulation, even if they're not square (or even rectangular).


Calculating phase response of filter with latency I have an LTI filter that I want to treat like a black box. It has a latency of 24 samples.

This is what I'm doing (which works for a filter with no latency):

Send unit impulse through my filter
Capture first 16,384 samples output
Perform FFT
For each of first 8192 complex values, take phase

Doing the above gives the red line shown below.

If I replace step 2 with:

Skip first 24 samples and then capture first 16,384 samples output

Then I get the blue line shown below.

The green line represents the actual phase of the signal.

Bode plot of phase vs log frequency (10Hz to 20kHz)

Thanks to everyone who provided helpful comments. As some suggested shifting the signal might be the way to go, but this results in the blue line where the green line is the phase I would expect.

It turns out to correct the phase after taking the DFT of the impulse response of the filter when the filter has a delay each bin needs the following angle adding to it:

$$
\frac{ 2\pi i L} {N}
$$

Where $L$ is the latency in samples, $N$ is the FFT size and $i$ is the $i$th FFT bin. So I was missing this as step 5.

For anyone else who stumbles upon this, the phase will end up being neither wrapped or unwrapped. So you will need to take the remainder of the final phase modulus $\pi$. This will give you the normal wrapped phase (as a correctly quadranted arctan would produce by just taking the phase of the complex valued FFT bin) and then you can unwrap as normal.


Complex filters used in analog domain Complex filters can be used in low IF architecture especially in DSP in order to phase out (reject) the image signal. How are complex filters implemented in the analog domain?

Vector Modulators and Single Sideband Mixer (frequency translator) are the common and relatively simple analog component for rejecting the image signal.

The top block diagram in this post is the architecture of a vector modulator and single sideband mixer and explains in further detail how it provides image rejection using complex signalling.

Frequency shifting of a quadrature mixed signal

Normally we see these components with a real input where an internal 90° splitter immediately converts the input to a complex signal and then it is further processed. However there is no reason these can't also be used in identical fashion with the quadrature splitter bypassed in the case of a low IF IQ signal already exists.

This would be the favored approach for eliminating the image signal.

However to implement a generic complex filter in the analog domain, you would need  to do the following with four real filters as shown in the diagram below. This is related to performing full complex multiplication with real signals as given by:

$$(s_I+js_Q)(c_I-jc_Q)$$
$$= (s_Ic_I+s_Qc_Q) + j(s_Ic_Q-s_Qc_I)$$

Corresponding to the multiplication of a signal (s) with a coefficient such as in an FIR filter (c) showing that we have an I and Q input, multiplied by and I an Q coefficient, resulting in a I (real part given by $s_Ic_I+s_Qc_Q$) and Q (Imaginary part given by s_Ic_Q-s_Qc_I).

Equivalently complex analog filtering (the filtering of a complex input with a complex impulse response) is formed with four real filters as shown in the graphic below. Each filter has a real impulse response given by $h_i(t)$ or $h_q(t)$ corresponding the real and imaginary portion of the complex impulse response desired.

Also to add, I did come across this: https://ieeexplore.ieee.org/abstract/document/4151912

Note that SAW filters are an analog equivalent of an FIR filter ("all=zero" filter, with poles at infinity, formed by the summation of weighted delays)


Relating bandwidth and wavelength of a carrier This came from a homework question, which I've solved. Just need some clarification.

I am asked to determine the bandwidth in a given amount of spectrum at a given wavelength. To find this, I know that;

$$\delta_f = \frac{c\cdot \delta_\lambda}{\lambda^2}$$

where $\delta_f$ is the bandwidth, $\delta_\lambda$ is the amount of spectrum, $\lambda$ is the given carrier wavelength, and $c$ the speed of light.

That's all well and good, but I can not find any reason for why this formula works, or where it comes from. Where does this formula come from, or what is it called, or how is it derived?

The question reads as follows;

  How much bandwidth is there in $0.1 \textrm{ micron}$ of spectrum at a wavelength of $1 \textrm{ micron}$?

I will try to derive your formula.  The formula you were given is technically not correct, but on a practical level it is.  I will try to explain below.

The formula for the frequency of any wave (sound, water, light, etc.) is $f = \frac{\nu}{\lambda}$, where $\lambda$ is the wavelength and $\nu$ is the wave's velocity.  This makes intuitive sense if you think about it because you will see more peaks the faster the wave travels, and you will see fewer peaks the longer the wavelength is.  In your case, of course, the wave velocity is $c$.

Bandwidth is the highest frequency minus the lowest frequency.

$$
\begin{align}
B &= f_h - f_l \\
&= \frac{c}{\lambda_l} - \frac{c}{\lambda_h} \\
&= \frac{c}{\lambda - \frac{\Delta\lambda}{2}} - \frac{c}{\lambda + \frac{\Delta\lambda}{2}}
\end{align}
$$
At this point to make the typing a little easier I will call $\frac{\Delta\lambda}{2}$ "$d$".
$$
\begin{align}
B &= \frac{c}{\lambda - d} - \frac{c}{\lambda + d} \\
&= \frac{c(\lambda + d) - c(\lambda - d)}{(\lambda - d)(\lambda + d)} \\
&= \frac{2cd}{\lambda^2 - d^2} \\
&= \frac{c\Delta\lambda}{\lambda^2 - \frac{\Delta\lambda^2}{4}}
\end{align}
$$
The formula you were given is, as a practical matter, correct, because usually $\Delta\lambda$ is much smaller than $\lambda$.  In your case $\Delta\lambda$ is 1/10 of $\lambda$, which is actually much bigger than any normal situation.  Usually $\Delta\lambda$ is orders of magnitude smaller than $\lambda$.  Even at 1/10, though, $\frac{\Delta\lambda^2}{4}$ is $\frac{\lambda^2}{400}$, which is negligible.  Thus, it is dropped.

I hope this helps.


What are the advantages of wavelet transform for image multiresolution analysis? Compared to gaussian pyramid and laplacian pyramid, what are the advantages and disadvantages of wavelet decomposition for multiresolution image analysis? Thank you.

Given the concise aspect of your question at the present time, a short answer is: wavelets are able to provide a critically-decimated representation (or non-redundant) of images, even orthogonal or close to. This is beneficial for applications such as compression (criticality), and for theoretical derivations (orthogonality). This can be a drawback for other practical uses.


Is this code for doing a moving average correct? Is this code for doing a moving average correct?

y=conv(x, ones(1,100)/100);

I want to filter my speech signal with moving average filter in code.  I think that I can use convolution for this.

Yes, this code is correct for implementing a moving average filter.

Nevertheless I do recommend to use the in-built smooth function in MATLAB.

y = smooth(x, 100, 'moving')


Why Do Most of The Papers Use the Frobenius Norm for Denoising? I have an noisy image and I want to remove noise from it; suppose $y$ is noisy image and $A$ is linear mask which makes my image noisy and $x$ is original image, so we have
$$
Ax + \eta = y
$$
and $\eta$ is additive noise. I have read in many papers that the equation below is a good cost function for removing noise:
$$
x^* = \arg\min_{x} \lambda\|A(x) - y \|^2_\mathcal{F} + \|x\|_{TV}
$$
First term aims at minimizing additive noise but why do most of the papers use the Frobenius norm? And second term wants preserve edges and features of images.

The Frobenius Norm has multiple equivalent definitions – the useful for error measure is probably this one:

$$\left\|M\right\|_\mathrm F = \sqrt{\sum_{p\in M}\left\lvert p\right\rvert^2}$$

That's a root square over all pixels.

Root mean squares are very useful cost functions, as they describe the power of a signal.


How to find correlation between 3 signals I have 3 signals, S1 S2 S3, with a same time domain and same time steps. I'm wondering if there is any correlation between these 3 signals. more specifically I want to know if it is possible to make Signal S3 by any combination of signals S1 and S2.

You need to process them jointly in pairs to get the correlation coefficients
$$\rho_{13}=\frac{\mathsf{Cov}(S_1,S_3)}{\sqrt{\mathsf{Var}(S_1)\mathsf{Var}(S_3)}}$$
$$\rho_{23}=\frac{\mathsf{Cov}(S_2,S_3)}{\sqrt{\mathsf{Var}(S_2)\mathsf{Var}(S_3)}}$$ that represent the linear dependencies between $S_1$, $S_3$ and $S_2$, $S_3$, respectively. So when $|\rho_{xy}|$ is close to one, it indicates high linear dependency (and values close to zero mean small linear dependency).

That said, the problem you mentioned can be solved using a least squares approximation. In MATLAB, assuming your signals are in form of column vectors S1, S2, and S3, you can do this to find the weights:

w = [S1 S2]\S3;

So the signal S3 can be approximated as

S3_hat = [S1 S2]*w;

and the mean squared approximation error is

mse_err = mse(S3-S3_hat);


Why RHP zero phase is not 180° to 90° The asymptotic phase behavior of an RHP zero is from 0 degrees to -90°, the mirror of an LHP zero. Graphically, I'm confused about why this is the case and the phase is not from +180° to +90°. See the below image for my reasoning/confusion. Is the angle definition not consistent?

You're right, a single (real-valued) zero in the RHP could cause a phase that moves from $\pi$ at $\omega=0$ to $\pi/2$ for $\omega\to\infty$. However, if you rewrite the corresponding contribution to the transfer function as $s_0-s$ (i.e., you change the sign), then you obtain a phase shift between $0$ at $\omega=0$ and $-\pi/2$ for $\omega\to\infty$.
The difference is just a sign inversion, i.e., an addition of $\pm\pi$ to the phase. The contribution of a single real-valued RHP zero in the form $s_0-s$ does not invert the sign for $s=0$, whereas if you write it as $s-s_0$ you get a sign inversion, hence the additional phase shift of $\pi$.


python rewrite lfilter (iir) with for It's a beginner question, but useful to users from python - signal.lfilter, I was using lfilter from
Find reverse one pole lowpass filter

The doubt now is, how to rewrite the iir filter to 'revert' the low pass filter, without using lfilter function:

low pass:

y[i] = y[i-1] * 0.8702 + x[i] * 0.1298
y[i] = signal.lfilter(b,a,x)

low pass revert:

y[i] = y[i-1] * ?? + x[i] * ?? + ??
y[i] = signal.lfilter(a,b,x)