processing gain in the receiver I am always confused about the processing gain in the receiver. Assume my receiver had the following:
$$\boxed{\text{LNA}}{\longrightarrow}\boxed{\text{Mixer}}{\longrightarrow}\boxed{\text{Filter}}{\longrightarrow}\boxed{\text{ADC}}$$

Does the filter in the analog section has any role to play with the processing gain?
I know if the ADC is 12 bit and 100 msps sampling its SNR is $$12\cdot6.02+1.76=74\rm \ dB$$
if I sample at 200 msps does processing gain increase?
If I am using FFT analysis the processing gain increases, but if there is no FFT how to increase processing gain in digital domain (by filtering)?
What is the significance of processing gain?

"processing gain" isn't something that a specific step in your signal processing brings you, for every possible signal, in every possible situation.

Instead: If you have knowledge about properties of your signal, you can use that to somehow discriminate signal from noise and thus increase the SNR. The increase in SNR is called a processing gain.

Notice how you didn't mention any properties? Especially none you could exploit for this?

Typical properties include narrow bandwidth (where a simple filter could reduce the amount of noise power and thus increase SNR), or being shaped with a specific filter to which you can use a matched filter (where you specifically design a filter so that it maximizes SNR in AWN), or you could know your signal is spread with a specific spreading sequence (where you can despread and thus exploit that you have power growing quadratically with the number of accumulated correlated signal samples, but only linearly with the number of uncorrelated noise samples).

But really, these aren't all properties one can exploit. There's things like known temporal behaviour, channel state information, and many nonlinear noise operations and much more that can have a processing gain for any specific system- and noise-model. So, you need to write that down and see how your processing step affects the SNR.

  if i am using fft analysis the processing gain increases

No! Or, only for signals that match the implicit assumption you're making (which I don't know, but which seems to assume your signal is very narrowband and fits exactly into one of your FFT bins; don't know how often that is the case in an unsynchronized system).

  but if there is no fft how to increase processing gain in digital domain

That depends on your signal model and your noise model, and we can't tell you that.


FIR filters: direct form / transposed FIR I am learning about FIR filters and I'm confused.
I am trying to find out different types of FIR filters.

Is direct form and n-tap FIR filter the same?
What does transposed FIR filter do?

A finite impulse response (FIR) digital filter implements the following convolution sum

$$y(n)=\sum_{n=0}^{N-1}h(k)x(n-k)\tag{1}$$

for each output sample $y(n)$, where $x(n)$ is the discrete-time input signal, $h(n)$ is the filter's impulse response, and $N$ is the filter length. The values $h(n)$ are also called filter taps, and $N$ is then referred to as the number of taps. The filter described by Equation (1) is also called an $N$-tap filter.

Direct-form and transposed direct-form are just different implementations, i.e. different ways to compute the sum in (1). In theory they are identical, but when computed with finite precision, there can be differences between the different implementations. The direct-form FIR structure is also called tapped delay line or transversal filter.

The two realizations below are the direct-form structure (transversal filter, tapped delay-line) and the transposed structure (from Oppenheim and Schafer, Discrete-time Signal Processing):


How to demodulation BPSK in real case(through a channel)? I am configuring a channel sounding system. Using PRBS to correlate for time delay. Some details as follow.

Generate PRBS.
Using BPSK to modulate PRBS and Tx transmit.
Rx receive signal and demodulation.
Correlation the signal that demodulation from BPSK and original PRBS.

But now, I am stock at step 3. The following block diagram is the system that I set and the signal representation of each node.

I would get the wrong demodulated results that if (delta1) minus (k1R1 or k1R2 or k1R3) >= 90deg .
How to demodulation BPSK in real case with multipath??

There will inevitably be frequency offsets and time offsets in the OP’s receiver, given the receiver is working on a completely independent clock, so cannot have the exact same carrier frequency (and in the case of mobile stations there is also Doppler offset). The typical BPSK receiver will use carrier and timing recovery loops to extract the estimates of the carrier and time offset from the signal itself.
A very simple carrier recovery for BPSK is done by simply squaring the received signal, which will produce a dominant synchronized tone at twice the carrier frequency. On its own this would usually be too noisy to use as the carrier, so is typically locked to a clean local oscillator using a PLL and then dividing the resulting frequency by two, which is then multiplied by the received signal to demodulate the BPSK signal.  If one has the flexibility of using a carrier that is an integer multiple of tour  symbol clock, then timing recovery can be done simply with additional frequency dividers; otherwise the Gardner Timing Loop described in the links below can be considered (an error discriminator is needed such as the Gardner Timing Error Detector; this error then gets accumulated and the accumulated error weighted by the desired loop gain controls the sampling location for the timing clock resulting in true recovered data).
The recovered carrier will have 0/180° phase ambiguity so further information to resolve this is encoded in the data itself. Data can be encoded with a known preamble that has many data transitions to allow for resolving this ambiguity and also to answer the OP’a multipath question, help train a channel equalizer. An example preamble could be a 11 bit barker code (10110111000) followed by (1010101010), or depending on acquisition SNR required could be a longer series of these or similar codes. The barker code has the excellent autocorrelation property that it is not highly correlated with shifted versions of itself so can help establish the correct start of a data packet. This autocorrelation property (ideally an impulse at $\tau$ = 0 ) is also ideal for training an equalizer since it is approximately white in frequency: An equalizer can only compensate for the frequencies it has measurements for. The 101010...sequence is very friendly to timing recovery given the transitions on every bit.
In addition the signal level will likely be variable so one would likely want to have an AGC (automatic gain control) loop as well.
Other posts where I provide more details on carrier and timing recovery and equalization all applicable to BPSK are summarized below:
Carrier Recovery:
Phase synchronization in BPSK
FFT-based coarse carrier recovery for QPSK
High modulation index PSK - carrier recovery
Recovering signal for psk
Timing Recovery:
High modulation index PSK - carrier recovery
Location of Matched Filter
Isn't Gardner's algorithm and Early-Late gate the same thing?
Channel Equalization:
Compensating Loudspeaker frequency response in an audio signal
non linear equalizer vs linear equalizer
Does zero-forcing equalizer need known channel impulse response?
The benefits of a fractionally spaced equalizer
Recursive Least Square Adaptive Linear Equalizer


Noise amplification with inverse filtering I'm trying to gain a better (mathematical) understanding of why inverse filtering is almost never the solution for correcting an image. From what I understand, we start with a discrete signal $s[n]$ that is the result of the true image $y[n]$ convolved by some kernel $A$ with some added noise $w[n]$. For discrete signals, this be converted to a matrix equation $s = Ay + w$, where $A$ is a circulant matrix representing the convolution. In Fourier Space, we have $s[f] = A[f]y[f] + w[f]$ which is much simpler.

Now, if we assume we know the kernel and its inverse exists, then we can approximately recover the true signal $y$ by inversion : $A^{-1} s = \hat s = y + A^{-1} w$. The problem becomes clear in the Fourier space, where inversion then becomes division, and if $A[f]$ is $0$ or very small at several frequencies this can lead to noise amplification. I am trying to get a better picture of this in the matrix representation.

I want to see if I can measure the noise amplification in this representation by trying to figure out an upper bound on the norm of $\hat s$ but I'm not sure if I'm going about it the right way. Calculating the norm of $\hat s$ I have $$|| \hat s|| = ||y + A^{-1} w||$$ $$||\hat s|| \leq ||y|| + ||A^{-1}||~||w||$$ Am I correct in assuming that the upper bound on the noise amplification is the norm of the inverse of the kernel?

The noise amplification would be the ratio of the norms:

$\frac{\text{Norm of noise after inverse filter}}{\text{Norm of noise before inverse filter}} \leq \frac{||\mathbf{A}^{-1}||*||\mathbf{w}||}{||\mathbf{w}||}=||\mathbf{A}^{-1}||$

It seems your reasoning is correct that the norm of the inverse filter is an upper bound on the amplification.


minimum sampling rate for very short duration signals I understand that sampling rate should be twice the signal bandwidth for successful reconstruction. But what I don't get is how this translates in the time domain.

For example, I have a sine wave pulse with frequency 100MHz and duration of say 50ns. This signal has a bandwidth of 20MHz (1/50ns = 20MHz). Theoretically, I can sample this signal at any frequency above 40Mhz. Let's say I sample at 50MHz, then I only get 2 samples per pulse!! Now what can you possibly do with those 2 samples is my question. In other words, although Nyquist criterion is satisfied in this example, there's not enough data points to proceed.

This kind of signals are common in radar and generally you want to run some sort of matched filtering/correlation once you have digitized the pulse. But you can't run any descent signal processing on just two samples? am I missing something here?

Any signal of finite extent (duration) in one domain (time) has infinite extent (bandwidth) in the other (frequency).  The bandwidth of a short pulse above some noise floor can be extremely high, due to the pulse window's transform convolved with the carrier frequency.

Sharp edges or cutoffs imply very high bandwidths.  Thus in practice, most designed pulses are shaped closer to Gaussian, with much wider tails than the main pulse lobe (in both domains).

Even so, if you know the exact frequency of an unmodulated fragment of a pure sine wave, then in zero noise, 2 or 3 non-aliased points are theoretically enough to solve the equations for the sinewave's  amplitude and phase.


Would it be possible to reconstruct 3d trajectory of ball projectile using only monocular video Assume ball location on the image frame in each timestamp is given.

Would it be possible to reconstruct the ball's 3d trajectory in the camera's coordinate frame?

The ball should be able to described using constant acceleration model, as it only subject to gravity and initial throwing force

Edit: I am hoping if someone could points out whether it be possible to deduce the unobservable perpendicular motion, after some flying time

I guess this is what Olli implies: Using a lens distortion matrix if necessary, the center of the ball in the picture gives you an elevation and azimuth of your ball relative to the focal point of your camera. The size of the ball image gives you distance.

That's a point in 3D space.

From there, it's really only a manner of grammar school level math to estimate the four (or only three – do you know the mass of the planet you're on?) unknowns in your problem (three components of initial velocity, and maybe magnitude of gravity), at least without considerations for measurement error from multiple observations of the ball.

If you need to take measurement error (finite shutter speed, finite frame rate, finite camera resolution, imperfect estimation of position from picture) into account, well, feels like the classical problem where you'd throw a Kalman at it and be done.

In fact, I think this nice little article "How a Kalman filter works, in pictures" has an explanation that nearly completely describes a solution for your specific problem.


Average Energy of modified QAM Constellation How to compute the average energy for the modified 16-QAM constellation below? I know that for regular/rectangular 16-QAM we can utilize the 4-PAM along each axes but I'm not sure how can we find $E_{avg}$ here.

The procedure is always the same. You need to compute the expectation $E\{|A_k|^2\}$, where $A_k$ are the complex symbols of the constellation:

$$E\{|A_k|^2\}=\sum_kP_k|A_k|^2\tag{1}$$

$P_k$ is the probability that the $k^{th}$ symbol is chosen. Usually you can assume that all symbols are equally likely, i.e., $P_k=1/M$, where $M$ is the number of symbols.


Sample rate changing audio pitch I am using rtl-sdr at an adjustable sample rate to get stuff from a WBFM radio station.  I then filter and demodulate it.  When I sample at 2.4Msps, the sound is 'choppy'.  I tied the sample rate to a slider on gnuradio, so I can vary the software radio sample rate by moving around the slider.  It turns out the pitch of the person talking/singing's voice gets lower as I move the sample rate down, and higher as I move the sample rate up.

