Where does $\frac{N}{2}$ came from in approximating an N-point DFT? I've came across the author saying that

  ... for a real cosine input having k cycles in the N-point input time sequence, the amplitude response of an N-point DFT bin in terms of the bin index m is approximated by the sinc function $$X(m)=\frac{A_{0}\color{red}{N}}{\color{red}{2}}\cdot\frac{\sin[\pi(k-m)]}{\pi(k-m)}$$ where $A_{0}$ is the peak value of the DFT’s input sinusiod.

I have no idea where $\frac{N}{2}$ came from. I'm wondering if you can give me some hints.

The common formulation of the forward DFT preserves energy (Parseval's theorem).  This means that a longer constant magnitude sinewave input to a DFT, which has proportionally more energy, must be represented by a proportionally larger value in the DFT result, preserving that greater energy.  Thus a factor of N, which scales with the length of the input.

For strictly real input, this energy is divided evenly between the positive frequency bin (k) and the negative frequency bin (N-k) of a DFT result (for a DFT of length N). These two complex result bins are always complex conjugates of each other (thus equal in magnitude) for strictly real input to a DFT.  Thus another factor of 1/2 is included to do this even energy split.

N and 1/2 results in a factor of N/2 which has to appear somewhere in the result of an energy preserving formulation of a forward DFT.

Note that the 1/2 does not appear for bin 0 or for bin N/2 of even length DFTs, as these have no complex conjugate image (or rather they include both of their complex conjugates summed together into a single result bin).

Also note that this Sinc function result is an approximation.  The actual result is closer to the sum of two Dirichlet functions (also called periodic Sinc functions).  See here for more details.


DSP Benchmarks for ARM Cortex-M4 processor Does anyone know where there are any benchhmark timings for the DSP functionality in the ARM CMSIS V2 running on a Cortex-M4 processor.  I am trying to work out whether this processor is powerful enough to meet my requirements, and to compare it to other microcontroller devices with DSP extensions.

edited to add

I will add a little background. I have an existing product that uses an old and relatively expensive DSP. I have been looking at alternative processors including the Renesas MX and ARM M4 both of which have MACC type instructions. There are benchmark figures for the supplied Renesas libraries that are of the form the x-point FFT at Y MHz takes Z ms, The FIR filter with W coefficients takes Vms etc. I was hoping that there were figures somewhere for the ARM library.

Attempting to evaluate the performance of a number of target platforms myself will be extremely expensive in tools (different compiler for each target), development boards and implementation time.

You're best off looking one level up, at datasheets or published benchmarks for actual devices that are available. One notable manufacturer of Cortex-M4 parts is Freescale. The Cortex-M4 is just a processor core design that is licensed by silicon manufacturers as the basis for their microprocessors. Each manufacturer designs their own peripherals and memory architecture and stitches them together with the core design. Therefore, what you might get from a core-level benchmark is the number of cycles required to execute your algorithm implemented in assembly.

What this doesn't capture are implementation-specific issues that can affect performance, such as memory latency, interrupt response latency, etc. These quantities will vary from manufacturer to manufacturer (and in some cases, across devices in a single family). Therefore, you really need to take the entire microprocessor into account if you're looking to do something that is pushing the envelope of what the chip might be capable of.

Even better: you can often find low-cost development kits for ARM processors (examples for Freescale's Cortex-M4 family). With that in hand, you can get a very good feel for how performant your algorithm will be once implemented in-system.


Sampling frequency of modulated signal A very naive question... A signal with bandwidth of 20 MHz is modulated using a carrier of 1 GHz. What will be the sampling frequency?
I was asked this question during an interview. I answered 2 GHz but the interviewer said it will be 40 MHz. Can anyone explain it?
Thanks.

Nyqvist-shannon sampling theorem  only says that your sampling frequency should be greater than twice the bandwidth of the signal and not twice the carrier frequency of the signal. Hence you can modulate your 20 MHz signal at any carrier frequency you need and still get back your original signal by sampling  at greater than 40 MHz. Hence you need not sample at 2GHZ for recovery of signal, you can do away with a sampling frequency greater than 40 MHz itself  For a more clear understanding refer to bandpass sampling.

http://en.wikipedia.org/wiki/Undersampling


Laser diode bitrate datasheet I'm a student studying communications engineering.  I'm trying to understand why this particular laser diode is made for 155 Mbps and 622 Mbps.  I'm thinking that the bitrate is entirely dependent upon the bandwidth of the photodiode.  I can see from the datasheet that $f_c = 2 \text{ GHz}$ and so what I'm thinking is that you could have, a 1 GHz carrier with around 2 GHz of data.  So wouldn't this allow for a 2 Gbps bitrate?  I'm just wondering if this thinking is correct.

The datasheet can be found here:
Datasheet

I guess I'll convert my comments above into an answer: apart from the theoretical signal processing concepts, physical devices are made to be used in real, physical systems. 155 Mbps and 622 Mbps are common data rates used on fiber-optic SONET networks, specifically OC-3 and OC-12 lines. Therefore, the diode that you referenced is specifically designed and manufactured for use on those systems; it would not be common for such a device to be used for some other nonstandard application (not many people would want to go to the expense of designing their own custom fiber-optic networking standard).

With that said, in theory you could potentially use the devices for higher-bandwidth transmission, barring any specific device-related limitations (I haven't read the full datasheet). As long as the carrier frequency is high enough to accommodate the bandwidth of the signal, then it would be possible (e.g. a 1 GHz carrier can't have any larger than 1 GHz one-sided bandwidth).


Absolute value based AM envelope detection viewed in the frequency domain In my understanding when a passband AM signal is demodulated with an envelope detector taking the absolute value restores the original baseband version of the signal and the low-pass filter removes all the (passband) high frequency content.

It is intuitively obvious when represented graphically that after we rectify the signal we get the envelope in the upper part of the graph.

But mathematically in the frequency domain it is not so obvious how taking the absolute value of an AM passband signal would recreate the baseband signal.

How does one derive mathematically that taking the absolute value of an AM pass-band signal would recreate the modulating signal at baseband?

$\DeclareMathOperator{\sgn}{sgn}$
The modulating signal in AM is

$$s(t) = C + a(t)\text,$$

where $a(t)$ is the (audio) amplitude, and $C$ is a constant so that $s(t) \ge 0 \;\forall t$. (Otherwise, your audio amplitude would just frequently "switch" the wave's sign, not really modulate the envelope.) That means, $C > - \min_t(s(t))$.

Therefore, the passband (around frequency $f$) signal is

$$r(t) = s(t) \sin(2\pi f t) = \left(C + a(t)\right) \sin(2\pi f t)\text.$$

As said, $C+a(t)$ is never negative; therefore the sign of $r(t)$ is

\begin{align}
\sgn(r(t)) &=\sgn(s(t) \sin(2\pi f t))\\
&=\sgn(s(t))\sgn(\sin(2\pi ft))\\
&=+1\cdot\sgn(\sin(2\pi ft))\\
&=\sgn(\sin(2\pi ft))
\end{align}

That tells us one way to write the absolute value of $r(t)$, because for any $x$:

\begin{align}
|x| &= x\cdot \sgn(x)\\
\implies\\
|r(t)|&=r(t)\cdot \sgn(r(t))\\
&= r(t)\cdot \underbrace{\sgn(\sin(2\pi ft))}_{:=w(t)}
\end{align}

This introduces us to a new function $w(t)$, which is always positive one when the carrier wave is, and negative when the carrier is negative. A square wave with the same frequency as the carrier! Thus, knowing that the square wave has a Fourier series representation $w_f(t)= \frac{4}{\pi} \sum\limits_{k\in\{1,3,5,\ldots\}} \frac{\sin\left(2\pi k ft\right)}k $

\begin{align}
|x| &= r(t)\cdot w(t) \\
&= r(t) \cdot \frac{4}{\pi} \sum\limits_{k\in\{1,3,5,\ldots\}} \frac{\sin\left(2\pi k ft\right)}k\\
&= \frac{4}{\pi} s(t) \sin(2\pi f t) \sum\limits_{k\in\{1,3,5,\ldots\}} \frac{\sin\left(2\pi k ft\right)}k\\
&= \frac4\pi s(t){\left[
  \frac11\underbrace{\sin(2\pi f t)\sin(2\pi \,1f\, t)}_
      {=\frac12(\cos(2\pi (f-f)t)-\cos(2\pi(f+f)t))}
+ \frac13\underbrace{\sin(2\pi f t) \sin(2\pi \,3f\, t)}_
      {=\frac12(\cos(2\pi (f-3f)t)-\cos(2\pi(f+3f)t))}
+ \frac15\underbrace{\sin(2\pi f t) \sin(2\pi \,5f\, t)}_
      {=\frac12(\cos(2\pi (f-5f)t)-\cos(2\pi(f+5f)t))} + \ldots
\right]}\\
&= \frac4\pi s(t)\left[{
      {\frac12(1-\cos(2\pi2ft))}
     +{\frac13\frac12(\cos(2\pi2ft)-\cos(2\pi4ft))}
     +{\frac15\frac12(\cos(2\pi4ft)-\cos(2\pi6ft))} + \ldots}
\right]\\
&=\frac2\pi s(t)\left[{
      {1-\cos(2\pi2ft)}
     +\frac{\cos(2\pi2ft)-\cos(2\pi4ft)}3
     +\frac{\cos(2\pi4ft)-\cos(2\pi6ft)}5 + \ldots}
\right]\\
&= \frac2\pi s(t) - \frac2\pi s(t)\left[{\frac23\cos(2\pi2ft)+\frac25{\cos(2\pi4ft)}+\frac27{\cos(2\pi6ft)} + \ldots}
\right]\\
\end{align}

That looks familiar: the first term in the square bracket, $\sin(2\pi f t)\sin(2\pi \,1f\, t)$ is simply mixing the carrier with the carrier, the second is mixing the carrier with a sine thrice its frequency and so on.

Thus, the first term yields, due to $\sin(a)\sin(b) =\frac12(\cos(a-b)-\cos(a+b))$, the baseband signal.

The rest of the terms mix up the passband to $2f$, $4f$, $6f$ and so on and are suppressed with a low-pass filter.


Describing a certain quality of an image Generally I guess I’m asking for a source of image processing terminology.

Specifically here I’m trying to find the word for the quality that would identify the main difference between the Mondriaan on the left and the Monet on the right. Graphical/strong edges vs no solid areas and no well-defined edges.

I know a good amount about image processing but this vocab part is a glaring omission for me.

The Monet image appears to have a much higher entropy.  Whereas the Mondrian image looks like it could very likely be compressed losslessly into a lower number of bits (depending on the compression method) due to a lower informational entropy.


What is the difference between (Sample Timing Offset), (Carrier Frequency offset), and (Sampling Frequency Offset)? Briefly, can anyone explain to me what is the difference between (Sample Timing Offset), (Carrier Frequency offset), and (Sampling Frequency Offset)?

and is Sampling Frequency Offset (SFO) like STO or like CFO?
I found that some papers say SFO is STO, while others say SFO is like CFO so I am a bit confused!
I would appreciate it if someone explains CFO, STO, SFO
thanks in advance

I would have the easiest time explaining these briefly with the aid of an eye diagram as depicted in the graphic below.

In this example we see the constellation pattern of a raised cosine QPSK waveform on the complex plane on the left with I as the real axis and Q as the imaginary axis, and the resulting eye diagram pattern of the resulting real (I) and imaginary (Q) waveforms each on the right.  The eye diagram is what we would see if we triggered an oscilloscope to the symbol clock with infinite persistence, and here we see a 3 symbols duration superimposed over many symbols.

Notice the red dots on the left-- these are samples with no Sample Timing Offset, corresponding to the time locations on the right where the vertical variation between all the samples is minimized-- the ideal time to sample the waveform in order to make the final decisions on what data was transmitted. As we shift that sample clock to the right or left (Sample Time Offset), the vertical variation will increase, reducing our margin against noise to make correct decisions.

This constellation as shown is with no Carrier Frequency Offset (The waveform is centered at DC, the baseband signal). If there was a Carrier Frequency Offset, the entire constellation would rotate at the rate of that Offset-- In fact we measure the rate of rotation of the red dots to determine and correct for carrier offset error.

Finally Sampling Frequency Offset would be if the sampling clock itself had an error in frequency. This will change the effective Sample Time Offset and Carrier Frequency Offset proportional to the sampling rate error. If the sampling rate error was 1ppm, then the Sampling Frequency Offset and Carrier Frequency Offset as observed would also be in error by 1ppm. (For example, a 10MSps sampling clock with a 1 ppm error would have a 10 Hz error. If the observed carrier offset was 1 Hz with the 10MSps clock with no error, the carrier offset as observed would have an additional error of 1e-6 Hz.)

Note that in OFDM, a Sampling Time Offset can produce results similar to Sampling Frequency Offset in that the individual constellations in each sub-carrier will rotate from sub-carrier to sub-carrier (so a rotation in frequency rather than a rotation in time). This is because a time offset is a linear phase in the frequency domain. So pay very close attention to what domain you are in: If you see a rotation in the frequency domain, this is due to an offset in time. If you see a rotation in the time domain (on a single sub-carrier from sample to sample in time) then this is due to an offset in frequency.


Spectrogram without using stft Can I use any kind of transform to make spectrogram, not just stft? Is it possible to get the spectrograms of a signal without using stft?
What difference does it make if I calculate the frequency components of the entire signal using fft instead calculating fft of few samples at a time (stft) and cover the entire signal by moving the window?

Can I use any kind of transform to make spectrogram?

If your input signal is a finite length sequence and you want to analyze how the frequency domain components are varying in that input sequence, you have to use one of the techniques used for time-frequency analysis. STFT is one such technique, there are others also like wavelet transform, quadratic time-frequency distribution functions, etc.
For the $2^{nd}$ question, if you take FFT of the entire signal as 1 window, then those FFT coefficients will show all of the existing frequency components at once, you won't have the time-information. Dividing the input into smaller chunks of overlapping windows and then taking FFT of each window provides you with time information too, as in what frequency component was significant at what time, during the whole signal collection.


how to get the location long and lat out of the raw measurements? I'am making a research about GNSS
I used GNSS logger to collect measurements
but couldn't understand how these measurements turn into location of long and lat?
the picture below of the measurements attributes.

I only know the vastly simplified form: all satellite nav systems work by the satellites telling you where they are and what time it is.  In GPS terms, that time is the "pseudorange", and it's all you need (in theory) to locate yourself.

With no satellites in view, you have no clue where or when you are.  With one satellite in view, you can place yourself on a 3-sphere* in time and space**.  With two satellites you can place yourself on a 2-sphere, etc.  When you have four satellites in view you can locate yourself to a point, and you're done***.

I suggest you look around and get a book on satellite navigation.  I'd love to recommend one, but first, I'm not current, and second, StackExchange discourages product recommendations.  But they're out there, and these days there should be one that covers GPS, GNSS and Galileo.

* I'm using the mathematical terminology here: a 3-sphere is the locus of points equidistant from a point in 4-dimensional space, a 2-sphere is the surface of a ball in 3-D space, a circle on a plane is a 1-sphere, and a point is a 0-sphere.

**  I.e., you know that if it's the same time as the satellite's time you're at the satellite's position, if it's 1ms of a second after the satellite's time you're on a 2-sphere 300,000 meters in diameter around the satellite's position, etc.

*** Except for noise, of course, and I can't recall if, even in the perfect case, you end up with just one possible location or a finite number of them.


How to Use Maximum a Posteriori Probability (MAP) in Classification Task I have a 2D image defined on a region $\Omega$. Let $I: \Omega \to R$ be a gray image. Assume that the region can be separated into $N$ sub-regions $\Omega_i$ such that $$\forall i,j=1,\ldots ,N:\Omega_i \cap\Omega_j=\emptyset$$ and $$\bigcup_{n=1}^{N}\Omega_n =\Omega$$
To simplify, we assume $N=3$ and it shows in below figure.

We denotes $p(x \in \Omega_i|I(x))$ is posteriori probability of region $\Omega_i$ given $I(x)$. Based on the Bayesian rules, we have
$$p\left(x \in \Omega_i|I(x)\right)=\frac{p\left(I(x)|x\in \Omega_i\right)p(x\in \Omega_i)}{p(I(x))}$$
The segmentation try to assign pixel $x$ to each region based on posteriori probability. That goal is that finds maximum a posteriori probability (MAP) of pixel $x$ in each region. It can express by formula
$$\left\{ \Omega_{i,x} \right\}_{i=1}^{N}=\arg \max_{\Omega_i} \left\{\prod_{i=1}^{N} p(I(x)|x \in \Omega_i)p(x \in \Omega_i )\right\} \tag{*}$$
Assuming that the pixels within each region are independent , the MAP will be achieved over all image domain $\Omega$
$$\left\{ \Omega_i \right\}_{i=1}^{N}=\arg \max_{\Omega_i} \left\{\prod_{i=1}^{N}\prod_{x \in \Omega_i} p(I(x)|x \in \Omega_i)p(x \in \Omega_i )\right\} \tag{**}$$

Are my formulas $(*)$ and $(**)$ correct?
If not, Could you help me to correct them?

The reference is Geodesic Active Regions and Level Set Methods for Supervised Texture Segmentation.

I will try to give you some intuition into it by a different example.

Think we have 3 machines which can generate the numbers 1, 2, 3.
The first machine generates the number 1 with 80% and the numbers 2, 3 with 10% each.
The second machine generates the number 2 with 80% and the numbers 1, 3 with 10% each.
The third machine generates the number 3 with 80% and the numbers 1, 2 with 10% each.

One of the machines is chosen, you don't know which and the generated number is 2.
What machine would you bet it would?

The question above is the likelihood with no prior knowledge.
Hence, given the number 2 the most likely is the second machine.

Yet,
What happens if you are told that the chosen machine would be chosen following this rules - 98% the first machine, 1% the second and 1% the third.

Now, what are the chances?

Seeing 2 from machine 1 are 0.1 * 0.98 = 0.098.
Seeing 2 from machine 2 are 0.8 * 0.01 = 0.008.
Seeing 2 from machine 3 are 0.1 * 0.01 = 0.001.

Now a better choice would be machine 1.
This is the MAP, because we took into account both the likelihood and the prior knowledge.


Get difference equation of a filter from MATLAB transfer function coefficients I need to get the difference equation of a specific elliptic filter.
I calculated the transfer function coefficients in MATLAB with:

  %% Low pass design
  n = 10;
  passband_ripple = 1;
  stopband_attenuation = 60;
  f = 4750;
  fs = 44100;
  [B,A] = ellip(n,passband_ripple,stopband_attenuation,f/(0.5*fs),'low');

However, when I use the difference equation, I calculated from the transfer function, the samples go to +-inf. Here is the code where I create the noise input sound and filter it with difference eq.

  xin = rand(1,fs/20);
  xin = xin - mean(xin);
  xin = [xin zeros(1,fs*2)];

  x = xin;
  y = zeros(1, length(x));

  %% Low pass
  b1 = 0.0027;
  b2 = -0.0134;
  b3 = 0.0362;
  b4 = -0.0651;
  b5 = 0.089;
  b6 = -0.0978;

  a1 = -7.953;
  a2 = 29.7485;
  a3 = -68.5714;
  a4 = 107.5494;
  a5 = -119.7162;
  a6 = 95.6848;
  a7 = -54.2043;
  a8 = 20.8340;
  a9 = -4.9109;
  a10 = 0.5401;

  for n = 11:length(y)
      y(n) = b1*x(n) - b2*x(n-1) + b3*x(n-2) - b4*x(n-3) + b5*x(n-4) - ...
      b6*x(n-5) + b5*x(n-6) - b4*x(n-7) + b3*x(n-8)- b2*x(n-9) + ...
      b1*x(n-10) - ...
      a1*y(n-1) - a2*y(n-2) - a3*y(n-3) - a4*y(n-4) - a5*y(n-5) - ...
      a6*y(n-6) - a7*y(n-7) - a8*y(n-8) - a9*y(n-9) - a10*y(n-10);
  end

The filter function does what I want, but I need to be able to filter the sound sample-by-sample.

y = filter(B,A,x);

How can I obtain the right difference equation of my filter?

There are two problems with your code. First, you should not invert every other $b$ coefficient, they should all have a positive sign, just like the $a$ coefficients that all have a negative sign. Second, - and that's the main problem -, why do you quantize the coefficients? Your actual coefficient $b[1]$ is 0.00266475348894505 (and even longer), and you rounded it to 0.0027. That's bad for the $b$ (numerator) coefficients, but detrimental for the $a$ (denominator) coefficients, because the filter's characteristic is significantly changed, and the filter can become unstable. You should just use the coefficients as stored in the vectors A and B.

Implementing the difference equation is a nice exercise to understand what's going on, but I think you'd be better off using filter.m.


Maximum phase with the same magnitude response I'm trying to solve an exercise in a general way and I can't find if my answer is correct.

Given this frequency response:

$$H(z) = 4 + 2\sqrt2 z^{-1} + z^{-2}$$

I need to find a frequency response that has maximum-phase response $H_{1}(z)$ such that $\lvert H(z)\rvert = \lvert H_{1}(z)\rvert$.

First of all, I write:

$$H(z) = \frac{4z^2+2\sqrt2 z + 1}{z^2}$$

So I see that I have zeros in:

\begin{align}
z_{01} &= \frac{\sqrt2}{4}+j\frac{\sqrt2}{4}\\
z_{02} &= \frac{\sqrt2}{4}-j\frac{\sqrt2}{4}
\end{align}

And poles in:

$$z_x = 0 \text{ (double)}$$

Now, I know that for a maximum phase system I need that all the poles and zeros have to be outside the unit circle, but from Oppenheim-Schafer I also know that I have to add a scalar factor for the minimum-phase case, so I suppose I also have to add that scalar factor now.

So I say that this:
$$H_1(z) = k \left(z- \frac{1}{z_{01}}\right) \left(z- \frac{1}{z_{02}}\right)$$
(I used that I have both poles in $z=\infty$, is that right?)
has to be the one that makes $\lvert H(z)\rvert = \lvert H_{1}(z)\rvert$.

But I can't find the value of $k$ (scalar factor), because when I write $z=e^{j\omega}$ and use the magnitude definition (magnitude = square root of the squared real part plus the squared imaginary part) I find that $k=1$ and $k=-1$.

Note that the zeros of the given transfer function have a negative real part:

$$z_{01}=\frac{1}{2\sqrt{2}}(-1 + j),\quad z_{02}=z_{01}^*\tag{1}$$

With $(1)$ you can write the transfer function as

$$H(z)=4(1-z_{01}z^{-1})(1-z_{01}^*z^{-1})\tag{2}$$

Reflecting a zero across the unit circle means replacing it by its conjugated inverse, i.e. $z_{01}\rightarrow 1/z_{01}^*$, from which you get

$$H_{max}(z)=k\left(z-\frac{1}{z_{01}}\right)\left(z-\frac{1}{z_{01}^*}\right)=\frac{k}{|z_{01}|^2}(1-z_{01}z)(1-z_{01}^*z)\tag{3}$$

The constant $k$ can be determined by evaluating $(2)$ and $(3)$ at a point on the unit circle $|z|=1$. If you choose $z=1$ you get

$$H(1)=4(1-z_{01})(1-z_{01}^*)\tag{4}$$

and

$$H_{max}(1)=\frac{k}{|z_{01}|^2}(1-z_{01})(1-z_{01}^*)\tag{5}$$

Comparing $(4)$ and $(5)$ you get

$$\frac{k}{|z_{01}|^2}=4\tag{6}$$

and with $(3)$

$$H_{max}(z)=4(1-z_{01}z)(1-z_{01}^*z)\tag{7}$$

Comparing $(7)$ with $(2)$ you can see that

$$H_{max}(z)=H(z^{-1})=z^2+2\sqrt{2}z+4\tag{8}$$

In the time domain this corresponds to a time reversal. If $h[n]$ is the impulse response corresponding to $H(z)$, and $h_{max}[n]$ corresponds to $H_{max}(z)$, then the following holds:

$$h_{max}[n]=h[-n]\tag{9}$$


Mathematical relationship between highpass and lowpass filtering Let $g, h_{HP}, h_{LP}: \mathbb{R} \rightarrow \mathbb{R}$ and $G, H_{HP}, H_{LP}$ denote their continuous Fourier transforms under the Fourier operator $\mathcal{F}$. Let $*$ denote the continuous convolution of two functions and $\cdot$ the pointwise multiplication.

Define the shorthand $1: \mathbb{R} \rightarrow \mathbb{R}, x \mapsto 1$ and let $H_{HP} = 1 - H_{LP}$ in the frequency domain, then I would like to show
$$
g * h_{HP} = g - (g * h_{LP}).
$$

In words:

  Filtering an image with a highpass is the same as subtracting the lowpass-filtered image, where the lowpass is the inverted highpass.

There is a proof (from some closed-source course notes) which goes as follows:

\begin{align*}
& g * h_{HP}\\
& = \mathcal{F}^{-1}(G \cdot H_{HP})\\
& = \mathcal{F}^{-1}(G \cdot (1 - H_{LP}))\\
& = g * \mathcal{F}^{-1}(1 - H_{LP}) \quad \text{(The next step is unclear to me)}\\
& = g * (1 - h_{LP})\\
& = g - (g * h_{LP})
\end{align*}

In fact, $\mathcal{F}^{-1}(1 - H_{LP}) = 1 - h_{LP}$ is equivalent to $1 = \mathcal{F}(1)$ by applying $\mathcal{F}$ on both sides, which is a contradiction with the first identity mentioned below.

Questions

  Is it true that the proof above lacks details?
  Is my proof below correct/more rigorous?

It is known that

$\mathcal{F}(1)=\delta$, from Constant function on MathWorld
$\mathcal{F}^{-1}(1)=\delta$, from Math Stack Exchange: How to prove that inverse Fourier transform of “1” is delta funstion?

I found another approach to possibly prove the proposition above:

First note: If $u, v: \mathbb{R} \rightarrow \mathbb{R}$ are continuous functions and $\int_a^b u(x) dx = \int_a^b v(x) dx$ for every $a, b \in \mathbb{R}$, then $u(x) = v(x)$ for every $x$, i.e. $u=v$.

\begin{align}
&\int_{a}^{b} (g * h_{HP})(x) dx\\
& = \int_{a}^{b} (g * \mathcal{F}^{-1}(1 - H_{LP}))(x) dx \quad \text{(as stated above)}\\
&= \int_{a}^{b} (g * (\delta - h_{LP}))(x) dx \quad \text{(linearity + identity above)}\\
& = \int_{a}^{b} \int_{-\infty}^{\infty} g(\tau) ((\delta - h_{LP})(x - \tau)) d\tau dx\\
& = \int_{a}^{b} \int_{-\infty}^{\infty} g(\tau) (\delta(x - \tau) - h_{LP}(x - \tau)) d\tau dx\\
& = \int_{a}^{b} \int_{-\infty}^{\infty} g(\tau)\delta(x - \tau) - g(\tau)h_{LP}(x - \tau) d\tau dx\\
& = \int_{a}^{b} g(x) - \int_{-\infty}^{\infty} g(\tau)h_{LP}(x - \tau) d\tau dx\\
& = \int_{a}^{b} g(x) - (g * h_{LP})(x) dx\\
& = \int_{a}^{b} (g - (g * h_{LP}))(x) dx
\end{align}

Especially, I am unsure whether the usage of the $\delta$ distribution is rigorous here.

Yes you are right. The following line (excerpted from your post) is wrong, below which I place the corrected version:
\begin{align*}
& g * h_{HP}\\
& = \mathcal{F}^{-1}(G \cdot H_{HP})\\
& = \mathcal{F}^{-1}(G \cdot (1 - H_{LP}))\\
& = g * \mathcal{F}^{-1}(1 - H_{LP}) &\quad \text{(The next step is unclear to me)}\\
& = g * (1 - h_{LP}) &\quad \text{(This line is wrong)}\\
& = g * (\delta(t) - h_{LP}) &\quad \text{(This is the corrected version)}\\
& = g - (g * h_{LP})
\end{align*}

I don't think so if you really need to go deeper than the above corrected proof to show that the inital claim was right. But it's up to you of course.


Question About Linear and Circular Convolution - 1D and 2D We know that applying a filter on an image is called correlation or convolution depending on the filter angle. In Gonzalez I have read that we can apply linear convolution on an image. Here I have a doubt.

Suppose we have a vector image of 4 elements and a filter/kernel also of 4 elements then we get the output of linear convolution which contains 7 elements (4+4-1). Does that mean now the image is represented by 7 elements? If not then what does these 7 elements represent?
Also I want to know where do we use Circular Convolution. I have read many posts but not able to get it fully.
Last question is from filtering in frequency domain. We know that convolution of f(x,y)* h(x,y) in the spatial domain is equal to the product of DFT's of the image and the kernel i.e F(u,v).H(u,v) in the frequency domain. Then if we have a 64 *64 image and 3*3 filter, please tell me what will be the size of the sequence obtained by applying Linear Convolution in the spatial domain. Also Do we take the DFT of the image f(x,y) as it is or do we need to pad the image first before taking its DFT and if padding is needed how much? I know filter size should also be 64 * 64, so we need to change filter size from 3 * 3 to 64 * 64. It means padding of 61 rows and 61 columns of 0's but where should we pad these 0's. These 0's cannot be equally divided top,bottom,left , right as 61 * 61 is an odd number. Please explain.

Let me answer you:

For a signal of size $ m $ and a filter of size $ n $ the output of Linear Convolution is $ n + m - 1 $. In case of 2D signal of size $ \left( m, n \right) $ and filter of size $ \left( p, q \right) $ the output size is $ \left( m + p - 1, n + q - 1 \right) $.
You can read about Circular Convolution in Wikipedia. Basically when a convolution is applied on finite discrete signals one should take care of the boundaries. In most cases the default is assuming the signal i padded with zeros which results in Linear Convolution. If you use padding which build a periodic / circular signal and then apply convolution you will get Circular Convolution. It turns out that frequency domain multiplication of discrete signals is equivelnt of Circular Convolution in spatial domain.
You need to pad it with zeros and line the axis origin to match the image.
Have a look at my answer for Kernel Convolution in Frequency Domain - Cyclic Padding. I also shared a MALAB code which shows how to build the kernel correctly.


Why are anticausal systems even defined? I guess the same question is usually asked for complex numbers too, but the fact is that complex numbers are used all the time practically. However, at least on a quick google search, I couldn't find any applications of anticausal systems. So why define it in the first place?

Of course they don't exist. But we can stop time and use systems that would be non-causal if we hadn't stopped time. Stop time? Yes, just store your data and work offline / non-realtime. Or work on data that are not temporal but something else, for instance spatial. Of course, non-causal is a misnomer in that case, but the term is still used, in analogy with temporal signals and systems.

One common "non-causal" system that is used in offline mode is filtering followed by reversed filtering using IIR filters. This is done in order to eliminate phase distortions introduced by the conventional (causal/realtime) use of the filter. See this answer for the details.


What is the difference between a range image and a depth map? https://en.wikipedia.org/wiki/Depth_map

https://en.wikipedia.org/wiki/Range_imaging

To me, it sounds like they are the same thing, except depth map implies only having the z coordinate info.

But then again, if you have a range image vs depth image, is this something one can tell the difference by examining the data or is the term mostly contextual?

Your understanding is correct, it's z-values vs distances.

A depth map contains the distances of points in the scene to a plane.

A range image contains the distances of points in the scene to another point. This point is very often the camera's position.

For a depth map the plane that the distances are computed with respect to can be any plane, but it is often parallel to the image plane. Examples of common planes:

The image plane
The focal plane
A plane closer to the object, such that the points closest to the camera have distance zero.

If you of the same scene have:

A range image with the distances to the camera
A depth map with  the distance to the image plane

these two images will look very similar, and it's not possible to conclude whether it's a range image or depth map, unless you have something you know is a plane in the image.


Cascading first order digital filters in C++ In a related question a probable solution was given to build a first-order digital filter and then cascade three of them in order to turn white noise into pink.  I have applied the C++ as follows but still the signal sounds and looks like white noise.
I would like to know what is wrong with my implementation and why I don't hear pink noise.  For reference, the poles and zeros come from Robert Bristow-Johnson's work here.
Header:
float *state = nullptr;
Implementation:
state = new float[0.0]; in Constructor.
Then in the loop, for (int i=0; i < numSamples; i++)
float first = first_order_filter(whiteNoise, 0.99572754, 0.98443604, state);
float second = first_order_filter(first, 0.94790649, 0.83392334, state);
float third = first_order_filter(second, 0.53567505, 0.07568359, state);
out1 = third;

Where first_order_filter is defined as in Robert's answer here:
https://dsp.stackexchange.com/a/70963/11391
I would love to know if this code is approximately correct/ where the problem lies.

My late night brain made foolish mistakes, for the record, if anyone needs this, the working code is as follows:
Header:
float state1;
float state2;
float state3;

Implementation:
In Constructor:
state1 = 0;
state2 = 0;
state3 = 0;

Robert's function:
// this processes one sample

float first_order_filter(float input, float pole, float zero, float *state)
{
    float new_state = input + pole*(*state);
    float output = new_state - zero*(*state);
    *state = new_state;
    return output;
}

Then in the loop, for (int i=0; i < numSamples; i++)
float first = first_order_filter(whiteNoise, 0.99572754, 0.98443604, &state1);
float second = first_order_filter(first, 0.94790649, 0.83392334, &state2);
float third = first_order_filter(second, 0.53567505, 0.07568359, &state3);
out1 = third;

Of course this works like a charm, producing Pink Noise and if the poles and zeros are swapped it produces Blue/ Azure Noise.


Error probability for 8-PSK How can I calculate bit error rate for 8-PSK modulation when I have $E_b/N_0$ given? I know that for QPSK/BPSK the matter is quite simple, but I cannot cope with 8-PSK.
Let's say that my $E_b/N_0$ = 63.1 [dB].
Now, can I use this formula?:

If yes, what can I substitute for n and Q? And if not, how can I calculate it differently?

First of all, you're understating your problem. The error probability depends not only on your noise power (density), but also on the actual distribution of noise. Also, you need to realize that bit error probability is not the same as symbol error probability, not even proportional, so at best, the formula you pasted (from an uncited source) is only an approximation for the special case of symbols being Gray-coded, and even then only for noise powers that are relatively small. I find it highly unlikely a general formula (for variable $M$, which you also didn't define, but I assume is the $M$ in M-PSK) is that simple.
I don't know the correct formula by heart, either, so me looking it up and check that my scenario makes exactly the same assumptions as the author of the thing I'm reading is as good as you doing it, with the difference that you actually know your assumptions, whereas I can only vaguely infer from your question statement what they might be.
One of the assumptions that I have infer is that you're considering additive, rotationally symmetric complex Gaussian noise to equiprobable symbols. And that you're assuming Gray coding. And that you're using a Maximum Likelihood detector.
It's likely you'll have to use the $Q$-function within your derivation at some point, but it's not as trivial as in the BPSK or QPSK case, since the decision boundaries aren't orthogonal to each other, so looking at the noise in I and Q separately simply doesn't work!
This is one of the classical problems with which stochastics tortures students (hint: that means there's a bazillion of resources out there). I'll outline the solution:

As in your previous question, you have to define a noise model before you start calculating anything related to noise. No shortcut there. Write it down.
Then, you calculate the probability density of your received value $R$ given that you've transmitted some specific symbol $s_T$ from the symbol alphabet $\mathcal S$: $f_{R|S_T}(r|s)$. In the smooth-pdf additive noise case, this boils down to the noise density $f_N$ shifted by the transmit symbol: $f_{R|S_T}(r|s)=f_N(r-s)$. Be careful, $R$ and $S_T$ are complex random variables w.l.o.g.
Then, you need a mathematical model of your detector, which maps received values $r$ to symbols $S_s\in\mathcal S$, is some function $S_r = d(s): \mathbb C \mapsto \mathcal S$.
Now you apply that model to your received symbols; this is often a coordinate transform $g(\cdot)$. For example, in PSK, you need to convert the received complex value from representation as real and imaginary part to polar coordinates. You will get a transformed random variable $\tilde R$, whose PDF you need to derive from $f_{R|S_T}$ and $g(R)$. If your stochastics for engineers course was anything like mine, you'll remember it here: Inverse of the mapping, Jacobi determinant, marginal density, done.
Now, based on that, it should be possible to derive the probability that the transmitted symbol is detected as each of the possible symbols (including itself, i.e. the correct one). Summing over the probability of each symbol multiplied by the bits that differ between that symbol's bit mapping and the transmitted symbol's bits gives you the bit error probability given that specific transmitted symbol. Do that for all transmit symbols. Often, you'll notice that the symmetry of the decision boundaries and the noise properties allow you to only calculate one (like in PSK in AGN) or a few (like in QAM constellations in AGN) probabilities instead of all.
Weigh each transmit symbol bit error probability with the transmit symbol's probability, sum, and get the overall bit error probability.

Notes

I'm sorry to pester you with 1., but you need to really really be more stringent about writing down models. I can't look into your head and tell myself what noise model you're assuming!
Point 2 is a dangerous game to play for many people: they infer the Maximum Likelihood decider is always based on the symbol with minimum distance to the received complex number. That's only the case for symmetric, concave noise PDFs.
If you didn't have my "stochastics for engineers" course, "random variable transformation PDF" is what you want to look up
Gray coding comes into play at 5., where a common simplification / approximation is to assume that the symbol errors that bring you into the neighboring wrong symbols are so much more likely than the symbol errors that take you multiple symbols away from what was transmitted, that you can neglect them altogether. Then, you can only consider the probability of the neighbors, which, due to Gray coding, have only a single bit error. However that's a strong assumption, especially if your noise isn't Gaussian or your SNR isn't good. Really note, whenever you're doing that simplification, that you're doing it, and why you think it's feasible.
Here's where you'll notice that while in PSK, all symbols have the same bit error probability, that's not the case for other constellations. Thus, you'll find modern transmission systems where the channel encoder's output transmit symbols aren't equally probable by design. (this gets complicated, because not making the symbols in a QAM equally probable changes the average transmit power, so you need to rescale the transmission, but then your E_b/N_0 chances, and you need to re-optimize and then your number crunching cluster goes brrrt and gets warm.)


The reasons for filter explosion Assume signal $x$, sampled at $f_s = 44100\; \mathrm{Hz}$. I tried to filter it using the Butterworth bandpass filter ( $30\; \mathrm{Hz} - 70\; \mathrm{Hz}$) of order $8$. However, as a result I get a vector with most elements being NaN (and some of them extremely small, approx.  $-2.5 \cdot 10^{306}$`).

If I try the same filter of order $6$, I get results as expected. What could be possible reason for order $8$ filter to 'explode'?

Here is the MATLAB code, just in case I made an error which I don't see:

[b, a] = butter(4, [60 / fs, 140/fs]);
x_filtered = filter(b, a, x);

Basically you never want to use the Transfer Function representation (with b and a) and rather use the Zeros-Poles-Gain (z,p,k). This will allow you to avoid the numerical errors. In your case you might design your filter in following way:

fs = 44100;           % Sampling frequency
Wp = [30, 70]/(fs/2); % Pass band frequencies (as normalized frequency)
Ws = [20, 90]/(fs/2); % Stop band frequencies
Rp = 3;               % Ripple at pass band
Rs = 50;              % Ripple at stop band

[n, Wn] = buttord(Wp, Ws, Rp, Rs);     % Get order and omega vector
[z, p, k] = butter(n, Wn, 'bandpass'); % Design filter accordingly
[sos, g] = zp2sos(z, p, k);            % Convert to state matrix
Hd = dfilt.df2sos(sos, g);             % Create the filter object

Which for some dummy random signal:

x = rand(1,100000);
y = filter(Hd, x);

Will produce a stable output:

And here is the filter frequency response (everything looks as requested):


Image processing based on wavelet transform in python I take the dwt2 for an image and saved it's coefficients  (LL,Lh, hl,hh) using pywt.dwt2 (image,'haar'),in my project I have to change them to uint8 but when I change their types and reconstruct the image an error appear prevent the show of the result image .can any body help me with this problem! My code is like:

import cv2
import numpy as np
import pywt
image=cv2.imread('1.tif')
image=cv2.resize(image,(256,256))
[b,g,r]=cv2.split(image)
x=b
coffes=pywt.dwt2(x,'haar')
ca,(ch,cv,cd)=coffes
ca2=np.uint8(ca)
ch2=np.uint8(ch)
cv2=np.uint8(cv)
cd2=np.uint8(cd)
coffes=ca2,(ch2,cv2,cd2)
origin=pywt.idwt2(coffes,'haar')
cv2.imshow('image',origin)
cv2.waitKey(0)
cv2.destroyAllWindows()

Three aspects should be taken into account. Suppose that the image is initially positive, coded on eight bits. Then:

standard wavelet oefficients become real, and are not converted to uint8 easily, because wavelet filter coefficients are often real (and not dyadic rationals)
since the filters are orthogonal, the coefficients risk to grow, and exceed the initial $[0,\ldots,255]$ range. So a mere rounding is dangerous
even if you use integer wavelets (by lifting), the same happen, and in 2D you may need $1$ more bit for each level
some detail coefficients will be negative, so uint8 will destroy the signs in a way you cannot recover afterward.

So, the inverse wavelet tranform is likely to be crappy. THe simplest way around:

divide all coeffciients by their maximal absolute value (and keep track of it),
code all the $[0,1]$ magnitudes on $7$ bit,
code the sign on the $8$th bit.


Filtering an acceleration signal I have a triaxial accelerometer giving me acceleration along x, y, and z. I then combine the 3 axes by calculating magnitude of acceleration over all time points (details in my question here)

I'm wondering about the correct time in the process to low pass filter my signal to eliminate some of the high frequency noise in the signals. I guess I have a few options:

Filter the raw xyz signals, calculate magnitude and use it as is
Filter the raw xyz signals, calculate magnitude, filter the resulting magnitude vector, and then use that
Use the raw xyz signals (no filtering), calculate magnitude, filter the result then use it

From a theoretical standpoint, does it matter which signals I filter and in which order?

I would imagine that if I dont filter the raw signals (#3) and directly calculate the magnitude I would end up amplifying a lot of noise? But is there a harm to double filtering (the raw signals, and then the magnitude) in #2?

I believe the best strategy is to filter prior to calculating the magnitude. To see this easily, consider the low pass filtering as an averaging process and consider the noise as a zero mean Gaussian white noise. The distribution of the noise prior to the magnitude computation as zero-mean Gaussian will average to zero if the noise is white. The distribution of the magnitude of a Gaussian process is Rayleigh which has a maximum at the variance $\sigma$ of the noise distribution with an average of $\sigma \sqrt{\pi/2}$.  With a signal present the distribution is Ricean with similar considerations.

Filtering after will help in removing the variability of your result to the extent the system is stationary but will not reduce the error due to the noise; simple case in point is as above with noise only (and stationary), the result converges to zero; while after the magnitude the result will converge to a constant proportional to the noise level.

I suspect for both filters you will be limited by the dynamics involved in how much you can filter, which is ultimately the trade you will make with desired SNR of your result.

For more details on this specific to filtering before or after computing the magnitude (or equivalently before or after decision) see this great app note by Agilent:

http://my.ece.ucsb.edu/York/Yorklab/Useful%20Stuff/Tutorials/SpecAnalyzerMeas-Noise.pdf

There are applications where it is theoretically better, such as when the phase of the information is not preserved.

Here is a decent reference on both cases:

http://www.acfr.usyd.edu.au/pdfs/training/sensorSystems/10%20Detection%20of%20Signals%20in%20Noise.pdf

If you don't want to read or understand the reference materials in detail, or have little time- I suspect filtering prior to taking the magnitude will be better but you can easily test this by establishing SNR metrics on your result; just pay attention to rate of change as well with regard to how fast your system will be changing and how fast you want to be able to track that versus averaging out those changes.


How to pad a windowed sinc filter in the frequency domain I'm creating windowed sinc filters to apply them to certain signals that I'm dealing with. To design the filters, I'm using the approach described in the book "The Scientist and Engineer's Guide to DSP". Here's a brief resuming:

\[h[i] =
\left\{
\begin{matrix}
Kw(i)\frac{\sin(2\pi f_c(i - \frac{M}{2}))}{i - M/2} & 0\le i \le (M+1)\\
\\
0 & \text{otherwise}
\end{matrix}
\right.
\]

\[
w(i): \text{A window function}\\ K: \text{DC gain}\\ f_c: \text{The cut frequency expressed as a fraction of the signal sampling rate}
\]

After calculating the filter kernel, I obtain the frequency response of the filter by taking the FFT of h[i] using an algorithm that approach the filter length to the nearest power of two(to get better performance). The same procedure it's also applied to the input signal in order to get its frequency response.

To my understanding, the next step would be multiply the two frequency responses to get the filtered signal in the frequency domain, but the problem is that both signals are not the same size and I haven't found an explanation on how to deal with this problem. So people, how should I pad the filter signal to make this point-wise multiplication? Should I fill the end of the filter with zeros to match with the signal size or should I do it symmetrically putting zeros to the left and right of the filter response?

Thanks in advance!

You don't pad it in the frequency domain, you pad it in the time domain (i.e. before you calculate the FFT).  You can put the zeros either symmetrically or after the filter.  Either way works, it just results in a shift in the output.


Math for projecting 3d world point to VR image coordinates Suppose, I've a 3d world point $P$ and I want to project it onto the image plane of camera, I can do it as follows
$$M = \begin{bmatrix}
R & t \\
0 & 1 \\
\end{bmatrix}$$
$$p = K \cdot M \cdot P$$
where $K$ is the camera intrinsic matrix and $p$ gives my image point in homogeneous coordinates.
Now, I want to do the same thing for VR images (360° or 180° degree). Can you provide me any links for the math involved here?
I've attached two sample images here where the camera moves forward. As you can see the change in the frames is not the same when camera moved forward in normal videos. I believe it uses some other coordinate system (may be spherical) and also the images are such that when viewed closely on a headset, it gives 3d feeling.

Your images look like being shot with an angular fisheye projection . Also, to account for epipolar geometry of stereo images, you should include the fundamental matrix into the transform engine of your geometry pipeline.
To render a 3D scene in OpenCV with a fisheye camera, one uses a fisheye camera model. Regardless of a graphics software you are using, you may find useful a description of this model's math in an OpenCV documentation.
EDIT:
When using an application like Blender, you need not construct matrices and call API methods. You just select values for dialog controls. For example, to select a fisheye lens, you do it in accord with the following ref.:
https://docs.blender.org/manual/en/2.79/render/cycles/camera.html :
CAMERA
.   Lens
.   Type : Fisheye
See also https://www.veer.tv/blog/blender-tutorial-how-to-render-a-3d-vr-video-from-blender/


Non Uniform DFT Algorithm by Dutt & Rokhlin I have implemented Dutt & Rokhlin's FFT algorithm [1] for non-equispaced data but for some reason I am getting very large errors ($E_2=0.04$ and $E_\infty=0.1$). I was wondering if there are any freely available implementation of this method so I can compare my implementation.

[1] Dutt, Alok, and Vladimir Rokhlin. "Fast Fourier transforms for nonequispaced data." SIAM Journal on Scientific computing 14.6 (1993): 1368-1393.

You may download the NFFT Library - Non Equispaced Fast Fourier Transform.
It has interfaces for MATLAB, Octave and Julia.


what does the order of filter mean? Recently, I read a paper. In it, it said "a $32^{nd}$ order 100MHz bandpass filter". So, how the filter was designed? Does it mean this filter have 32 taps?

The order is classically the maximum delay (in samples) that is  used to produce each output (Filter order by JOS). For FIR filters, the length is "one plus order":

$$y[n] = x[n]-x[n-1]+x[n-2]$$

would be 3-tap and order 2.
If the filter is IIR, it would have "more coefficients" than "one plus order". I am not sure one can say "taps" for IIR filter. This recursive  equation would be order 2:

$$y[n] = x[n]-x[n-1]+y[n-1]+y[n-2]$$

A $32^{\mathrm{nd}}$-order with infinite impulse response (IIR) is likely to have stability issues. Even if it does not, the technique that goes "with each tap randomly varied by up to 5% to simulate physical differences" would be likely to endangered the stability randomly.


confidence interval estimate for coherence Is there any solution for estimating Confidence interval of coherence derived by Welch method with overlapping segments?

Moreover, the window length is not same for all the observations, only the number of segments are fixed for all observations

You might look at

G. Carter, C. Knapp and A. Nuttall, "Statistics of the estimate of the magnitute-coherence function," in IEEE Transactions on Audio and Electroacoustics, vol. 21, no. 4, pp. 388-389, Aug 1973.

doi: 10.1109/TAU.1973.1162487

Abstract: Expressions for the statistics of the estimate of magnitude coherence between two random stationary Gaussian processes are presented. These statistics include the probability density function, the cumulative distribution function, the bias, and the variance of the estimator. The expressions presented are in convenient and accurate forms for digital computer evaluation. Graphical examples of the bias and variance are included. Simple approximations are also given for the maximum bias, variance, and mean-square error.

keywords: {Density functional theory;Distribution functions;Gaussian processes;Mean square error methods;Probability density function;Random processes;Random variables;Statistical distributions;Statistics;Upper bound},

URL


Frequency domain samples to time domain without knowing the sampling rate? I have some frequency response data from 802.11a OFDM communication in channel 8 of the 5 GHz band, and I would like to inverse transform this to produce the corresponding time domain response. That is, I have $N$ equally spaced samples $H(f_1), H(f_2), ..., H(f_N)$ (except the middle sample is missing, see below) of the radio channel transfer function, and I wish to obtain the time domain impulse response $h(t)$, or in practice a discrete approximation $\hat{h}[n] \approx h[n] = h(nT_s)$.

Now, if I inverse FFT the samples I get a sequence $g[n]$, but I do not know what time values the elements of the sequence represent. The ways I know of to obtain the time resolution are all based around knowing the sampling frequency $f_s = 1 / T_s$ that was used to produce the original FFT, but that is unknown here. I do however know the frequencies $f_1, f_2, ..., f_N$, and intuition tells me that it should be possible to infer the time resolution from this knowledge. I am not very familiar with discrete transforms however, so I haven't been able to figure it out myself.

To be specific, the $N = 52$ frequencies in question are centered around $5040 \; MHz$ and spaced $312.5 \; kHz$ apart. The center frequency is missing (unused), so there are 26 below center and 26 above, making the full range $f_1 = 5031.875 \; MHz$ to $f_{52} = 5048.125 \; MHz$.

Generally the relationship between the sampling frequency $f_s$ and the frequency spacing of each bin $\Delta f$ is given by:

$$\Delta f= f_s/N$$

For example if you have 1000 bins and the sampling rate is $f_s = 1$ KHz, then each bin is spaced by 1 Hz given by $f_s/N$. So if the frequencies $f_1, f_2, ..., f_N$ (using the OP's indexing) were associated with bins 0 to $N-1$, each spaced by 1 Hz. In this case the frequencies would be $0, 1, 2... 999$ Hz and the sampling rate is $999$ Hz + 1 bin or $999$ + $1$ Hz = 1 KHz.

Update due to additions to the question:
In the specific case that the OP is using, the center unused bin would most likely represent the carrier frequency which would map to DC (and it is unused due to the possibility of interference due to carrier feed through). So here the center is likely to be bin 0, the frequencies to the right of center are likely to be bins 1 to 26, and the frequencies to the left of center are likely to be bins 27 to 52 (bin 27 is to the far left, representing the negative frequencies as would be mapped with fftshift in MATLAB/Octave). Regardless due to the cyclical nature of the FFT we can still use the guidance I provided to determine the sampling rate such that it would be N=53 total bins (including the nulled bin) and with a frequency spacing of 312.5 KHz the sampling rate would be 312.5 KHz x 53= 16.5625 MHz.


Recovering a signal after filtering A signal $x[n]$ is to be transmitted over a communication channel. The communication channel is described by an FIR filter of length $M$ such that the received signal
  is given by
  $$r[n]=\sum_{k=0}^{M-1} x[k]h[n-k]\text.$$

  The data are transmitted in blocks where each block contains N samples of a data
  stream $d[n]$, and in addition each transmitted block also contains a so-called cyclic
  prefix where the last values of the block are duplicated at the beginning of the
  block. Consider for example the first block of data ${d[0], d[1], . . . , d[N − 1]}$.

  The
  corresponding block to be transmitted with the cyclic prefix is
  ${x[n]}_{k=0}^{N+K-1}= {d[N − K], d[N − K + 1], . . . , d[N − 1], d[0], d[1], . . . , d[N − 1]}$
  $i.e. x[0] = d[N − K], x[1] = d[N − K + 1], . . ..$
  At the receiver end, the first $K$ values are discarded and the next $N$ values are saved in a buffer  ${y[n]}_{k=0}^{N+K-1}={r[K], r[K + 1], . . . , r[K + N − 1]}, i.e. y[0] = r[K], y[1] =r[K + 1], . . ..$
  a) Show that at the receiver end the original data stream d[n] can be recovered
  using the following algorithm if the cyclic prefix is sufficiently long.

  $Y [k] = FFT({y[n]})$
  $D'[k] = Y [k]/H[k]$
  $d'[n] = IFFT (D'[k])$

  Show all workings, and in particular determine the length $K$ of the cyclic prefix such
  that $d'[n] = d[n]$,$ n = 0, . . . , N − 1,$ and show exactly how $H[k]$ should be computed
  from $h[n]$.

To solve this question for the given info do we have to convert the sequence into odd and even sequence respectively and then arrive at the FFT or is there any alternative? I don't need the solution. A hint would be good to start.

don't need the solution. A hint would be good to start.

This is bog-standard cyclic-prefix OFDM, and that's really well-covered in literature.


Why is the high-pass filter result in a discrete wavelet transform (DWT) downsampled? From Wikipedia's description of the Discrete Wavelet Transform, a signal yields a set of:

approximation coefficients (low-pass: by averaging + downsampling)
detail coefficients (high-pass: convolving with wavelet + downsampling)

This process is cascaded over approximation coefficients:

Why is the high-frequency part of a signal downsampled and how does it retain information? Can it be reconstructed? Is information lost due to the Nyquist theorem?

Sources I looked at:

Mathworks: Wavelet decomposition algorithm
Guide for wavelet transforms for machine learning

Assume the original sampling rate is $F_s$ and that the filters are perfect brickwall Halfband filters. After the low pass filtering the frequency content is from $0 \rightarrow \frac{F_s}{4}$. Because of the new bandwidth, the sampling rate can be reduced to $F_s/2$.

For the high-pass filter, the frequency content lies in the range $\frac{F_s}{4} \rightarrow \frac{F_s}{2}$. There is no content in the range $0 \rightarrow \frac{F_s}{4}$. We can downsample by a factor of 2 in this case as well. The downsampling causes the frequency content in the $\frac{F_s}{4} \rightarrow \frac{F_s}{2}$ to be aliased into the range $0 \rightarrow \frac{F_s}{4}$, but because there is no frequency content in this range, due to the high-pass filter, we have not lost any information.


Transmit bitstream from USB over radio I am extremely new to this and I'm learning by building stuff myself so please let me know if this is too vague and Ill elaborate if I can. Im sure this is something very simple, but cant figure out how to do it.

I have built radios before and know how to modulate a signal as ASK but now I'm trying to send a simple bitstream from a USB port to a radios input.
My idea was to use a serial program to dump bytes into the radio input signal but all the research I have done indicates using rs232 or UART or something. I have no experience with this and I dont want some protocol sticking parity bits in all over the place.
I simply want to be able to send, for example, 0x1234 as
00000001000000100000001100000100 at a baudrate of 2048.

My intention is to monitor the frequency with an sdr and then learn how to demodulate the bits from the IQ data.

Is there a simple method for doing this or something I dont know about which will prevent me?

Any help or direction would be greatly appreciated.

Check out PySerial - a python library for serial communications.


FFT of a stretched vector Lets say I have a small size vector x=[a b c d]. Now I stretch this vector 3 times and I got x3=[a a a b b b c c c d d d]. What would be the relation between fft(x) and fft(x3)?

Conceptually you should split this into two steps

Up-sample by a factor of 3, i.e. x = [a 0 0 b 0 0 ...]. This results in a 3 times periodic repetition of the spectrum.
Convolve with a rectangular pulse of length three, i.e. h = [1 1 1]. This crates time stretched sequence you want. Convolution in time is multiplication in frequency, your three copies of the original spectrum get multiplied with a $sinc$ function.


Finding the equivalent filter H(u,v) in the frequency domain of a 3x3 spatial mask I'm trying to find the equivalent frequency domain filter, $H(u,v)$, of a 3x3 spatial mask that averages all neighbours of a point $(x,y)$ in said 3x3 neighbourhood excluding the point itself. So far, I have been able to calculate:
$$g(x,y) = 1/8 \cdot [f(x+1,y) + f(x+1,y+1) + f(x+1,y-1).... f(x-1,y-1)]$$  -(no $f(x,y)$)- but here is where I get stuck.
I've read that I can use linearity and translation properties to get $H(u,v)$ after performing DFT on $g(x,y)$ but for me that is currently a 'draw the rest of the owl' situation.
How do I get $H(u,v)$ from $g(x,y)$ and what form does it take?
I.e., what will the DFT of $f(x+1,y-1)$ look like and how do I then get $H(u,v)$ from that?

Hint: $g(x,y)$ averages all pixel values around  $f(x,y)$, at 8-neighbour connectivity distance 1, except $g(x,y)$. The classical uniform $3\times 3$ average is:
$$ m(x,y) =  \frac{1}{9}\sum_{i\in[-1,0,1]}\sum_{j\in[-1,0,1]}f(x-i,y-j)\,.$$
Thus: $$9m(x,y) = 8g(x,y)+f(x,y)\,.$$
Knowing the Fourier transform of the uniform $3\times 3$ mask, and that of the identity, and using linear superposition, you should be on track.


Frequency domain confusion I'm a bit confused with some details about frequency domain or studying systems with Laplace and Fourier analysis.

Let's say we have a system with a transfer function $H(s)$ with no poles and zero's at origin. We know that the transfer function is a ratio between output and input. Also, an input of $X(s)=1$ represents a single impulse (Dirac) in time domain, which can be viewed as a signal containing all frequencies in frequency domain. Also, applying this input to the system will lead to an output $Y(s)$ which is equal to $H(s)$ which is like the system's unique characteristic (its transfer function) and also can be viewed as a response to all frequencies. This signal $H(s)$ can be transformed back to time domain.

Let's say we have a system whose impulse response is something like this:

When I studied Fourier Transform we used to study system's response and Bode plot by analysing variations in $\omega$. So, for example, making $\omega=0$ we could find the system's DC gain. Making $\omega$ low values we can observe the system response to low frequencies and the same for high frequencies with high values in $\omega$. It was very useful to study filters, for example.

But note that all of this is only using the transfer function. And the transfer function may have an appearance like this above. And this does not make sense to me. If we think of the time->frequency transform, we should think that in frequency domain we are finding all frequencies which compose our signal in time domain, correct? This explains the Dirac containing all frequencies, for example. So, when making s=0 we should find H(s)=0 since this signal above does not have a DC level.

When we think about the transfer function as a system response to $X(s)=1$ it seems ok. But if we think about the transfer function as the same representation of this signal above, it does not make any sense. Also, if we apply a step signal, our system may respond like this:

We saw that for this system (with no poles at origin), we can find something like $H(0)=K$ which represents the system's DC gain. Applying a step response will lead us to a signal $\frac{H(s)}{s}$. So now making now s=0 will lead us to $+\infty$. Again, if we think of this signal being the one above but represented in frequency domain, making s=0 should left us a DC finit level since the time domain representation show us that the signal contain a limited DC level.

I know the final and initial value theorems but still I want to understand these relations between these domains.

Note that if we want to study a filter response, we should apply some signal to the input and see what is the output. Mathematically we should convolute these signals in time domain or multiply in frequency domain and see what is the resulting signal. In time domain it all makes sense since we just convolute and see the system's response. In frequency domain it also makes sense just like the case where you use $X(s)=1$ which represents all frequencies and multiply by the system and see the resulting signal. But when thinking about the same signal being represented in both domains, I think I'm messing things up.

I just want a help to understand relations between 's' variable and time domain. Making s=0 shouldn't lead me to the DC component of the signal in time domain? Generically, each value of 's' shouldn't lead me to 'how much' of this frequency is present in this signal?

in short:

A signal like the first one does not have any DC component so it goes to zero when $t=\infty$. So why its representation in frequency domain have a DC component when we put s=0?
Also, the second signal does have a DC level when $t=\infty$. But in frequency domain the signal goes to $\infty$ when s=0 (which I thought it should represent DC level).
It happens with Heaviside function too. Its representation in frequency domain is $\frac{1}{s}$ altough when s=0 it goes to $\infty$. But it has a DC limited component. How do we explain all of that?

The impulse response and the transfer function of an LTI system are different, sometimes extremely different (except for a specific shape of Gaussian).

The DC component of any impulse response is it's average value.  Your example time-domain squiggle seems to have a non-zero average value (integrate to determine the area under the curve), thus is has a non-zero DC component in the frequency domain.


Align low-res images of a moving feet I have a pressure platform that provides 64x128 pixel grayscale images (1px aprox 5mm) of a subject walking on top of it.
If I capture the whole step, I can obtain the average distribution of pressure for all the frames. Just by averaging (in fact the operator used is the max of each pixel) I obtain what I need. That is because at each timestep, each point in the foot always touches the same sensor.
The difficult case is when I put this pressure platform below a treadmill. From the point of view of the pressure platform, the user is "sliding", i.e. each sensor captures information from different parts of the foot at different times.
What this means is that there is an offset $o_i$ between $image_i$ and $image_{i-1}$.
Right now I have two algorithms: one local (that works in real time, providing a rough estimation) another global (very slow, though unoptimized).
The local one is just to check the offset (with subpixel accuracy of, say, 0.1px, so linearly interpolating the images) that minimizes the SSD between both images, that is:
$$
o_i = argmin_{o} ||translate(image_i, o) - image_{i-1}||
$$
which results in images like this:

The global method is to maximize a focus measure over the aggregated image, that is, if we take $O$ as the vector of all the alignments for a whole step consisting of $N$ images then:
$$
O^* = argmax_O = focus(\sum_i^N translate(image_i, o_i))
$$
Here I use the variance of the laplacian as focus measure, and it results in images like this:

As you can see, the alignment is much better, particularly noticeable in the front of the foot (the left part of the image).
To implement it, I initialize $O$ with the alignment from the local algorithm, then iterate through each element assigning the value that maximizes the global focus of the resulting image (leaving the rest fixed). I perform multiple passes until a maximum is reached (no change of global focus respect last iteration).
The following link has a video of a whole sequence, aligned using the local algorithm: https://imgur.com/a/CMolyUt
The displacement is only in the X direction, so only a positive real number is estimated. In the case of the video, that step is 182 images, so the parameters are 182 positive real numbers.
For the global method, I just do a brute force search over all the parameters, and reduce the search space by optimizing over a finite list of parameters (e.g. $0.3, 0.35, 0.4, 0.45, 0.5, 0.55$) which I know are in the range of possible offsets (just take min/max of the local algorithm output and add some tolerance).
The question is two-fold:

Is there any other algorithm/measure that I could try in order to solve the alignment problem
Is there another way of optimizing my current global algorithm instead of the brute force approach I am following?

Very interesting problem
I will focus on the second aspect of your question, how to improve it.
What I would do next, based on the experience you shared is to combine the two methods.
One problem of your first approach is that it accumulates the errors in the alignment between multiple images. A possible approach would be to use your focus measure to align the images iteratively.
$$ F_{i} = F_{i-1} + translate(image_{i}, o_{i}) $$
$$ o_i = \arg\min_{o}\quad  focus\left[ F_{i} + translate\left(image_{i}, o\right) \right]$$
This approach is maybe not the best at the beginning, but after the first step it have the advantage of not accumulating errors in the same way as if you compute $o_i$ as using only $image_i$ and $image_{i-1}$. Notice that $F_{i-1}$ aggregates all the previously aligned images, so if the one image is not aligned very well it will make so much difference since it sees how it was aligned before.
This is maybe not going to be as good as your global method. What I would try, to improve your global method is to make the displacement continuous with some type of interpolation, because a continuous optimization problem is much easier to solve than a discrete optimization problem, because can use the gradient information[1]. I am sure you are not doing the search for all the offsets independently that would be $s^N$, where $s$ is the number of possibilities in your search list.
Also instead of limiting the the displacements to a given range I think you could simply add some cost to the differences of the displacements $k \sum ||o_i - o_{i-1}||^2$, this makes sure that if have a high relative displacement between two frames the cost function will indicate that this is unreasonable.


Can I implement a LP filter using DMA? This questions seems to be a little strange but bear with me.

I have an application that needs a digital filter to be implemented in a microcontroller.

What I want to do

I want to use the DMA to get the data from the ADC every 300us (approximately). and store it in a buffer and I'm using FreeRTOS. I want to have a task check the DMA for half transfer completion and whenever it does reach that point, then the microcontroller would execute the LP filter on half the array, while the DMA keeps on filling the rest of the buffer. Due to my microcontroller's specs I'm sure I can run this filter in those points way before it has filled the other end of the buffer.

After this it would go to a waiting mode until it finishes the later half  and when it's done it will execute the filter on the second half.

Does this make sense?

What I have done in the past

Program an interruption , when the interruption occurs put the data in a FIFO and keep filling the buffer until you reach its limit while doijg the calculations

Question

Is "what I want to do" ? Of course it can be done, but will it actually be a proper filter? It sounds strange in my head when I think about it, but it somehow makes sense (ever get that feeling?). Or do I have to do it the way I've done it in the past?

Another option would be to do it in a mixture kind of way? Fill the buffer with DMA and execute FreeRTOS every milissecond and do the filter the same way I used to do it (fill FIFO, do calculations, etc)?

EDIT:

Thank you @Ben for your response. I'm adding the same info I put as a comment here in the edit so it can be better visualized.

Desired filter - IIR filter initially

Sampling rate- 10kHz. An ADC acquisition every 100us (waaay above the Nyquist frequency, but it's necessary for what I want .. a power meter)

Tap - between 2 and 10, it doesn't need to be a high order filter

Acceptable latency- Well I'd need to use the output to do al the calculations related to power and then store it and send it to storage somewhere, so the latency can be high

Basically you want to implement a filter with batch processing, right?

Is this an FIR filter or an IIR filter?

What is the sampling rate? Number of taps? Acceptable latency between output and input as batch processing adds some processing latency.


Determine the Signal Curve from Parameters of a Power Curve by Noisy Measurement I have a class of signals described by function:
$$
f(inc,d,t)=inc\cdot t^d
$$
where inc and d have a finite set of values like 1, 2, 3, i.e.
$$
inc, d\in \left \{1,2,3  \right \}
$$
and
$$
0\leq t<1
$$
Example plots:

I need to determine parameters inc and d given a discrete signal of f obscured by additive white Gaussian noise. Preferably an algorithm capable of working at real time.
Some candidates are:

sole differentiation (this amplifies noise)
filtering + differentiation (loses information)
curve fitting (least squares?)

Polynomial fitting seems to be an overkill, because function is known and is monotonic. Isotonic regression also seems too general, because it is applicable to all increasing/decreasing functions. Maybe there is a better solution when we have a concrete function and a set of possible parameters.

Since your domain for the parameters is limited (Only 9 options) the best way for White Noise would be going through them and pick then one with the least Mean Squared Error (MSE) which is the parameter to minimize for AWGN.
in MATLAB it will be something like:
function [ paramAlpha, paramBeta ] = EstimateModelParameters( vT, vY )

vParamAlpha = [1, 2, 3];
vParamBeta  = [1, 2, 3];

bestMse = 1e50;

vX = zeros(size(vY, 1), 1, class(vY));

for ii = 1:length(vParamAlpha)
    currParamAlpha = vParamAlpha(ii);
    for jj = 1:length(vParamBeta)
        currParamBeta = vParamBeta(jj);
        vX(:) = currParamAlpha * (vT .^ currParamBeta);
        currMse = mean((vX - vY) .^ 2);
        if(currMse < bestMse)
            bestMse     = currMse;
            paramAlpha  = currParamAlpha;
            paramBeta   = currParamBeta;
        end
    end
end

end

In a simple performance simulation I got the following:

Which means that for STD in reasonable value for this case (Since your maximum value is around ~3) your performance will be great.
The code is easily adaptable to any values of the parameters you'd like.
Though I think such brute force methods are reasonable for up to ~100 combinations.
Another thing you might try is use Affine Model on the log of the values.
Then apply some "Rounding" to the values. Though this will change the properties of the noise, in real world it seems to be effective way.
The full code is available on my StackExchange Signal Processing Q70753 GitHub Repository.


How can I compute a time-frequency cross-spectrum in MATLAB? MATLAB has a spectrogram function for the time-frequency analysis of a single signal. It also has a cpsd function for estimating the cross-frequency spectrum for two signals. However, cpsd averages across windows, collapsing the time axis into a single estimate.

Is there a function in this same family that does not average, but returns a time-frequency cross-spectrum instead?

Because I want to use this in a step-by-step tutorial for teaching purposes, I'd like to avoid the following two possible solutions:

The cohgramc function from the Chronux toolbox (uses a multitaper spectral estimation method I don't want to get into)
The wcoher function from the Wavelet toolbox (uses wavelets, idem)

Both of these methods return nice time-frequency cross-spectra, but I would like just to have a basic version with the familiar windowing parameters that spectrogram and cpsd use. Does this exist?

I don't know of a function that does all of what you ask, but it is easy enough to write.  They key step is to create short-time integrations on which the FFTs may be performed.  As you say, cpsd then averages all these STIs, but you can write your own to skip that step.

The key function to do this is $y = buffer (x,fftlen,overlap)$.  You input your FFT processing parameters, similar to cpsd, and get a matrix of overlapping time series.  I apply a Hanning window to each STI before performing FFT. Then you can simply apply $S = FFT(y)$ to get a matrix of overlapped spectra.

So it would look something like this:

% Assuming x is longer than fftlen:
overlap = fftlen/2; % 50% overlap
win = hanning(fftlen);
X = buffer(x,fftlen,fftlen/2,'nodelay');  % Matrix of overlapping STIs
numSTIs = size(X,2);
winX = X.*win(:,ones(1,numSTIs); % Time-domain windowed STIs
S = fft(winX,fftlen,1)/fftlen;  % Double-Sided Complex Spectrum Matrix
SdB = 20*log10(2*abs(S(1:fftlen/2+1,:)));  % Log Scale Single-Sided Real Spectrum Matrix

You can then calculate your auto-spectrum using either of these matices.  Of course you would need an x2 to repeat and calculate a cross-spectrum matrix.

for i = 1:numSTIs
    Sx1x2(:,:,i) = S1(:,i)'*S2(:,i);
end

(If anyone knows of a way to do this without a for loop, please let me know!)


A kind of Phase Retrieval problem I know there are lots of papers proposing algorithms for the problem of reconstructing a signal from modulus of its Fourier Transform (so-called Phase Retrieval Problem).

Also, recently it is studied for the case of sparse signals.

I want to know whether there is any paper or survey working on "reconstruction of signal from modulus of its DFT" or equivalently, reconstruction of signal from its circular auto-correlation ?

Or is there a way we can relate these two problems?

Edit: The prior information I have is :

Time-domain signal is positive-valued.
Time domain signal is zero in some points, i.e. signal has a support $P \subset \{0,1,2,\cdots,N-1\}$

Problems such as those are closely related to "inverse problems", whereby one attempts to reconstruct a signal from an under-determined system of equations, (ill-posed problems) given some apriori information. As @Hilmar mentioned, you cannot reconstruct phase information from the absolute magnitude of a DFT, unless you have prior information that constrains you. For example, in the paper you mention, the author states:

  Both the problem of phase
  retrieval from two intensity measurements (in electron microscopy or
  wave front sensing) and the problem of phase retrieval from a single
  intensity measurement plus a non-negativity constraint (in astronomy)
  are considered, with emphasis on the latter. It is shown [...]

The bold is mine. Even in the paper you allude to, there are constraints that then allow for a Maximum Likelihood framework to be applied, or one where sparsity is rewarded, etc.

  Edit: The prior information I have is : Time-domain signal is
  positive-valued. Time domain signal is zero for in some points, i.e.
  signal has a support P⊂{0,1,2,⋯,N−1}

All your signal being positive valued is irrelevant to the phase information of the entire DFT spectrum. To see this, consider this following:

Note that having an all positive signal $x[n]$ means you are guaranteed to have a DC bias. This means that the first bin in the absolute magnitude of you DFT corresponding to $|X(0)|$ will be strictly positive. Let us split the input signal into its zero-mean partition, ($z[n]$) plus the mean value, $m$. Then:

$$x[n] = z[n] + m$$

The DFT of the result is, (exploiting the linearity property of the DFT):

$$
X[k] = \sum_{n=0}^{N-1} x[n] \ e^{-j \frac{ \Large 2 \pi n k }{N}} = \sum_{n=0}^{N-1} \Big[z[n] + m \Big] \ e^{-j \frac{ \Large 2 \pi n k }{N}} = \sum_{n=0}^{N-1} z[n] \ e^{-j \frac{ \Large 2 \pi n k }{N}} + \sum_{n=0}^{N-1} m \ e^{-j \frac{ \Large 2 \pi n k }{N}}
$$

Recall that the DFT is nothing but a projection of your signal $x[n]$ onto an orthogonal basis set composed of complex exponentials. Save for the DC component at $k=0$, all the complex sinusoids have mean $0$. Thus the second term $\sum_{n=0}^{N-1} m \ e^{-j \frac{ \Large 2 \pi n k }{N}}$ completely vanishes for all $k \ne 0$. Hence, we can now write the DFT of $x[n]$ as being:

$$
X[k] = \left\{
  \begin{array}{l l}
     \sum_{n=0}^{N-1} z[n] \ e^{-j \frac{ \Large 2 \pi n k }{N}} + mN = mN & ,k = 0\\
    \sum_{n=0}^{N-1} z[n] \ e^{-j \frac{ \Large 2 \pi n k }{N}} & , k \ne 0
  \end{array} \right.
$$

This then shows that your signal always being positive is completely irrelevant to the phase signal spectrum.

Hence, this prior information ($x[n] \ge 0, n \in 0, 1, ... N-1$) is irrelevant, and no recovery is possible given $|X[k]|$, with only this constraint.

You are right in that these types of problems have some solutions to them, but it is also important for you to also remember that you have to define your constraints before you can solve for your unknown entities. That is, you have to give prior information that you can utilize.

As a simple example, given the absolute magnitude vector $|X(f)|$, there is an infinite number of possible solutions to what values $\angle X(f)$ can take on. However, if you "know" before hand that, say, all the imaginary co-efficients of $X(f) \in \mathfrak{C} $ are zero, then we can easily recover what the real co-efficients are up to a sign ambiguity.

The take home message however, is what prior information do you have, and how is it related to the DFT. Hope this helped.


FFT of size not a power of 2 My question is regarding the input size of a signal which is not a power of 2 and we have to take the fft of it. Some solutions say that suppose  if we want to take the fft of 1800 we should zero pad it till the length of 2048 to make it power of 2 and then apply the radix 2 algorithm. But there are other solutions as well which applies a combination of different algorithms without zero padding and then calculating the required FFT. My question is Does zero padding a signal to length of 2048 in case we have to take fft of 1800 makes any difference in results , if we use a combination of different algorithms to calculate the fft of size 1800. Would there be any difference or the result would be same.

The resulting FFT will be different: instead of calculating the FFT at the frequencies $2\pi n/1800\ \ $ for $n=0\ldots 1799\ \ $, you will be calculating them at $2\pi n/2048\ \ $ for $n=0\ldots 2047\ \ $. However, there is no degradation of the information.

Both approaches are correct: using 1800 or 2048. I'd use the "minimum energy solution" (i.e. the simplest, the laziest solution).  This would usually be using the 2048 length transform.

People tend to use radix-2 transforms because they don't know any better. There seems to be much mis-information about FFTs HAVING to be power-of-2. There is no such constraint.  Also, they probably don't know about about decent non-radix-2 algorithms, such as those available in FFTW and other libraries.

To see that the FFT of any length is information-preserving:

  Assume your data of length 1800 is in x.  Form X = fft(x,2048); (or some other  length different from 1800).
  Find xx = ifft(X);.
  See what sum(abs(x-xx(1:1800))); is.

See also this question and its answers.


Window functions and the Fourier transform I want to clarify the 'correct' way to use windowing and Fourier transforms.

My question is somewhat related to this one, but I have a few additional queries.

I have some real, non periodic signal in the time domain $y(t)$. Since this is a real signal, the data is naturally discrete.

I want to calculate the Fourier transform of this signal, $Y(f)$, within an observation time, $T$. I understand that

$$ \tilde{Y}(f) = \int^{+T/2}_{-T/2} y(t) \exp(-2 \pi i f t) dt $$

$$ = Y(f) \ast T \operatorname{sinc} (\pi f T)$$

where $Y(f)$ is the unwindowed FT of the infinite signal

To avoid the effects of the $\operatorname{sinc}$ function, windowing is recommended to taper the signal within the observation period $T$.

My question is: If I want to calculate $\tilde{Y}(f)$ computationally (e.g. using python/numpy),

Is it required for the time signal to be zero at the edges? Why so?
Does spectral leakage still occur if so?
Aren't we distorting the true time signal $y(t)$?
If we increase the observation time $T$ to e.g. $2T$, how does $\tilde{Y}(f)$ change?

So called spectral leakage always occurs when transforming any finite length observation.  That's because it's an artifact of windowing, and any non-infinite length observation (-t/2 ... +t/2) requires some sort of window, inherent rectangular or otherwise.  A window function allow you to choose your artifact: from window "distortion" to Sinc effects, and/or many other choices.


Low-Pass FIR that keeps all values positive Let's say, that I want to low-pass the norm (hence positive values)  of an analytic signal (related to Hilbert transform). This allows me to perform experiments with the envelope of audio.

For low-passing, I apply a Gaussian blur. It's technically the same as applying multiple time a moving average.
The advantages of a Gaussian blur are:

easy to implement, and enough fast.
good attenuation in the "stop band" (I just have to iterate more moving average to improve this).
all values remain positives. This is important, because  I use the low-passed norm to perform divisions.

But I'd like to test a low-pass with a narrower transition band.

Could you suggest me a simple solution, for a low-pass FIR (or something "linear phase" anyways) , with a narrower transition band than a straightforward Gaussian blur, and that allows to keep all values positive. For instance, with a sinc based FIR, there are oscillations around the zero line, and I want to avoid that. Adding an offset is out of question, as I want a result that is representative of the amplitude of the audio.

You can't, at least not strictly. Low pass filters with an all positive impulse response will always have a very wide transition band.

However, if you want to do audio envelope detection both Hilbert Transform and low pass filters are the wrong tools for the job: Hilbert Transformers don't work well on broad band signals and low pass filters mess up your transients and onsets.


Compute absolute power spectrum in MATLAB I have EEG signals. Is there a way to compute the absolute power spectrum in MATLAB within a frequency band? I know I have to use the fft function but I can't find any details about the absolute power spectrum.

If you have the Signal Processing Toolbox, you can use pwelch. If you do >>help pwelch   There will be some other choices. Try a few. Only you have your data so only you can determine if your results are acceptable

If you don't have the Toolbox, using FFT isn't hard. It's a matter of the length of your data, window choice, overlap, and a few other parameters


How do you find the frequency and amplitude from a DFT? I generated a time history of 200 points using the equation
$$x(t) = A\cos(2\pi f t)+ B\sin(2\pi f t) $$

with $A = 1$, $B = 0.1$ and $f = 17.2 \mathrm{Hz}$ and a sample rate of 400 samples per second. The numerical Fourier transform of the time history has a modulus plot of:

The question is how to get the frequency, amplitude and phase from this plot? This is the typical case where the frequency falls between bins. As there are three unknowns I would imagine you could get the precise values from, say, the three largest points which I call premax, max and postmax. However, complex values are available for all the points so many methods are possible.

This question arose as an output from a previous question where it was suggested I ask afresh and that several answers were available.  I have also asked this question on Mathematica Stackexchange here where it was suggested I ask on Signal Processing.

Many answers are available. Exact rather than approximate answers are preferable. Some, like using a Prony series or curve fitting a sine wave, do not use the Fourier transform. These are also interesting but I would like to use the Fourier data.

If there is noise on the data then some sort of fitting is probably needed depending on the extent of the noise. All suggestions welcome. Thanks

Edit

A bit more in response to the answers which I still have to fully investigate

There is an exact analytic formula for the Fourier transform of a sampled sine wave. For my data it plots as

This plot could also be obtained by calculating in-between values for the original data. I also calculated the frequency at the maximum of the analytic curve. To 25 significant figures it is

17.22726210702718945079653

which is significantly different to value of 17.2 used in the calculation of the  time history.

Some of the answers suggest using parabolic interpolation. I show the results of this here

I have tried parabolic interpolation using the complex values of the Fourier transform and the absolute values. The absolute values show a better maximum than the complex values but they are still well off. I still have other suggestions to try.

Methods based on time domain analysis look more promising than using Fourier which is disappointing.

For a lower quality frequency estimate, but still more precise than just using the (center) frequency of the peak magnitude DFT result bin, try using parabolic interpolation (using 3 neighboring points) of the spectral peak.

You can't get the phase from a DFT magnitude-only plot.  For phase you will need an estimator that uses the complex DFT results.

For a higher quality estimate, including a phase estimate, try using windowed-Sinc interpolation of the complex DFT result components, perhaps followed by successive approximation to find a better estimate of the magnitude peak near the peak DFT bin.  The quality of the estimate will vary with the width of the window.

Equivalent to Sinc interpolation would be to use a much longer zero-padded DFT or FFT to get more closely spaced frequency points.

There are a few added comments regarding the above frequency estimation methods on my blog:
http://www.nicholson.com/rhn/dsp.html#1

Note that near DC or Fs/2 (closer than the width of your interpolation kernel), you may need to use a mirrored interpolation kernel.  Note that hybrid methods are possible: e.g. using Sinc interpolation to up the resolution by 4X or 10X, then using parabolic or higher-degree polynomial interpolation to estimate the peak location, then using Sinc interpolation to estimate the magnitude and phase at that estimated peak frequency.  etc. Note that sinusoidal phase estimation in noise will often be more accurate when the phase is referenced to the center of the DFT window (N/2), not to sample point 0.  Referencing phase to the center of a data vector can be done using some form of fftShift.

The question does not ask about frequency estimators that do not involve a DFT, but note that in zero noise, only 3 or 4 unaliased points are required to solve for 3 unknowns (frequency, magnitude and phase of a pure unmodulated sinusoid).  See: http://www.claysturner.com/dsp/3pointfrequency.pdf and http://www.claysturner.com/dsp/4pointfrequency.pdf


Signal reconstruction in FMCW radar I am doing fmcw radar simulation. I have an output of mixer and from that I want to reconstruct the amplitude and phase of the signal. Therefore, I am taking hilbert transform and by the abs() and angle(), I can get the phase and amplitude information. Since the hilbert transform is in time domain, I want to map this time domain axis into frequency domain so I can see the value of amplitude and phase at different frequencies. FFT is not an option since by taking the fft, I will get the beat frequency.Therefore, how to map this time domain axis to frequency domain.

here are some of the relations between time and frequency. t = f/K, where

K = Bw / Tm
f = frequency,
Bw = bandwidth,
Tm = transmit time

I guess here, deltaf = 1/ t(2)-t(1) and fmax = 1/deltat, I am trying to map the time axis into frequency, but it gives dimension mismatch. Still there is something missing( some multiplyting factor or diving).Acoordig to my SV, its a demultiplexing process. That woud be really good if anyone can help me.

mixed1=(y.*(Y2));
hil = hilbert(y_filt);
figure, plot(t,abs(hil));
figure, plot(t,angle(hil));
deltaf = 1/t(2) - t(1);
fmax = 1/delta t;
f =?            %%% how to define frequency axis so i can plot amplitude & phase
figure, plot(f,abs(hil));
figure, plot(f,angle(hil));

I am pretty sure you will not be able to do what you are after without a FFT. I belive you exclude the use of FFT due to a common misunderstanding of the 'time' vectors involved in a FMCW radar.

First, your modulated FMCW chirp signal is in t = 0:T_m, if your channel is a single delay $\tau$, you now have

$y(t) = \sin(2\pi(f_l t + K/2\cdot{}t^2)$

$y_2(t-\tau)=\sin(2\pi(f_l (t-\tau) + K/2\cdot{}(t-\tau)^2)$

mixed = y.*y_2

Note that you can use a spectrogram on both the $y$, $y_2$ and mixed signals to get an idea of the frequency and amplitude.

Your only sensible option at this point is to take a FFT of the mixed signal to obtain the beat frequency (from matlab fft doc):

n = 2^nextpow2(length(t)*5); % zero pad to 5 times the length
Y = fft(mixed,n);
f = Fs*(0:(n/2))/n;

you may want to do a hilbert transform and also apply a window before the fft, but the principle stays the same.

It is trivial to derive that $y(t)\cdot y_2(t-\tau)$ will have a frequency term

$f_\text{beat} = K\tau$

hence we can plot the above in terms of two-way travel time $\tau$:

two_way_t = f/K;
plot(two_way_t, abs(Y)); % expect a peak at \tau
plot(two_way_t, angle(Y)); % expect a wrapped version of 2*pi*t*K*tau + 2*pi*(f_l*tau - K*tau^2)

This two_way_t time vector is now your 'compressed time' vector.

I really hope this clears something up, because I realise I have not really answered your question.


Why do we say that spectral efficiency is measured in [bit/s/Hz], not just [bit]? The link spectral efficiency of a digital communication system is measured in bit/s/Hz.

We have that  1 Hz = 1 s$^{-1}$.

So, $\left((1 \,\text{bit}) /1  \text{s} \right)\, / \left({\frac{1}{1 \text{s}}} \right) = 1 \, \text{bit}$.

Thus, why don't we just say then that spectral efficiency if measured in bits, but wherever I read it's bit/s/Hz? Is there some special meaning to this?

The units of spectral efficiency should be thought of as (bits per second) per Hertz.  It measures the ratio of the data rate to the bandwidth.  Yes, dimensional analysis shows that the dimension is equivalent to just bits but using that simplified unit for spectral efficiency does not show what the meaning of the term is. So, it is better to use (bits per second) per Hertz as the unit, or, for those whose (US) keyboards have broken shift keys and cannot type parentheses (the ( and ) come out as 9 and 0 respectively), as bits per second per Hertz or bits/second/Hertz.  (Typing / does not need the shift key).


What is matched filter in QPSK signal? Consider we have received a QPSK signal. We than make I (In phase) and Q (quadrature component) which I understand is nothing but multiplying the received signal by the cos(2*pi*f) and -sin(2*pi*f) Where f is the carrier frequency. Is this creation of I and Q is matched filter or there is some other filtering we do which is termed as matched filtering. I will appreciate any clarification?

Edit: Is it that we use some matched filter on I and Q data. What matched filter we use depends on the transmitter characteristic?

Well, you're not going to send and receive those symbols directly, you know?

What you'll need to convey the information over a channel is a waveform. Your receiver will then typically use the same waveform to correlate the signal. Say you are using root raised cosine pulses to transmit, then your receiver will employ a root raised cosine filter as a pulse matched filter to maximize the SNR.

There may be other things you can matched filter against. For instance, if you have an idea what your channel looks like, you might use a channel matched filter to maximize the SNR conveyed over your channel. Or maybe you are using spread spectrum techniques to transmit your information or things like code multiplex? You might want to use a code matched filter.

So, matched filtering can be used for a lot of things, the term alone doesn't tell it...


Fourier series of $cos(\omega_0 t)$ in continuous time Can any one please help me with understanding how we can calculate the Fourier series of Cos(w0t) using the formula:

I saw that they did the following calculus, but I Don't really understand how we calculate the integral of Xs[k] in this case:

It'll be so helpful if someone can solve it step by step.

They did not actually solve the integral, they just used a version of Euler's formula:

$$\cos(\omega_0t)=\frac12\left(e^{j\omega_0t}+e^{-j\omega_0t}\right)$$

which already is the Fourier series of $\cos(\omega_0t)$. Comparing $(1)$ to the formula for the Fourier series, you see that all coefficients are zero except for the ones with indices $k=1$ and $k=-1$, and they equal $\frac12$.