Why?  2.4 million samples / second is much greater than 2*20,000Hz.  Why would there be a difference between sampling at 2.4 and 2.8 Msps?  All the frequency content is already covered in both cases...

Also: in gnuradio I can adjust the width of the low pass filter that is placed before the FM demodulator.  I thought that making the filter very narrow would remove the high frequency content of the audio.  Well, it doesn't.  It just makes the sound 'staticy'.  Why?

A couple of things jump out when looking at your diagram:

The first low-pass filter has a decimation of 4, so the sampling rate at its output is 256 kHz. However, you have set the sampling rate of the FM receiver and the second LPF to 224 kHz. These should be 256 kHz.
The cut-off frequency of the first LPF is 150 kHz. This means a total passband of 300 kHz. This is larger than necessary, so you'll end up picking up more noise than the optimum.
You multiply the audio output by a constant factor of 2. This may be too little, or too much. I'd suggest using a slider to control the volume. If you connect a scope sink to the audio sink input, you'll be able to easily get the optimum gain: you want the signal to be withing [1,-1].


Wavelets / Identifying bursty signals I thought my network was experiencing outages. I collected bandwidth (Mbps) measurements every 30 min.

A cumulative probability distribution shows the probability that bandwidth will exceed any given value (e.g. P(x > X)). However, this was not very insightful.

A continuous wavelet transform using PyWavelets did not show any strong peaks (i.e. heat map of cwt in scale and time). That meant to me that there was not a consistent duration for the outages.

What would be some ways to analyze the signal to see the time scale of outages?

Sample Data

Cumulative Probability

CWT using PyWavelet

I thought my network was experiencing outages. I collected bandwidth (Mbps) measurements every 30 min.
. . .
What would be some ways to analyze the signal to see the time scale of outages?

An "outage" is a random event. Therefore, looking for regularity in outages might be hinting at something more widely already known.
The answers to the questions in the comments are very much dependent to the objectives. But I will try to make this a bit more constructive.
For example:
The CDF presented is over 1400 instances at 30 minute intervals and it therefore gives you the probability of sustained download at a particular speed over 1400*30 minutes or about 30 days (If the horizontal axis depicts the index of the measurement). If you were asking if your line would sustain a download speed over the next 15 minutes, the answer would be "I don't know", because you have not measured anything more frequently than 15 minutes.
On their own, wavelets and discrete fourier analysis cannot tell you much. They can tell you that a particular segment of the signal suddenly requires many more components to describe it and this usually happens at discontinuities but this helps in recognising a "spike" which you can already do here by thresholding (for example).
Rather, if you suspect that there is some regularity in the outages and you want to estimate it, then what you can do is some form of time series decomposition.
You can run autocorrelation for example and this will already give you a handle on periodicity. Or a full blown discrete fourier transform (or even wavelet transform) to determine periodic components more accurately.
By determining the Trend and Seasonality components of your signal, you can then create "filters" to discover areas where it veers off its usual course. (And then of course, having discovered the timestamps of those events, you can analyse (or try to predict) the time that passess between them and see if there is anything useful in there).
But notice here, we always need a time reference. For example, when you assess Trend, a trend over what time scale? This will then determine how often you need to observe the signal.
Hope this helps.


Stitching together audio buffers with potential discontinuity at the boundary I am writing a sound library for my c++ hobby game engine. The sound library mixes music and sound effects into a single buffer that is sent to the audio hardware.

While mixing, any number of source sounds may appear or disappear between two output buffers causing a discontinuity at the border between the buffers.

The discontinuity can be heard as a disturbing high frequency click when the mix is played by the audio hardware.

I have solved this by doing internal bookkeeping of what sound sources are playing and neatly fading sound sources in and out. However this is more involved than it would first seem. The mixing code becomes difficult to read and maintain.

So I would prefer to instead simply handle any discontinuities in the final mix buffer before sending it to the sound card.

I have experimented with a wavelet transform solution which takes an area around the border between two buffers and adjusts any excessively large coefficients until they fit in with the overall distribution.

However that technique seems tricky and introduces a latency to the output because potentially quite old samples have to be updated if the discontinuity carries all the way down to low frequencies...

It would be preferable if the solution would only alter the new buffer and potentially carry over any residual errors at the end of the new buffer to be resolved in the next buffer.

Are there any known techniques for this?

Any help is greatly appreciated!

In 1998 I made a "zero level adjustment click removal" system that shifted the zero level to fix the discontinuity, and this "zero level" was gradually moved back to 0, by decreasing it by 1/2048 of its value every sample. It worked well for low-frequency sounds, but when a high-frequency sound ended abruptly at a large sample value, that created a low-frequency click. I ended up removing the system and simply adding more channels in order to fade out sounds. For this I used a one-pole lowpass filter on all volume envelopes, and when a sound was killed the volume envelope was set to zero and some time was allowed for the fadeout before making the channel available to other sounds. If too many channels were reserved for fading out, the one with the smallest amplitude value was freed to use.

But the thought kept bugging me that if the zero level adjustment would be done only to a lowpass filtered portion of the output, then it might sound better. Purely high-frequency sounds don't really need click removal.


Plotting the discrete time signal Consider the following analog signal

$$x(t) = 2 \sin(100\pi t)$$

The signal $x(t)$ is sampled  with a sampling rate $F_s=50\textrm{Hz}$. Determine the discrete time signal. Plot the discrete time signal.

Also determine the total number of samples.

I don't understand how to approach this question.

$x(t)$ is a model for an analog signal. After sampling, the real $t$ variable is turned into an integer variable $k$, corresponding to a regular splitting of the real time interval, starting from some origin time $t_0$, evenly spaced by a sampling period $T_s$ being the inverse of the sampling frequency: $T_s = \tfrac{1}{F_s}$.

A discrete signal is thus:

$$ x[k] = x(t_0+k.T_s) = 2\sin(100\pi (t_0+k.T_s)) = 2\sin(100\pi t_0+k.100\pi/50)) $$

Thus, since $$2\sin(100\pi t_0+k.100\pi/50)) = 2\sin(100\pi t_0+2k\pi))= 2\sin(100\pi t_0))$$

the discretized signal is determined by a constant value, depending only on the initial sampling time $t_0$. You can plot if with many plotting software (even Excel).

The total number of samples is ambiguous. Since $k$ take any integer value, the number of discrete samples is infinite. Since it is constant, $x[n]$ takes only one value.

Explanation: you start from a continuous sine, and only get a constant signal. This is an illustration of the Nyquist-Shannon sampling theorem (or aliasing). If the sampling frequency is not consistent with signal maximal frequency, the discrete signal might be a poor version of the analog signal.


Negative frequency in the Fourier Transform Why use this instead of a positive frequency? I asked on the math forum but it just wasn't helpful.

I understand the idea of a negative frequency is important in general since many real signals like cosine can be constructed using complex exponentials of two opposite rotating frequencies. But why in the FT? Also would it be correct to view this as 'frequency' since in the complex domain it can be negative? I recall someone arguing it is not frequency but index: https://www.quora.com/What-is-the-meaning-of-negative-frequencies-in-the-Fourier-transform

An interesting explanation I had never seen was from wikipedia:

But can any signal input be modeled as a complex exponential form as above? For example, if your input was a rectangle function.

I think you may be confused about the use of negative frequencies and the sign of the integration coefficient.

Let's take a look at the inverse Fourier Transform

$$ x(t) = \int_{-\infty}^{+\infty}F(\omega)\cdot e^{ 2 \pi i  \omega} d \omega $$

That basically means that you can construct any time domain signal as the from a set of complex exponential. So you can break down a complicated signal into a number of constituent parts that are hopefully easier to work with.

That's typically called an orthogonal base. Complex exponentials are not the only choice, you could use a orthogonal polynomials, wavelets, etc. Complex exponentials are popular since they are easy to work with, specifically when you are dealing with linear time invariant systems

Note two things

The integration goes from $-\infty$ to $+\infty$, i.e you need both negative and positive frequencies to generate any possible time domain signal
The integral variable is positive, but that does only mean than you first integrate the negative frequencies and then the positive ones.

Of course you need to know what the weights for the complex exponentials, so you need to invert the formula above. You get

$$F(\omega) = \int_{-\infty}^{+\infty}x(t)\cdot e^{ -2 \pi i  \omega} d t$$

That's just how the math works out. The integration variable is negative, but that doesn't mean anything. The formula needs to produce results for both positive and negative values of \omega since we need both positive and negative frequencies to represent a time domain signal.


Zone Plate Reconstruction I am trying to reconstruct the Zone Plate image and am struggling to remove the last remaining aliasing.

ShaderToy: https://www.shadertoy.com/view/wdGGWK

In the shader above you can see that I am sampling each pixel independently 1000 times with [-3, 3] jitter, and filtering those samples with (currently) Lanczos2. Why is there still aliasing?

My understanding is that the aliasing just at the edge of the visible pattern is from the filter and can't be removed without a different choice of filter, but the subtle aliasing in the middle of the axis seems like an error, and one that I have been unable to remove.

Any ideas?

Lanczos2 has a frequency response peak of about -40 dB or 0.01 at about 0.9 times the sampling frequency, see this question. Your peak aliasing outside the main lobe is just left from where the aliasing goes to zero frequency so I think it is the same peak. The amplitude of the aliasing is also roughly as expected.

So I don't see any error by you.


What measurement errors occur in an inteferogram? This question is bordering on physics, but it is the signal side that I am interested in so I hope this is the correct forum.

I am very familiar with instrumental noise sources in CCD based spectrographs and what their signal properties are. I have come to realise how important understanding noise sources are in order to effectively mitigate for them.

Now I am taking an interest in FTIR which uses an interferometer to do the measurement and would like to understand how measurement errors manifest in the interferogram and how they manifest after the FT.

My suspicion is that the variations due to the technical act of measurement will share some aspects with CCD instruments but also present some sources that are very different. For example, spectrograph-CCD systems are affected by e.g. pixel-to-pixel shot noise, fixed pattern response, read out noise, spectra of optics (blank signal), dark signal, vibrations, alignment etc.

How do the following present in both the interferogram and the FT?

random noise (i'm guessing its very different to shot noise in CCDs, not only pixel-to-pixel)
vibrations (physical or thermal)
mechanical whiplash

As I am not familiar with the technological errors associated with an interferometer I'm sure I've not anticipated some sources, so please educate me  on any other sources of variation due to the technology (not the sample or sample-technology interface)

I've tried some searches but not found any useful resources so any good resources would be appreciated

You did an excellent job in presenting the background and context of your question and your listing of error sources is fine. But FTIR also presents some other issues, particularly given the relatively poor SNR performance of IR detectors compared to UV-Vis regime detectors. So two things. First, I recommend this book: P.R. Griffiths, J.A. de Haseth, Fourier Transform Infrared Spectrometry, 2nd Ed., Wiley-Interscience, John Wiley & Sons, Inc., NY, ©2007.

Second, when I simulate FTIR using my optical calculus software, I generally include light source temperature fluctuations, jitter on the moving mirror and, of course, detector and preamplifier noises. The next few figures show the simulation model (with the 3 noise souces I stated), which looks somewhat like a block diagram, but is actually the simulation program, and some results for ethyl acrylate in the 400 to 4000 wavenumber range. Figure 1 follows:

Above is Fig. 1, the FTIR simulation program with ethyl acrylate as the test substance. The simulation conditions are shown in Fig. 1.

Figure 2 is the last of the ten 16k interferograms generated. It takes about 2.5 minutes per interferogram on my 2006 iMac 24. (It would be much faster on a modern Windows-based PC, but speed isn't everything and I am currently holding off buying a new PC.) Figure 2 follows:

The next figure (Fig. 3) is an expanded view of the center of the interferogram:

Now Fig. 4 shows the pseudo-transmittance of the ethyl acrylate test substance:

Then Fig. 5 shows the pseudo-transmittance of the background:

Finally, dividing the two pseudo-transmittances and taking the negative log yields the absorbance spectrum of ethyl acrylate, as shown in Fig. 6:

My software runs on PCs (and my real old Macs) and is free, with full commented source code. However, it is an add-on to a commercial simulation program called ExtendSim® (with which I am not financially associated: I am just a very early user).

If I can find any relevant and detailed papers discussing FTIR considerations in particular, I will edit this answer to include them. And it would benefit me as well, even though the spectroscopy I used to do was not in the IR.


How to derive the answer to this convolution problem? I came across this below question, (which was a homework assignment question for Signal Processing class, which my friend mailed me for help solving), mulled over it for an hour and had no idea how to proceed with solving it.

Let $C(x) = A(x)B(x)$ where:

$$A(x)=\sum_{n=0}^{N_1}a(n)x^n$$
$$B(x)=\sum_{n=0}^{N_2}b(n)x^{2n}$$
$$C(x)=\sum_{n=0}^{N_3}c(n)x^n$$

Find expressions for $N_3$ and $c(n)$ as functions of $N_1$,$N_2$, the $a(n)$ and $b(n)$.

Apparently, it has something to do with convolution, as in $c(n)$ is the convolution result of $a(n)$ & $b(n)$ or something like that.
But I still can't figure out how that is. Can anybody please explain me the answer of $N_3$ and $c(n)$

HINT: what is the highest power of $x$ after multiplying $A(x)$ and $B(x)$? This gives you directly the value of $N_3$. Then rewrite $B(x)$ as

$$B(x)=\sum_{n=0}^{2N_2}\hat{b}_nx^n$$

and you can use normal convolution of $a(n)$ and $\hat{b}(n)$ to derive $c(n)$. Now you just need to express $\hat{b}(n)$ in terms of $b(n)$.

EDIT: OK, so here's the solution:

We have

$$C(x)=\sum_{n=0}^{N_3}c_nx^n$$

Since $C(x)=A(x)B(x)$, the highest power of $C(x)$ must be $N_3=N_1+2N_2$. You can rewrite $B(x)$ as

$$B(x)=\sum_{n=0}^{2N_2}\hat{b}_nx^n$$

with

$$\hat{b}_n=\begin{cases}b_{n/2},&n \text{ even}\\
0,&n \text{ odd}\end{cases}$$

Now the coefficients $c_n$ can be written as the convolution of $a_n$ and  $\hat{b}_n$:

$$c_n=\sum_k\hat{b}_ka_{n-k}$$

For the index $k$ in the above sum we have the following constraints:

$$0\le k\le 2N_2\quad\text{and}\quad 0\le n-k\le N_1$$

which results in the summation limits

$$c_n=\sum_{k=\max\{0,n-N_1\}}^{\min\{n,2N_2\}}\hat{b}_ka_{n-k},\quad 0\le n\le N_3$$

If you sum only over even $k$, you can replace $\hat{b}_k$ by the coefficients $b_{k/2}$:

$$c_n=\sum_{k=\max\{0,n-N_1\},k\text{ even}}^{\min\{n,2N_2\}}b_{k/2}a_{n-k},\quad 0\le n\le N_3$$

which can be rewritten once more as

$$c_n=\sum_{k=\lceil{\max\{0,n-N_1\}/2}\rceil}^{\lfloor{\min\{n,2N_2\}/2}\rfloor}b_{k}a_{n-2k},\quad 0\le n\le N_3$$


Typical number of samples per symbol Assuming a pulse shaping of a source of symbols (like $s=\{3,-i,2+5i\}$), and a SRRC waveform (or any other waveform). I know that it depends on many factors like modulation scheme, but, in practice, are there any typical value or range of "number of samples per symbol" required for the pulse shaping process?

2 to 4 samples per symbol is common in my experience. The incentive for using fewer samples per symbol is an overall lower sample rate (and thus less processing load). It depends what point in the system you're talking about, but you also need to take the excess bandwidth of your pulse shape into account to ensure that the sample rate is high enough to avoid aliasing.

In some parts of the processing chain, you can get away with fewer samples per symbol than you might expect given the Nyquist rate and the signal's bandwidth. If the signal is perfectly synchronized, for example, then you don't need additional oversampling. In that case, you can use just one sample per symbol (basically, the symbol values themselves). You can then use that and the known pulse shape to reconstruct the signal to whatever fidelity you desire.

In the absence of perfect synchronization (for example, in the various receiver blocks that perform frequency, phase, and/or timing synchronization), you'll have some amount of oversampling. In the literature you'll find 2 samples per symbol as a desirable target rate, because it's really the smallest that you're realistically going to achieve. You can always go higher if you need to. It just results in a higher sample rate (so more computational work)


So is there such a thing as a DSP engineer nowadays? I just read the article linked to below:

http://www.eetimes.com/document.asp?doc_id=1275520

Now, I'm planning on doing a Masters in Signal Processing. I'd like to know what the career prospects are for someone graduating in this specialization. Note that I live in the high-tech capital of India, so I should be in about the best place in this part of the world for any tech related positions.

I really like the idea of applying mathematics to solve practical problems, so that's why I picked DSP. Am I too late to the party, as this article seems to imply?

Yes there is. There definetely is there such a thing called DSP  Engineer, but it is not exactly what the article described, i think it even goes more high level.

Nowadays, most companies relies on already established Knowledge such as standards, or papers from IEEE or AES, or such places. So a DSP engineer requires now more than ever, a lot of knowledge about Programming, in a low level fashion, programming in Assemly, C/C++ to producte robust, reliable, mantainable, understandable and expandable code. It is true what the artcile says that the are cookbooks, so if you don't remember exactly how LMS algorithm goes, you can pick up the book and refresh it in 5 minutes. But the concept needs to be there. If you are completely new to DSP, but you know a lot about programming and Kernel and operating systems, you will have a hard time implementing algorithms if you don't know the basics.

There is no real answer to your question but i state again my first comment. Yes there is, but it has incorporated a lot more about embedded Systems that what it had bewteen 1960 and 1980.


"Low-pass filter" in non-EE, software API contexts I am an experienced software engineer and am working on smartphone sensors. I've taken fundamental EE classes in DSP and am trying to apply my knowledge. I believe that I understand convolution, transfer functions, z-transform, etc. I know a little bit about FIR and IIR filters.

Now, when reading through software APIs and documentation, I see people are applying a LPF to sensor data in the time domain. I know that you do that through the use of difference equations (e.g. y[i] = y[i-1] + 2*x[i]), but I learned in my EE class that LPF are typically applied through the convolution operation where you convolve the time signal with the coefficients of a sinc wave (for example) and with a specific cut-off frequency. So the colloquial use of "low-pass filter" is not exact enough for me.

For example, the Google Android API has this documentation:
http://developer.android.com/reference/android/hardware/SensorEvent.html#values

 public void onSensorChanged(SensorEvent event)
 {
      // alpha is calculated as t / (t + dT)
      // with t, the low-pass filter's time-constant
      // and dT, the event delivery rate

      final float alpha = 0.8;

      gravity[0] = alpha * gravity[0] + (1 - alpha) * event.values[0];
      gravity[1] = alpha * gravity[1] + (1 - alpha) * event.values[1];
      gravity[2] = alpha * gravity[2] + (1 - alpha) * event.values[2];

      linear_acceleration[0] = event.values[0] - gravity[0];
      linear_acceleration[1] = event.values[1] - gravity[1];
      linear_acceleration[2] = event.values[2] - gravity[2];
 }

How do I interpret that low-pass filter? What is the cut-off frequency? What is the transition bandwidth? Are they using this LPF solely to do averaging?

The filter in your example is a first-order infinite impulse response (IIR) filter. Its transfer function is:

$$
H(z) = \frac{1 - \alpha}{1 - \alpha z^{-1}}
$$

which corresponds to a difference equation of:

$$
y[n] = \alpha y[n-1] + (1-\alpha) x[n]
$$

where $x[n]$ is the filter input and $y[n]$ is the filter output.

This type of filter is often used as a low-complexity lowpass filter and is often called a leaky integrator. It is favored because of its simple implementation, low computational complexity, and its tunability: its cutoff frequency depends upon the value of $\alpha$. $\alpha$ can take on values on the interval $[0,1)$. $\alpha = 0$ yields no filtering at all (the output is equal to the input); as $\alpha$ increases, the cutoff frequency of the filter decreases. You can think of $\alpha = 1$ as a boundary case where the cutoff frequency is infinitely low (the filter output is zero for all time).

You can think of this intuitively by noticing that the filter input is weighted by $\alpha$, so as the parameter increases, the quantity $1-\alpha$ decreases, so each input sample has a smaller proportional effect on the value of any particular output sample. This has the effect of smearing out the filter's impulse response over a longer period of time. Summing over a longer period of time is similar to computing a long moving average. As the length of a moving average increases, the cutoff frequency of the average decreases.

For your example, where $\alpha = 0.8$, the frequency response of the filter is as follows:

From the example, I would guess that this filter is being used to smooth high-frequency noise out of a time series of measurements from a sensor, trying to tease out a comparatively low-frequency signal of interest. This would be a very typical application for this sort of filter.

On your other sub-question, you are correct that filtering is often implemented via convolution of the input signal with the filter's impulse response. In most cases, this is only done with finite impulse response (FIR) filters. IIR filters such as this one are typically implemented using the filter's difference equation; since an IIR system's impulse response is infinitely long, you must truncate it to some finite length to make convolution with it tractable, at which point the filter is no longer IIR. The difference equation format is almost always cheaper to implement computationally, although the feedback inherent in that structure can lead to numerical issues that must be addressed (such as internal overflow and roundoff error accumulation).


Solution for FFT size 1536 I have already asked this question 2 times, but still i am confused because this is not available anywhere on the internet and explained properly. I found some information in this document and i am going to apply their method , but confused to understand. here is the document http://www.altera.com/literature/an/an480.pdf. If anyone have some time to explain me kindly explain me where is radix 2 stage in this document to calculate the fft of size 512 because it has been broken into 3 small sequences of length 512 and then applied the fft on small lengths and then calculated it with combining and multiplying it with twiddle factors respectively. Kindly help me i am stuck in here for more than a week.

It's hard to see what could be added to Peter K.'s answer to make it any clearer. Here is a snippet of MATLAB code that shows the whole process in all it's glory. Just copy and paste in MATLAB or octave and you can step through it one line at a time.

%% 1536 FFT based on three FFTs of 512 each
n = 1536;
% Create a piece of noise
x = randn(n,1);
% calculate FFT using MATLAB native fft() function.
% We'll use this as a reference to prove it works
fx = fft(x);
% Break down into three signals of 512 points each
p = x(1:3:end);
q = x(2:3:end);
r = x(3:3:end);
% FFT each of those. This is  a 512 power-of-two standard FFT
fp = fft(p);
fq = fft(q);
fr = fft(r);
% Do three times periodic extention (just repeat it three times)
fp3 = [fp; fp; fp];
fq3 = [fq; fq; fq];
fr3 = [fr; fr; fr];
% calculate the 1536 twiddle factors
k3 = (0:n-1)';
W3 = exp(-j*2*pi*k3/n);
% assemble the result
fy3 = fp3 + W3.*fq3 + W3.^2.*fr3;
% calculate the error
ferror = fy3-fx;
fprintf('Error = %6.2f dB\n',10*log10(sum(ferror.*conj(ferror))./sum(fx.*conj(fx))));


Method for determining probe angle by analyzing skewed sine wave I've got a fun problem and would be curious to get feedback on how some of you would go about solving this.

Imagine I have a probe and am scanning the surface of some material. This material surface is described by a perfect cosine signal. If my probe is perfectly normal to the surface, my probe will output a "1 to 1" match of the surface signal (excluding some DC offset).

However, if my probe is not perfectly normal to the surface, the collected signal would appear to be a skewed sinusoid. I've put together a simple diagram of set up.

Currently, my method for determining the probe angle would be to determine the frequency content of the signal. In the case of the probe being perfectly normal to the surface, I would expect that frequency domain would show a single impulse corresponding to the surface frequency.

In the case of the angled probe, I would expect to a find a range of frequencies. Obviously the peak would still correspond to the dominant spatial frequency, however, there would now be a certain bandwidth associated the frequency domain. This bandwidth is associated with the fact that the collect signal is a skew sinusoid.

My question would then be, how would I correlate signal bandwidth to probe angle?It would be nice to get an analytical solution to signal bandwidth and probe angle.

Furthermore, if anyone else could see a way of determining probe angle that isn't frequency domain related, I would be interested in hearing it.

Edit: As I think about this problem more, I'm becoming less confident in whether or not probe angle would create a skewed sinusoid. Assuming the angle doesn't exceed the max slope found in the signal, I'm starting to think the resulting signal would simply be a phased offset sinusoid...

I will assume:

The probe is like a line, infinitely narrow.
We move the probe along the $x$ dimension, parallel to the large-scale surface, and from each position, probe the surface:
The probe is extended to the direction it is pointing to until it meets the surface. (An alternative assumption would be that probing is done to the direction $y$, that of the normal of the large-scale surface, in which case the waveform measured would be sinusoidal unless the probe stem meets the surface before the tip. I wonder if your edit describes this scenario.)

Equation of line starting at coordinates $(x_0, 0)$ and with a positive slope $k = \frac{\Delta y}{\Delta x} \ge 1$ (where $1$ is the maximum slope of the sinusoid of Eq. 2, to ensure continuity):

$$y = k (x - x_0)\tag1\\
\Rightarrow x = x_0 + \frac yk.$$

Equation of the first period of the constant-biased sinusoid, spanning $0 \le x \le 2\pi$:

$$y = 1 - \cos(x)\tag2\\
\Rightarrow x = \begin{cases}
\arccos(1 - y)&\text{if }0 \le x \le \pi\\
2\pi - \arccos(1 - y)&\text{if }\pi \le x \le 2\pi.
\end{cases}$$

By Eq. 2 the period starts at $(0, 0)$ and ends at $(2\pi, 0)$, which makes the corresponding period of the skewed sinusoid also start and end at these points. We obtain the waveform by solving $y$ as function of $x_0$ from the system of Eqs. 1 and 2. Unfortunately there is no closed form solution, but we can still make some plots:

Figure 1. The skewed waveform for different values of the constant $k$. At $k=\infty$, the waveform is equal to the constant-biased sinusoid. The waveform shown represents the probe tip $y$ position as function of the coordinate $x_0$ at which the probe shaft crosses the $y$ axis. If you desire the "probe length" waveform, you need to multiply this waveform by a factor $\tfrac{\sqrt{k^2 + 1}}{k}$.

I doubt that you will be able to find an analytical solution in frequency domain. What we can do is determine $k$ based on the waveform peak location $x_0 = x_p$ by setting $y = 2$, for which the system of Eqs. 1 and 2 gives:

$$k = \frac{2}{\pi - x_p}.\tag{3}$$

To get the angle $\phi$ by which the probe angle deviates (moderately) from the normal of the large-scale surface:

$$\phi = \operatorname{atan2}\left(1, k\right) = \operatorname{atan2}\left(\pi - x_p, 2\right).\tag{4}$$


convolution matrix for image scaling Is there any way how to compute convolution matrix for Nearest Neighbor (bicubic, bilinear) image scaling (upscale/downscale)?

This can't be done as a convolution because some of the operations you mentioned are not linear. Additionally, convolution by definition (multiplication by a toplitz matrix) shouldn't ever change the input and output size.

For bicubic:

http://en.wikipedia.org/wiki/Bicubic_interpolation#Bicubic_convolution_algorithm

Notice the high order terms.


pdf of the sum of gaussian distributions using fft Ask I am trying to derive the pdf of the sum of independent random variables. At first i would like to do this for a simple case: sum of gaussian random variables. I was surprised to see that i don't get a gaussian density function when i sum an even number of gaussian random variables. I actually get:

which looks like two halfs of a gaussian distribution. On the other hand, when i sum an odd number of gaussian distributions i get the right distribution:

below the code i used to produce the results above:

import numpy as np
from scipy.stats import norm
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
%matplotlib inline

a=10**(-15)
end=norm(0,1).ppf(a)
sample=np.linspace(end,-end,1000)
pdf=norm(0,1).pdf(sample)
plt.subplot(211)
plt.plot(np.real(ifft(fft(pdf)**2)))
plt.subplot(212)
plt.plot(np.real(ifft(fft(pdf)**3)))

Could someone help me understand why i get odd results for even sums of gaussians distribution? Thanks in advance

The FFT algorithm expects the origin of the signal to be on the leftmost sample of the signal. You are convolving two shifted Gaussian together, yielding a Gaussian that is even more shifted. Because of the FFT imposes periodicity, when the curve is shifted past the right edge it comes back in on the left edge.

Whit an odd number the shift is exactly the full width of the signal, so you don’t notice it.

A solution would be to use a zero-mean Gaussian distribution, which you can obtain from yours by ifftshift. After performing the convolution, you can shift it back to where it was using fftshift.

np.fft.fftshift(ifft(fft(np.fft.ifftshift(pdf))**2))


explanation on the mimo technique In MIMO systems we have:
$y=Hx+n$ with $H$ is a matrix $N_r \times N_t$.

But the previous equation is only true for non-frequency selective channels.

What is the formula for frequency-selective channels?

Normally, $y=h \star x +n$. But if H is a matrix, we can not the convolution. How to solve this problem

The input-otput relation for frequency selective channel:

$y_i=\sqrt{\frac{P}{M_T}}H_is_i+n_i$

More information can be found here.


From a Waveform to a .wav file I would like to make an autoencoder for music generation using PyTorch.
In input, I would like to have a Waveform (from a .wav file) and I want the output to be a Waveform.
My problem is that I don't know how to get a .wav file from a waveform that is the output...
Does anybody have an idea?

In MATLAB you can create a .wav file waveFileName.wav from a waveform y with sampling frequency Fs using MATLAB's audiowrite as follows:
filename = 'waveFileName.wav';
audiowrite(filename,y,Fs);

You can also read and output a waveform y and its sampling frequency Fs from the .wav file waveFileName.wav using MATLAB's audioread as follows.
[y,Fs] = audioread('waveFileName.wav');

I am sure there must be corresponding functions in Python.


Wavelet transform in MATLAB Suppose I have a wave with $20 \textrm{ kHz}$, $100 \textrm{ kHz}$ and $300 \textrm{ kHz}$. Sampling frequency used is $1000 \textrm{ kHz}$. I apply the discrete wavelet transform on the wave like dwt(wave,'db2'). I will get one level of approximation and detail coefficients. According to the basics, the detail coefficient will have  $300\rm k$ component and approximation coefficient will have the $100\rm k$ and $20\rm k$ components

But when I did fft on the output(on approximation and detail), I didn't get what I expected.

Could anyone post some MATLAB code with which I can verify it?
And also tell whether I am doing the procedure correctly?
Also if anyone could explain the practical side of this tool?

What you might be forgetting is that dwt does downsampling. After filtering (low-pass or high-pass), the filtered signal is subsampled. In order to keep the number of samples. So aliasing may occur, depending on the quality of the filters.

And db2 is a quite poor filter. So if you idwt either approx or details, by replacing the other coefficients by zeros, you almost get your coefficients in order: on the top, the original signal,  the two low frequencies in the second plot, the high frequencies in the third plot. Of course you have some aliasing, due to the poor resolution.

Try the attached code with db24, you will get much less artifact.

What happens now with the approximation and the details? My interpretation is that the approximation is now sampled at $500$ kHz, and the details too, except that it is now aliased to the base-band:

and they need to be unmixed by reciprocal filters of the idwt.

The FFTR.m code should be downloaded.

But why study sines with discrete  wavelets?

f = [20/1000  100/1000  300/1000];
nSample = 1024; time = (1:nSample)';
data = cos(2*pi*time*f);
dataSum = sum(data,2);
waveName = 'db2';
[a,d] = dwt(dataSum,waveName);
z = zeros(size(a));

dataSumD = idwt(z,d,waveName);
dataSumA = idwt(a,z,waveName);

timeAD = 1:2*length(a);
timeAD=timeAD(1:2:end);

figure(1)
subplot(3,2,1)
plot(time,dataSum);axis tight
subplot(3,2,2)
plot(FFTR(dataSum));axis tight

subplot(3,2,3)
plot(time,dataSumA);axis tight
subplot(3,2,4)
plot(FFTR(dataSumA));axis tight

subplot(3,2,5)
plot(time,dataSumD);axis tight
subplot(3,2,6)
plot(FFTR(dataSumD));axis tight

figure(2)
subplot(2,2,1)
plot(timeAD,a);axis tight
subplot(2,2,2)
plot(FFTR(a,2));axis tight

subplot(2,2,3)
plot(timeAD,d);axis tight
subplot(2,2,4)
plot(FFTR(d,2));axis tight

So you have the code, I do not know whether you did stuff correctly, and wavelets are quite useful if used widely on signals that are worth  it. Sines typically are not on the  list.


Is any time-domain vector stored on a digital machine always intrinsically causal? Reading up on causality, I understand the mathematical definition, in so far as that a causal system, is one where the output depends only on the current time, and possibly the past time, but never the future times.

I have also seen that we can 'make a (transfer) system causal' by 'shifting it in time by an appropriate amount'.

I do not necessarily understand the above statement, and would like to see a simple example of two of such a thing happening. (That is, plot it when its non-causal, and re-plot it when it is causal).

For example, if I have a length $N$ impulse response, stored in a vector on my machine, how does it 'know' that this is causal or non-causal? It seems as though the time-axis would always start at $n=0$, where the first point in my vector is.

Causality is not so much a characteristic of a signal as it is a characteristic of a system.  For example, a non-causal system can have an output at time $t$ which depends on the input at time $t+1$.  When thinking in terms of time, a non-causal system breaks our intuition because it has to "see the future" in order to operate.

Let's say that I want to create an audio effect that is a sort of "reverse echo."  In other words, I want to hear an echo of the sound before the sound event actually occurs.  This would be an example of non-causal processing because at any given time, the output depends on input which has not yet occurred. This is not a problem in the case of a recorded audio signal because we already have all of the time samples available to us, so we can "look ahead" and use this "future information" right now.

But what if I wanted to implement this "reverse echo" effect in a live performance?  I obviously can't "look ahead" to grab samples out of the future, but I can wait until all of the samples I need are available to me, and then apply the processing.  This would produce the exact same output signal as if I had done the non-causal processing mentioned above, but with one significant difference: my output would be delayed.

Given a signal $y$ which depends on input $x$ as

$$y[n] = a_2 x[n-1] + a_1 x[n] + a_0 x[n+1]$$

we can make a causal version of $y$, which we will call $y^\prime$, by simply delaying $y$ by one time step:

$$y^\prime[n] = y[n-1] = a_2 x[n-2] + a_1 x[n-1] + a_0 x[n]$$

We have not altered the signal in any way, besides to shift its time indexing. Any shift in phase can be seen as a direct result of the delay, and cannot be reversed by multiplying some phase term.


Why do we use window in time domain rather than do FFT modify the spectrum and than inverse FFT I thought that DSP would be done by using FFT of portions of a signal, modify the samples that result from FFT (since they represent the spectrum of our signal + noise) and remove any unwanted signals and than do an inverse FFT to get a time domain representation of the filtered signal (noise has been removed now). This however is not done, instead we do all the work in time domain using window functions. Why?

If we multiply window function in time domain than we are convolving the window function's frequency response with spectrum of our signal in frequency domain, how does that work out? I mean if we just did all the work in frequency domain by multiplying our signal with the frequency response of the filter, that would be like filtering right? But here we do all the stuff in time domain instead using window.

->Lets see where my confusion comes from. For analogue filters e.g low pass filter, we have this pulse like frequency response. When we filter a signal, we are effectively multiplying our signal's spectrum with the pulse like frequency response of the filter. This will reduce all frequencies in our signal above a cut-off to 0. This is how a low pass filter works in essence. Why not do the same with digital filters as well?

Windowing is used because the DFT calculations operate on the infinite periodic extension of the input signal. Since many actual signals are either not periodic at all, or are sampled over an interval different from their actual period, this can produce false frequency components at the artificial 'edge' between repeated intervals, called leakage. By first multiplying the time-domain signal by a windowing function which goes to zero at both ends, you create a smooth transition between repeated intervals in the infinite periodic extension, thus mitigating the creation of these artificial frequency components when we then take the DFT.

This paper gives a more in depth look at this phenomenon, as well as some insight into the effects of different windowing functions.


Is the Laplacian Filter an High Pass Filter (HPF)? Is this kernel name as mean difference kernel or Laplacian filter?

Is mean difference just a method?

There are many approximations for the Laplacian Filter (See The Hypermedia Image Processing Reference - Laplacian/Laplacian of Gaussian):

Indeed this is an High Pass Filter. Namely it will remove low frequencies (Specifically it will remove the DC Value, namely the output image will have mean value of 0).


Determine a sequence whose N point DFT is given without computing inverse For example, I'm given a sequence x(n) = { 3, 5, 1, 4, -3, 5, -2, -2, 4 } whose 9 point DFT is X[k]. Without computing IDFT, i need to find the sequence y(n) whose 9 point DFT is given by $Y[k] = (-1)^k ~ e^{j \frac{\pi}{3} k}$ . How do I approach this question?

Hints only:

Think about time-frequency duality of DFT pairs. What is the DFT of sequence $-1^k$?

And then there is a phase term in the DFT, so there should be a shift in discrete time sequence. Use time-delay property of DFT.

These 2 properties should give you the solution.


Designing high-pass filters using impulse invariance One of the main disadvantages of realizing digital filters using impulse invariance is aliasing. According to the Nyquist sampling criterion, in order for the frequency response of the digital filter to match that of the corresponding analog filter with impulse response $h(t)$, $h(t)$ has to be band-limited. That is, the analog filter has to be low-pass. Does this mean we cannot map high-pass analog filters to high-pass digital filters using impulse invariance? Or is it just a matter of mapping a low-pass analog filter to a low-pass digital filter with frequency response $H(e^{j\Omega})$, and then finding the frequency response of the corresponding high-pass filter simply by substracting it from 1?

You're right, you can't use the impulse invariance method for transforming continuous-time filters with a transfer function that is not strictly proper, i.e., for which the degree of the numerator is not strictly smaller than the degree of the denominator. High-pass and band-stop filters have transfer functions with numerator and denominator polynomials of the same degree, which means that the corresponding partial fraction expansion has a constant term. The latter corresponds to a (scaled) Dirac delta impulse in the impulse response, which obviously cannot be sampled.

Take a first-order high-pass filter as an example:

$$H(s)=\frac{s}{s+a}=1-\frac{a}{s+a}\Longleftrightarrow h(t)=\delta(t)-ae^{-at}u(t)\tag{1}$$

Now you can't determine $h(nT)$, $n\in\mathbb{Z}$, at $n=0$.

In sum, for the impulse invariance method to be applicable we require

$$\lim_{s\to\infty}H(s)=0\tag{2}$$


Why do FIR systems need to be antisymmetric about M/2 to be linear phase? Why do FIR systems need to be antisymmetric about M/2 to be linear phase?
where M=delay.Also im talking about discrete time.And why are hilbert transformers antisymmetric on n=τ?

Consider the $M$ th order FIR system, where $M$ is even
$$\begin{aligned}
H(z) &= \sum_{m=0}^{M-1}b_m z^{-1} \\
     &= b_0 + b_1z^{-1} + \ldots b_{M-1}z^{M-1} \\
H(e^{j\omega}) &= b_0 + b_1 e^{-j\omega} + b_2 e^{-2j\omega} + \ldots b_{M-1}e^{-j(M-1)\omega} \end{aligned}$$
Now if our coefficients are antisymmetric, then $b_m = -b_{M-1-m}$, or, $b_0 = -b_{M-1}, b_1 = -b_{M-2}, \ldots, b_{\frac{M}{2}} = -b_{\frac{M}{2}-1}$. This gives us
$$\begin{aligned}
H(e^{j\omega}) &= b_0 - b_{M-1}e^{-j(M-1)\omega} + b_1 - b_{M-1}e^{-(j(M-1) + 1)\omega)} + \ldots \\
&= e^{-j(M-1)\omega} \left[ 2b_0 + b_1(e^{-j\omega}-e^{j\omega}) + \ldots b_{\frac{M}{2}} (e^{-j\frac{M}{2}\omega}-e^{j\frac{M}{2}\omega}) \right] \\
&= -2je^{-j(M-1)\omega}) \left[jb_0 +b_1 \sin \omega + \ldots + b_{\frac{M}{2}} \sin{\frac{M}{2}\omega} \right] \\
&=
-2je^{-j(M-1)\omega})\left[jb_0 + \sum_{m=1}^{M/2} b_m \sin(m\omega)\right]
\end{aligned}$$
We see that $H(e^{j\omega})$ is linear phase, with its phase $ = -(M-1)\omega - \frac{\pi}{2}$. This would have not worked out if the filter coefficients were not antisymmetric or symmetric.


The Fourier transform of sinusoids' products with possible other components I know that in general it transforms the signal from the time to the frequency field but these specific cases seem pretty demanding. Do I calculate each part separately and then just leave them with convolution between them? Or do I have to calculate any integrals?
\begin{align}
&\frac{1}{t^2}\cdot\cos(2πt)\cdot\cos(2πt)\\
&8\cos(20πt)\cdot\cos(40\pi{t})\cdot\cos(80\pi{t})
\end{align}
For example for the second one will the result be
$$\bigg(4\big[\delta(f-10)+\delta(f+10)\big]\bigg)\star \bigg(\frac 12\big[\delta(f-20)+\delta(f+20)\big]\bigg)\star\bigg(\frac 12\big[\delta(f-40)+\delta(f+40)\big]\bigg)$$
Is that correct?

HINT:
It is easy to see that things can be simplified using the trigonometric product-to-sum identity in Equation $(1)$ below:
$$
\cos(\alpha)\cos(\beta) = \frac 12\big[\cos(\alpha+\beta) + \cos(\alpha-\beta)\big]\tag{1}
$$

In the first example
$$
\frac{1}{t^2}\cdot\cos(2πt)\cdot\cos(2πt) = \frac 1{2t^2}\big(\cos(4\pi t) + 1\big)
$$
From there you would want to visit the cosine modulation frequency-shift and differentiation properties of the Fourier transform, here in $(2)$ and $(3)$ respectively.
\begin{align}
\mathcal F\big\{x(t)\cos(2\pi f_0 t)\big\} &= \frac 12\big[X(f - f_0) + X(f + f_0)\big]\tag{2}\\
\mathcal F\left\{\frac{d^n x(t)}{dt^n}\right\}&=  \left(j2\pi f\right)^nX(f)\tag{3}
\end{align}

In the second example, using Equation $(1)$ you then have
\begin{align}
8\cos(20\pi t)\cos(40\pi t)\cos(80\pi t) & = 8\bigg(\frac 12\big(\cos(60\pi t) + \cos(20\pi t)\big)\cos(80\pi t)\bigg)\\
& = 4\big(\cos(60\pi t) + \cos(20\pi t)\big)\cos(80\pi t)\\
& = 2\big(\cos(140\pi t) + \cos(20\pi t)\big)\\
&\quad + 2\big(\cos(100\pi t) + \cos(60\pi t)\big)\\
\end{align}
With this you simply have sums and you don't have to think of convolutions, you have individual sinusoids at frequencies $10\ \rm Hz$, $30\ \rm Hz$, $50\ \rm Hz$, and $70\ \rm Hz$.


Effect of noise and entropy Relation of Entropy and SNR : Based on this question and answer, I had another question that struck me and I am curious to know, if somebody can shed some light, on the following situation: $y= desired_{signal} + noise$ is received by the receiver and the $error = y - \hat{y}$ where $\hat{y}$ is simulated at the receiver by using guessed parameters with the knowledge of the kind of model used at the transmitter.

For increasing SNR (attenuating noise) would Entropy of error decrease? In the link provided, the entropy of the transmitted signal decreases with increasing SNR since the uncertainty decreases when noise is getting attenuated. But what about the entropy of error? Should error entropy increase with increasing SNR? I am confused about this.

EDIT: The way I calculated entropy of error:
Let the system model be AR(2): $y(t)= a1y(t-1) + b1y(t-2) + noise$. At the receiver end, I have $z(t) = y(t) - (a2y(t-1) + b2y(t-2))$ where $(a2,b2)$ are close guesses to ($a1,b1$). I calculated 2 entropies H1 for y(t)  and H2 for z(t).

z(t) will be close to zero when a2,b2 will be equal to a1,b1 in which case entropy of z(t) is found to be minimum among all different z(t) for a particular SNR value of noise.
 So, for a particular noise value, I have 10 pairs of (a2,b2) and for each I get 10 values of H2(z(t)). I chose the minimum among all H2 for a particular SNR level. But, for different noise levels, as I increase the SNR at the transmitter end, the minimum entropy, H2(z(t)), increases with increasing SNR.
I found that H1 decreases with increasing SNR but the reverse happens for H2.
Is this trend correct?

According to the link you provided:

$H(error) = H(y-\hat{y}) = H(Var[y]+Var[\hat{y}])$. While $Var[y]$ decreases with increasing of SNR, $Var[\hat{y}]$ stays the same (because it's not dependent of the noise, it depends on the signal alphabet only). So we may guess $H(error)$ decreases with increasing of SNR.


What's wrong with this Whittaker-Shannon-Kotel’nikov interpolation implementation? I tried to implement Whittaker-Shannon-Kotel’nikov interpolation formula but I get unexpected results: the reconstructed signal lags with respect to the original.

I know that I can not expect a perfect reconstruction because I am sampling a finite time signal (just one cycle of a sine wave) and this finite time signal is not bandwidth limited and so the Whittaker-Shannon-Kotel’nikov theorem does not hold... but I do not know whether to expect a lag.

I suspect there is something wrong in my implementation (for example in a previous implementation I was erroneously using $\frac{\sin x}{x}$ instead of $\frac{\sin \pi x}{\pi x}$) but I can not find it.

f=1;
T=1/f;
mysignal = @(t) sin(2*pi*f*t);

f_s=5*f; % sampling rate
T_s=1/f_s; % sampling period
t=0:T_s:T; % times at which mysignal is sampled
t1=0:.05*T_s:1.5*T; % times to plot mysignal "continuosly"

% plot mysgnal as a "continuos" signal
hold off;
plot(t1,mysignal(t1),'r')
hold on;

function y = Whittaker_Shannon_interpolation(samples,T_s,t)
% samples are the samples of the original signal;
% T_s is sampling period;
% t are the times at which interpolate the original signal;
% y are the interpolated values (reconstructed) of the original signal at times t.
  y = zeros(1,length(t));
  normalized_sinc = @(x) sin(pi*x)./(pi*x);
  for k = 1:1:length(samples)
    % https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula
    % the above formula uses so called "normalized sinc", i.e. sin(pi*x)/(pi*x)
    y+=samples(k)*normalized_sinc((t-k*T_s)/T_s);
  endfor
endfunction

samples=mysignal(t);
stem(t,samples,'b')
hold on;

y=Whittaker_Shannon_interpolation(samples,T_s,t1);
plot(t1,y,'k-o')

legend('original signal','samples','reconstructed')

Replace the line

y+=samples(k)*normalized_sinc((t-k*T_s)/T_s);

by

y+=samples(k)*normalized_sinc((t-(k-1)*T_s)/T_s);

because the indices inside the sinc function must start from $0$, not from $1$. The first sinc must have its peak at $t=0$, not at $t=T_s$.


What happens to the standard deviation of a white Gaussian noise after low pass filter? If I apply an ideal box-car low pass filter on a random signal with white Gaussian noise, what happens to the Gaussian distribution's $\sigma$?
I know that the auto-correlation function becomes a sinc(), but not sure where to go from here...
Thanks!

I assume you are working with discrete-time, since continuous-time white noise has infinite power ($\sigma^2$).

First, remember that the power of a stationary process is always equal to the autocorrelation at 0 ($P_x = R_x[0]$); and the variance is the autocovariance at 0 ($\sigma^2_x = \rm{Cov}_x[0]$). These 2 expressions are equal for processes with 0 mean, which is the case for white noise. So, if you already calculated the autocorrelation, you have the power.

But you don't need to calculate the autocorrelation of the filtered signal. You can use Parseval's property and calculate the variance in frequency domain:

$\sigma_y^2 = \frac{1}{2\pi}\int_{2\pi} |H(e^{j\theta})|^2 G_x(e^{j\theta}) d\theta = \frac{1}{2\pi}\int_{-\theta_0}^{+\theta_0} 1 \times \sigma^2 d\theta$ since the ideal lowpass filter is 1 between $\pm\theta_0$, and the spectral density of the input signal is $\sigma^2$ (white noise).

The result is very easy to calculate, and gives $\sigma_y^2 = (\theta_0/\pi) \sigma^2$.


Scale of Amplitude in FFT (Matlab) Using Matlab's FFT function, I am getting different Amplitude values in a particular frequeny, when I use different range of Bandpass filters centering around that particular frequency. What's wrong here? How can a particular Amplitude scale be selected?

In the following code, I calculated two FFTs of the same signal over the frequency range of (59Hz to 61Hz) and (59.98Hz to 60.02Hz) after bandpassing them in the said range.

[x_noise,fs] = audioread(h);
x_noise = x_noise(1:length(x_noise)/360);
xlen = length(x_noise);
x_main=x_noise;

i=1;
%%
%Noise
f1=59;
f2=61;
fn=fs/2;
w1=f1/fn;
w2=f2/fn;

[b,a]=butter(5,[w1 w2],'bandpass');
xn_har = filter(b,a,x_noise);

%FFT of Noise
X_Noise = abs(fft(xn_har));
 p = linspace(0,1000,length(X_Noise));
 figure(1)
 plot(p,X_Noise)

%%
%Signal
f1_sig=59.98;
f2_sig=60.02;
w1_sig=f1_sig/fn;
w2_sig =f2_sig/fn;

[b,a]=butter(5,[w1_sig w2_sig],'bandpass');
x_sig_har = filter(b,a,x_main);

%FFT of Signal
X_sig = abs(fft(x_sig_har));
q = linspace(0,1000,length(X_sig));
figure(2)
plot(q,X_sig)

But after doing FFT, the Amplitude in 60Hz seems to vary a lot. As seen in the image, one has a value of around 1600, while the other has around 1.2. Why is this happening, they are at the same frequency?

The problem is likely to be your filter. Designing a very narrow filter is very hard; a very large order is required and even then, some responses may simply be unrealizable. Note that "narrow" here is relative to the Nyquist frequency. I don't know what your sampling frequency is, but let's assume 8000 Hz. Run this code:

fn=4000;f1=59;f2=61;w1=f1/fn;w2=f2/fn;
b=butter(5,[w1 w2],'bandpass');
freqz(b);

The command freqz plots the actual frequency response of the filter that was returned by butter. It is very important always to verify that your filters' response resemble what you had in mind -- that is, that they meet your requirements.

I got slightly better results using fir1 instead of butter. However, in your case, what I would try to do is to downsample the signal so that, relative to the Nyquist frequency, the filter is not too narrow. To see the effect, run this code:

fn=120;
f1=59;f2=61;w1=f1/fn;w2=f2/fn;
b=fir1(200,[w1 w2],'bandpass');
freqz(b)

You'll see a nice, narrow response with 0 dB gain at the center frequency and around -70 dB rejection out of band.


Why does gamma correction of image use two different definitions of gamma? I've noticed that the gamma value used in level adjustment of Photoshop and Canvas is 'opposite' of gamma values defined in Wikipedia or MATLAB imadjust documentation.

In Wikipedia or MATLAB, high gamma (>1) will make an image darker, whereas in Photoshop high gamma (>1) will make an image brighter.

I've also noticed that, even among image processing software, Zeiss' Zen software uses high gamma (>1) to make an image darker.

There seem to be two definitions of gamma value. This is just confusing... I assume that this is related to 'gamma encoding' and 'gamma correction'. But can anyone explain the rationale of these two defintions and when to use which one?

Vout = Vin ^gamma % Wikipedia or MATLAB

Vout = Vin ^1/gamma % Photoshop or Canvas

The concept of the gamma (and its correction) has emerged as a means of controlling and correcting the nonlinear voltage-brightness characteristics of the legacy input cameras and legacy display CRTs which are both analog devices.

Both of them had nonlinear relationships between the optical brightness levels being captured (or displayed) and the corresponding voltage levels that result of (or driving) it.

In order to reproduce the image brightness levels (and RGB channel intensity levels as well) as faithfull as possible to the original content, the exponential distortions introduced by the CRT displays are pre-corrected by the studio processing equipment before the content is being transmitted for broadcasting. The correction being applied at the studio is the inverse of the distortion being created at the CRT display, hence they will cancel each other, yielding a faithful, brightness neutralized, image at the CRT output.

Furthermore, for adequately displaying video content recorded with different standards and/or captured with different gamma correction parameters, the TV sets and monitors are available with gamma selection options as well.

Hence you can apply a gamma correction either at the camera input stage before it's transmitted or gamma adjustment at the display (CRT) stage during it's played back. So the difference you percieve possibly refers to the stage at which it's meant to be applied. So as you have stated, the same gamma parameter value will create two opposing results on the same image according to which stage it was applied to. Or stated in other words, a given value $\gamma$ will create the same effect with $1/ \gamma$ when applied at opposite sides.


Are SURF descriptors normalized with respect to orientation? I know that both SIFT and SURF descriptor vectors are normalized with (unit vectors) to deal with illumination changes, but as far as I can see, the SIFT descriptors are also normalized with the major orientation of the keypoint.

Does anyone know if the SURF descriptor does this, as I cant find anything about this in the SURF paper by Herbet et al?

I would appreciate any help with this!
I am currently using SURF descriptor as included with matlab computer vision toolbox.

The SURF descriptor implicitly accounts for orientation of the interest points based on the distribution of edge intensities around it. So yes it is invariant to rotation. This is also mentioned in the description of the SURFPoints class:

  Orientation

  Describes the orientation of the detected feature. This value must be
  specified as an angle, in radians. The angle is measured from the
  X-axis with the origin at the point given by the Location property.
  Typically, this value gets set during the descriptor extraction
  process. The extractFeatures function modifies the default value of 0.
  Do not set this property manually. Rely instead on the call to
  extractFeatures to fill in this value. The Orientation is mainly
  useful for visualization purposes.

  Default: 0.0


Aliasing Square Wave I'm having a really weird issue with generating waveforms from frequencies C1 to C8.

When I generate sine waves, everything sounds good; no problems whatsoever. However, when I generate a square wave... BAM weird aliasing issues on the higher frequencies (around C6 to C8). I have no idea what I'm doing wrong at this point.

I'm going to take a wild guess here but I'm thinking it might be because square waves have a lot of harmonics so it's not producing the higher harmonics correctly when I reach a higher frequency... thus the really high harmonics are getting aliased. Would like someone to confirm with me/guide me to a solution to this. Is it possible to try to implement a filter or is there a better solution ?

Thanks!

I'm setting my sample rate to 44.1 kHz and here are my functions for waveform generation:

float generateSin {
    sample = sinf(_phs);

    _phs += _phs_incr;
    _phs = wrapPhase(_phs);
}

float generateSqr {
    sample = sinf(_phs);
    if (sample > 0) sample = 1;
    else if (sample < 0) sample = -1;

    _phs += _phs_incr;
    _phs = wrapPhase(_phs);
}

My guess is there's a problem with your wrapPhase function.

If I use the R implementation of it below (C8 shown), then I get all sorts of funny results.  That's the second image.  Removing it from the iteration (but including it in the calculation), makes for much cleaner signals (the first image).

R Code Below

#30582

wrapPhase <- function(phi)
{
  return(phi %% 2*pi)
}

# https://en.wikipedia.org/wiki/Piano_key_frequencies
cs <- c(32.7032, 65.4064, 130.813, 261.626, 523.251, 1046.50, 2093.00, 4186.01)
c1 <- cs[1]
c8 <- cs[8]

fs <- 44100

phs <- rep(0,T + 1)
phs_incr <- c8 / fs * 2 * pi

T <- 1000

sample <- rep(0,T)
sqr_wave <- rep(0,T)

for (t in seq(1,T))
{
  sample[t] <- sin(wrapPhase(phs[t]))

  if (sample[t] > 0)
  {
    sqr_wave[t] <- 1
  }
  else
  {
      sqr_wave[t] = -1;
  }

  print(sample[t])

  phs[t+1] <- phs[t] +  phs_incr
  #phs[t+1] <- wrapPhase(phs[t+1]);
}

phi <- seq(-pi, pi, 0.001*pi)
phi_wrap <- rep(0, length(phi))

for (k in 1:length(phi))
{
  phi_wrap[k] <- wrapPhase(phi[k])
}

plot(sample[1:50], col="blue", type="l")
lines(sqr_wave[1:50], col="green")


Proving Convergence of a limit related to Cross-Spectrum Estimation Proving asymptotic covariance of Cross-Spectrum Estimation I came across this mathematical problem:

Consider two constants $M,L$ such that  $M\to \infty, L \to \infty \hbox{ and } \frac{L}{M} \to 0$, And $\ell,s,r \in \mathbb{Z}$ $\nu_1,\nu_1 \in \mathbb{R}$

$$ \lim_{M\to +\infty} \left\{\frac{1}{L}\sum_{\ell=-M}^{M} U_M(r,\ell,s) w\left(\frac{\ell+s}{L}\right)\  w\left(\frac{\ell}{L}\right) e^{-2\pi i\ell(\nu_1-\nu_2)}\right\}$$

where  $r$, $s$ are constants and:

$$ U_M(r,\ell,s)= \left\{\begin{matrix}
0 &  r<-M+\ell+s\\
1-\frac{\ell+s+r}{M}& -M+\ell+s\leq r \leq 0 \\
1-\frac{\ell+s}{M}& 0 < r \leq s\\
1-\frac{\ell+r}{M}& s<r \leq M-\ell\\
0 &  r>M-\ell
\end{matrix}\right.$$
And $w(x)$ is a function which is even, bounded,  square-integrable and such that $w(0)=1$.

It has been proven [1] that this limit exists and it is equal to  0 for $\nu_1\neq\nu_2$ and :
$$=\int_{-\infty}^{\infty}w^2(x)dx $$
for $\nu_1=\nu_2$.
I am trying to prove that replacing in the limit $U_M(r,\ell,s)$ by

$$V_M(r)= \left(1-\frac{|r|}{M}\right)$$

Produces the same result but I don't have the background to understand the mentioned proof, is there a simple way to prove this claim? I am only interested in the case of $V_M(r)$.

[1] Hannan, E. J. Multiple time series, vol. 38. John Wiley & Sons, 2009. Proof of Theorem 9 chapter V, page 314.

A first thing to try is to reduce the unwanted variables. What happens if you set  $\ell=s=0$?

$$ U_M(r,0,0)= \left\{\begin{matrix}
0 &  r<-M\\
1-\frac{r}{M}& -M\leq r \leq 0 \\
1-\frac{0}{M}& 0 < r \leq 0\\
1-\frac{r}{M}& 0<r \leq M\\
0 &  r>M
\end{matrix}\right.$$

Luckily, this definition seems to fit with $V_M$, since the third inequality never happens. So your $V_M(r)$ is just a special case. Nothing more to do.


What determines the accuracy of the phase result in a DFT bin? What are the factors that affect the accuracy and precision for the phase that is given by the DFT?

Just thinking medium-hard about this, it occurs to me that it must have something to do with the window length $L$ and the sampling frequency $F_s$, and maybe the fundamental frequency that I'm interested in extracting?

I'm going to keep working on some of the math, but I'd like to start a discussion and get some insight for which way to look.

If we refer to "phase resolution" as how many decimal places we can use for an estimate of the phase (such as the phase reported by one bin of the FFT), then
this phase resolution is a time domain characteristic, and we would need to use the statistics of the signal involved in order to determine the useful precision. We cannot provide this experimentally based on the having one sample result alone; this would be equivalent to providing an average and then from that one average result alone, being able to determine its accuracy and precision, which is not possible without knowledge of the statistics of the signal.

Each bin in the frequency domain reports the magnitude and phase for that particular frequency. The OP used this as helpful analogy in his comments to help explain what it is he is really looking for, and specifically being able to present a precision on that phase estimate, and looking for the relationship between the time observation and the ability to resolve phase (increase precision).

For context with the rest of my response, let us first review the DFT computation and specifically for a single frequency bin which provides a magnitude and phase result. One bin in an FFT is the result of a time-domain correlation of the signal with the specific frequency of that bin. A time domain correlation is simply an average of the complex signal, and in this case after we have frequency shifted the signal back to DC (this is what occurs in the DFT computation).

For one bin $k\omega_o$ which is at a multiple $k$ of the fundamental frequency $\omega_o$ which is $2\pi/T$ with T being the duration of the time domain signal (or in Hz simply $1/T$) given as:

$$G(k) = \frac{1}{N}\Sigma_{n=0}^{N-1} g(t)e^{-j k \omega_o n} $$

We see that to compute the complex value of that bin we first frequency shift (by de-rotating) the spectral content of that bin and then do a simple low pass filter by averaging.

Thus if we want to know the accuracy of that result, we want to understand the accuracy of any time domain signal after the result of averaging. With that I respond:

The OP well understands the relationship between frequency resolution and the duration of the time domain sample that was used in the DFT, which is simply $1/T$ where T is the duration when we use an equivalent noise bandwidth in the frequency domain (to state simply, two distinct frequency tones will appear as one DFT bin if the frequency spacing is significantly less than $1/T$, and distinct tones when the spacing is significantly more than $1/T$.) This does suggest loosely how many decimal places we could use to report what frequency a single tone is at in the DFT. That frequency result is the correlation over that duration to the frequency in that particular bin, and is an averaged result for that duration used.

The phase, being the derivative of frequency, will have a slope vs time as derived directly from the frequency resolution previously described, so we can readily quantify the precision of that slope, and this result is the mean slope over the duration of observation.

If we want to consider the accuracy of the magnitude and/or phase of any one sample (one bin in the frequency domain), we must then resort to statistical measures and understand how that phase may have deviated from the linear slope over the observation interval.

The magnitude and phase accuracy of that estimate are similar quantities describing the complex value for that bin, so the phase accuracy is more consistent with the magnitude accuracy and perhaps easier to conceptualize in how we define the accuracy of the magnitude. The accuracy of both are similarly driven by the duration of observation, and we would typically approach this as a statistical problem based on the temporal characteristics of the signal and the noise. Quantization of the vertical axis is a contributor, but other noise sources may be present and we also must consider the stationarity of the signal in deriving such estimates.

A VERY useful and graphical statistic demonstrating phase and frequency accuracy vs observation interval is the Allan Variance (and Allan Deviation, or affectionately known as ADEV). This along with a similar metric TDEV (time deviation) is a widely used tool in the clock community (atomic clocks, frequency standards, GPS, etc) and I will provide a brief example from the example graph of ADEV shown below that I grabbed from http://leapsecond.com by googling ADEV images (a great site by the way if you want to dive deaper in into ADEV).

To attempt a brief explanation of this plot for my salient points, the horizontal axis shows the averaging time $\tau$, and the vertical axis shows an accuracy metric in terms of fractional frequency (the frequency error divided by the frequency of the clock). If the signal was stationary and the noise white, the accuracy would improve at a rate of $1/\sqrt{\tau}$ which the plot is closer to on the left side with shorter averaging time. This does immediately tell us over what observation intervals we could use such simplifying statistics as a white and stationary process. The plot tells me I could average such a signal for up to 500 seconds and continue to improve my accuracy of the estimate through averaging (processing gain etc), and in a fixed point system as long as I increased my precision accordingly (think extended precision accumulators!). It also immediately tells me that if I averaged for 500 seconds, I would not get the full SNR processing gain that I would expect (the famous $10Log(BW_{1}/BW_{2})$ equation since the result is not 1E-12 at 100 seconds and therefore the trend came short of $1/\sqrt{\tau}$. But look at this interesting point- If I continued to average past that, my result would be less accurate!! This has such wide ramifications across all signal processing, such as using an equalizer that has a span that exceeds the delay spread of the channel, or too many taps in an FIR filter, or the limit of processing gain in a spread spectrum system such as GPS, and many other implementations where we are ultimately doing some sort of weighted average.

That said, the point is to determine an accuracy of the estimate of a sample (whether we use phase or magnitude parameters or both when describing complex samples), we must understand the statistics of the signal and noise involved, and then proceed with statistical approaches in communicating precision and confidence. For example, if we knew the the process was white and stationary, and we knew what the phase over all the frequency bins should be, we could then use the error terms for every bin to estimate the accuracy of any one bin. Without knowing that, if the system was ergodic and stationary, we could use successive fft captures to determine a statistical estimate of accuracy. (And as explained, there is often a shorter duration of time in which it is valid to make such an estimate, therefore in our favor to make the captures as soon as possible with minimum delay in between). If the system is white, then each sample is independent of the next, and therefore shifting by one sample and taking a new FFT would be a completely independent result; we can do a sliding FFT to be able to get many more captures in less time and further improve the confidence of our result). The point is in all cases we have no way of defining the accuracy from one sample alone unless deduced from known signal statistics.


Image compression using SVD in MATLAB I have to compress an image using SVD. I did a few things, but all of them produced a black & white image. However, the output must be RGB. How can I compress an image without changing its original format?

The typical thing to do is the low-rank approximation on separate channels. Assume that $C$ is a channel of the RGB image $I$:

rank = 10;
[U,S,V] = svd(C);
L = U(:,1:rank) * S(1:rank, 1:rank) * V(:, 1:rank)';

Now, L should be the compressed image. If you do this operation and compose the channels back, you should get a compressed RGB image.

However, such a method in my opinion is only good for mathematical understanding. It is not very practical due to the fact that it depends on image dimensions too much and cannot generate high SNR results. Or better stated, as you decrease rank, you drastically reduce the quality of the reconstruction. For this reason, one has to retain around ~50-100 components. Better methods doing this blockwise, or formulating the compression as an energy functional exist. Image denoising is a huge area of its own.

You might want to check some results with varying rank in this presentation.

Also, a similar question has already been asked; please follow the Stack Exchange discussion here.


FFT the whole input vs FFT segment suppose I have a 4 seconds signal captured at $f_s$ = 1024Hz, which gives 4,096 data array. What is the different between these 2 FFT methods (MATLAB code)

[1] FFT the entire input with 4096-points FFT

N = length(input);
Y = fft(input)/N;

[2] Divide the input into 8 segments (4096/8=512 data point each), and feed each segment to 512-points FFT, and then average the output:

N = length(input);
Y = zeros(1,512);
for i = 1:8
  segment = input(512*(i-0)+1 : 512*i);
  Y = Y + abs(fft(segment,512))/512/8;
end

Question #1: Are these 2 methods equivalent?

Question #2: Is the frequency resolution for both methods equal to $\frac{f_s}{N}$ = 0.25Hz per bin?

Question #1: Are these 2 methods equivalent?

No, method #1 is taking an FFT of the entire sequence resulting in frequency bin widths of fs/N = 1024/4096 = 0.25 Hz. Method #2 is summing shorter FFTs at a frequency resolution (bin width) of 1024/512 = 2.0 Hz. For #2, if you divide by the number of FFTs after you sum them you are then averaging them. This is often times used to reduce the variance of the result.

  Question #2: Is the frequency resolution for both methods equal to 0.25Hz per bin?

No, they are 0.25 & 2.0 Hz. The resolution, or bin width, of the FFT is fs/N.


Proving Nyquist Sampling Theorem for Strictly Band Limited Signals (Whittaker Shannon Interpolation Formula) I understand that the Nyquist sampling theorem dictates that the minimum  sampling frequency, $f_s$, be s.t. $f_s > 2B$, where $ B $, is the bandwidth of the signal. I have read the explanation for what happens when the input signal contains an impulse/non zero ferquency spectrum at $f_s = f_m$. But is there any way to do a mathematical proof that sampling at $f_s = 2B$ works when the signal is strictly bandlimited.
A little background : This is a HW problem

The given hint was to use the fact that "one way to interpret Nyquist
Sampling theorem is to note that any bandlimited signal can be represented as a superposition of bandlimited signals that are orthogonal to each other."

Approaching The Sampling Theorem as Inner Product Space
Preface
There are many ways to derive the Nyquist Shannon Sampling Theorem with the constraint on the sampling frequency being 2 times the Nyquist Frequency.
The classic derivation uses the summation of sampled series with Poisson Summation Formula.
Let's introduce different approach which is more similar to function analysis - Building an orthogonal space and using projection for analysis and synthesis.
Forming Orthonormal Basis
In this section we'll define an orthonormal base and derive the decomposition and composition process.
Definitions
First one should define the space of Band Limited Functions. The space of Band Limited Functions is defined by:
$$ \mathcal{B}_{ {W}_{s} } = \left\{ f \left( x \right) \mid F \left( \omega \right) = \mathscr{F} \left\{ f \left( x \right) \right\}, \; F \left( \omega \right) = 0, \; \forall \, \left| \omega \right| > {W}_{s} \right\} $$
In words, it means that for each function $ f \left( x \right) \in \mathcal{B}_{ {W}_{s} }$ its fourier transform $ F \left( \omega \right) = \mathscr{F} \left\{ f \left( x \right) \right\} $ vanishes for frequencies $ \left| \omega \right| > {W}_{s} $.
The inner product in this space is given by:
$$ \langle f \left( x \right), g \left( x \right) \rangle = \frac{1}{T} \int_{- \infty}^{\infty} f \left( x \right) g \left( x \right) dx, \; T = \frac{2 \pi}{ {W}_{s} } $$
One could easily show that this is a indeed an inner product space with valid inner product.
The main claim is the orthonormal basis of this space is given by:
$$ f \left( x \right) = \operatorname{sinc} \left( \frac{ x - n T }{T} \right) $$
Where $ \operatorname{sinc} \left( x \right) $ is the Normalized Sinc Function given by $ \operatorname{sinc} \left( x \right) = \frac{ \sin(\pi x) }{ \pi x } $.
The basis functions are parameterized by the parameter $ n $. Basically we have shifted and scaled function as the basis.
Proof of The Orthonormal Property
One must show the orthonormal property of the basis under the defined inner product:
$$ \begin{aligned}
\langle f \left( x \right), g \left( x \right) \rangle & = \frac{1}{T} \int_{- \infty}^{\infty} f \left( x \right) g \left( x \right) dx = \frac{1}{T} \int_{- \infty}^{\infty} \operatorname{sinc} \left( \frac{ x - n T }{T} \right) \operatorname{sinc} \left( \frac{ x - m T }{T} \right) dx && \text{} \\
& \overset{1}{=} \frac{1}{T} \int_{- \infty}^{\infty} \left( \frac{1}{2 \pi} \int_{- \infty}^{\infty} T \, \Pi \left( \frac{ \omega }{ {W}_{s} } \right) {e}^{-j \omega n T} {e}^{j \omega t} d\omega \right) \operatorname{sinc} \left( \frac{ x - m T }{T} \right) dx && \text{} \\
& \overset{2}{=} \frac{1}{T} \int_{-\infty}^{\infty} \frac{1}{2 \pi} T \, \Pi \left( \frac{\omega}{ {W}_{s} } \right) {e}^{-j \omega n T} \left( \int_{- \infty}^{\infty} {e}^{j \omega x} \operatorname{sinc} \left( \frac{x - m T}{T} \right) dx \right) d\omega && \text{} \\
& \overset{3}{=} \frac{1}{T} \int_{-\infty}^{\infty} \frac{1}{2 \pi} T \, \Pi \left( \frac{\omega}{ {W}_{s} } \right) {e}^{-j \omega n T} T \, \Pi \left( \frac{-\omega}{ {W}_{s} } \right) {e}^{j \omega m T} d\omega && \text{} \\
& \overset{4}{=} \frac{T}{2 \pi} \int_{ - \frac{ {W}_{s} }{2} }^{ \frac{ {W}_{s} }{2} } {e}^{j \omega \left( m - n \right) T} d\omega && \text{} \\
& \overset{5}{=} \begin{cases}
1 & \text{ if } m = n \\
0 & \text{ if } m \neq n
\end{cases}
\end{aligned} $$
Where:

Since $ \mathscr{F} \left\{ \operatorname{sinc} \left( \frac{ x - n T }{T} \right) \right\} = T \, \Pi \left( \frac{ \omega }{ {W}_{s} } \right) {e}^{-j \omega n T} $.
Changing order of integration for converging integrals.
Applying $ \mathscr{F} \left\{ \operatorname{sinc} \left( \frac{ x - m T }{T} \right) \right\} \left( - \omega\right) $.
Integration boundaries according to the Rect function (Multiplication).
Integration over a cycle or over a constant.

Since we proved the suggested base is indeed an orthonormal basis of the space the following holds:
$$ \forall f \left( x \right) \in \mathcal{B}_{{W}_{s}} , \; f \left( x \right) = \sum_{n = -\infty}^{n} \langle f \left( x \right), {g}_{n} \left( x \right) \rangle {g}_{n} \left( x \right) $$
Where $ {g}_{n} \left( x \right) = \operatorname{sinc} \left( \frac{ x - n T }{T} \right) $ and $ \langle f \left( x \right), {g}_{n} \left( x \right) \rangle $ is the projection of $ f \left( x \right) $ on $ {g}_{n} \left( x \right) $.
Projection Process
As written above, using the result of a projection of a function in the space onto the basis one could reconstruct it as:
$$ f \left( x \right) = \sum_{n = -\infty}^{n} \langle f \left( x \right), {g}_{n} \left( x \right) \rangle {g}_{n} \left( x \right) $$
The question is, what's is the projection of a general function in this space? Well, it turns out it can be shown in a sloed form way:
$$ \begin{aligned}
\langle f \left( x \right), {g}_{n} \left( x \right) \rangle & = \frac{1}{T} \int_{- \infty}^{\infty} f \left( x \right) {g}_{n} \left( x \right) dx = \frac{1}{T} \int_{- \infty}^{\infty} f \left( x \right) \operatorname{sinc} \left( \frac{ x - n T }{T} \right) dx && \text{} \\
& \overset{1}{=} \frac{1}{T} \int_{- \infty}^{\infty} \left( \frac{1}{2 \pi} \int_{-\infty}^{\infty} F \left( \omega \right) {e}^{j \omega x} d\omega \right) \operatorname{sinc} \left( \frac{ x - n T }{T} \right) dx && \text{} \\
& \overset{2}{=} \frac{1}{2 \pi T} \int_{- \infty}^{\infty} \left( \int_{- \infty}^{\infty} \operatorname{sinc} \left( \frac{ x - n T }{T} \right) {e}^{j \omega x} dx \right) F \left( \omega \right) d\omega && \text{} \\
& \overset{3}{=} \frac{1}{2 \pi T} \int_{- \infty}^{\infty} T \, \Pi \left( \frac{-\omega}{ {W}_{s} } \right) {e}^{j \omega n T} {F} \left( \omega \right) d\omega && \text{} \\
& \overset{4}{=} \frac{1}{2 \pi} \int_{ - \frac{ {W}_{s} }{2} }^{ \frac{ {W}_{s} }{2} } F \left( \omega \right) {e}^{j \omega n T} d\omega && \text{} \\
& \overset{5}{=} f \left( n T \right)
\end{aligned} $$
Where:

Since $ \mathscr{F} \left\{ f \left( x \right) \right\} = F \left( \omega \right) $.
Changing order of integration for converging integrals.
Applying $ \mathscr{F} \left\{ \operatorname{sinc} \left( \frac{ x - m T }{T} \right) \right\} \left( -\omega \right) $.
Integration boundaries according to the Rect function.
Applying Inverse Fourier Transform at $ x = n T $.

Wrapping it yields:
$$ f \left( x \right) = \sum_{n = -\infty}^{n} \langle f \left( x \right), {g}_{n} \left( x \right) \rangle {g}_{n} \left( x \right) = \sum_{n = - \infty}^{\infty} f \left( n T \right) \operatorname{sinc} \left( \frac{ x - n T }{T} \right) $$
Which is known as the Whittaker Shannon Interpolation Formula.
Conclusion
In the process above the analysis and synthesis of Band Limited Functions is shown using Orthonormal Basis. If one set $ T $ to be the Sampling Interval, usually denoted by $ {T}_{s} $ then the Sampling Frequency is given by $ {F}_{s} = \frac{1}{ {T}_{s} } $.
Now, since the function in frequency domain stretches in the range $ \left[ -\frac{{F}_{s}}{2}, \frac{{F}_{s}}{2} \right] $ we indeed have the known relation that the sampling frequency has to be twice the one sided support of the function in frequency.


Can Carrier Offset cause Image Problems The main reason why we use the Direct conversion receiver with zero IF is that we don't have to deal with the problem of images. For achieving zero IF we need to ensure that the receiver Local Oscillator (LO) is working at the same frequency as the carrier. But I was wondering if there is an offset between the carrier and the receiver oscillator (due to reference mismatch or Doppler shift etc.) do we still have the image problem?
For instance, if the offset is equal to the channel spacing $f_{off} = f_{ch}$, then both the $f_c \pm f_{ch}$ fall at $f_{off}$ after  downconversion mixer. Basically, the offset acting as an IF.
If so, do we still need image reject filters for direct conversion architecture and what then is the benefit of direct conversion architecture.

What I mean is explained in the figure below. We have a carrier frequency $f_c$ but due to some offset we send our data at $f_c+\Delta f_c$. Some other user is sending his data at $f_c - \Delta f_c$ (the image sideband). But after demodulation both the desired and undesired signal components fall back at $\Delta f$

Having images of a signal is very different from having a frequency offset in the signal. Having an image means the signal spectrum is replicated at two or more places.Having a frequency offset in a signal means the spectrum is just shifted in frequency.

First let us look at it from the point of view of the innocent receiver who does not know that such a catastrophe has occurred that the offset itself is equal to the channel frequency. Simply speaking the systems would break down with IQ phase reversed.

If however we have realistic frequency offsets like that of Doppler shifts which for typical cellular frequencies in current wireless is of the order of few hundred Hz then this is usually measured at the receiver and the clock is adjusted. This offset is usually measured by correlating time domain samples after ADC in baseband, since a linear shift in frequency domain equals phase shift in time domain.


When channels are in a cascade and are independent is the capacity the sum? Just wondering if two channels are independent then is the overall capacity, the sum of the individual capacities?

For example suppose I have two BEC where the probability of erasure is $\varepsilon$ and $\varphi$ respectively.

The individual capacities are then  $1- \varepsilon$  and $1-\varphi$ respectively. Now if these channels were independent and in a cascade would the overall capacity be
$1- \varepsilon$ +  $1 - \varphi$?

Best regards and thanks for indulging me.

Let the first channel, channel A, have capacity $C_A=1-\varepsilon$ and the second channel, channel B, have capacity $C_B=1-\varphi$. Now put them in cascade, that is, first channel A then channel B. The resulting capacity of the cascaded channel will be $\text{min}(C_A, C_B)$. You can figure this out by thinking about the information flow like water through pipes. For example, if you're able to push $X$ bits through channel A but only $X-10$ bits through channel B, well then the total capacity will be limited by the smallest capacity channel in the cascaded system. Hope this helps!