Range-Doppler matrix I am working on a radar simulation to calculate the range and the Doppler frequency of a target. I would like to build a range-Doppler map or range-Doppler matrix. I do not know the steps to develop the range-Doppler map or matrix. I hope there is someone out there kind enough to help me.

Its very simple: The time domain signal per range gate should be windowed (e.g. hamming, blackman-harris, etc.) to avoid ringing, then a FFT per range gate should be calculated. All that is left is to arrange everything in a matrix where each row is a frequency bin and each column is a range bin. (or vice versa).


Kalman Filter - Position/Velocity from Accelerometer and periodic position measurements I am trying to implement an EKF to estimate my position and velocity states by using accelerometer measurements as well as periodic GPS (position) measurements.  Basically I want to use the constantly use the accelerometer and re-calibrate it for every position measurement that is available.  When there is GPS dropouts, the estimates rely solely on the accelerometer readings.

I believe I have it working if all measurements are available.  Since my acceleration is not one of my states, my accelerometer measurements are used in the predict side of my filter.  The position measurements are in the update side like normal.

If a position measurement is unavailable, what do I do during my update?  Without a position measurement, am I basically just saying that a priori estimate is equal to my posteriori? (i.e. my "pretend" position measurement is just my estimated position).

Or is there something fancier that has to be done?  Is there issues with the system not being "observable" (i.e. observability matrix not full rank) when the position measurements are not available?

Thanks in advance!

I'm not sure how well it works, but one way to tweak the model is to set the measurement noise variance for the missing position measurement to a large value (meaning high uncertainty of the measurement), just for that missing time instant.

Setting it to be "infinite" should have the same effect as what you're doing now: just equate a priori and a posteriori estimates.

What setting the variance to some smaller value might give is a way to keep "momentum" in the position change.

I've not tried this approach, so please take with a grain of salt!


How to sample $f:=f_{1}\cdot f_{2}\cdot...\cdot f_{K}$? Let $f_{k}$ be $\sigma_{k}$-bandlimited, $k=1,2,...,K$ and define $f:=f_{1}\cdot f_{2}\cdot...\cdot f_{K}$; then how should I choose $T>0$ such that we can recover $f$ from samples $f(nT)$?

According to the Shannon sampling theorem, we must choose $T$ such that $\frac{1}{T}>2\cdot B$, where $B$ is the bandlimit. Hence the question boils down to how we determine $B$.

Now, we have to choose $B$ such that $\hat{f}(\xi)=0$ for $|\xi|>B$. If we compute

$$\hat{f}(\xi)=\mathscr{F}(f_{1}\cdot f_{2}\cdot...\cdot f_{K})(\xi)=\hat{f}_{1}(\xi)\ast\hat{f}_{2}(\xi)\ast...\ast\hat{f}_{K}(\xi),$$

then this still seems inconclusive to me, so I'd appreciate some help in this regard.

The definition of bandlimited is that $\hat f_k$ has compact support of size $\sigma_k$. To make our argument easier, we'll also note that we can decompose every $\hat f_k$ into $\tilde f_k * \delta_k$ (a convolution with a dirac impulse, just a shift) so that $\tilde f_k$ is "0-aligned", ie. $\sup f_k = \sigma_k$ and $\inf f_k=0$.

For $f_1\cdot f_2$ follows directly that the support of $\tilde f_1*\tilde f_2$ has the property

$$\text{supp}(f_1*f_2) \subset [-\sigma_2;\sigma_1]$$

from the way we calculate the convolution:

$$(f_1*f_2)(x)= \int\limits_{-\infty}^\infty f_1(y)f_2(x-y)\,dy$$

Now, that support has "length" (to be more exact, an upper boundary for measure) $\sigma_1+\sigma_2$.

It's very intuitive to now apply induction to show that this applies to all $k$.

So, in the end, $\text{supp}\left({\Large*}\left(\tilde f_k\right)_{k=1}^{K}\right)\subset [0, \sum\limits_{k=1}^K \sigma_k]$.

Since the shifting around in frequency domain doesn't theoretically (measure of $\delta_k0\,\forall k$) nor practically matter (compare: undersampling), we can just sample sufficiently by observing the sum bandwidth.

Note that this has very practical effects, for example, when sampling signals that are reflected by moving targets.


"Complex Matched Filter" What does a complex matched filter mean? . A complex bandpass filter can be formed by multiplying a low pass filter by exp(jwt) getting an analytic filter (one sided in the frequency domain)  which is complex, the filter will have complex coefficients but it is implemented as two filters with real coefficients (coefficients multiplied by cosines for the real and the same coefficients multiplied bu sines for the imaginary) but how about the matched filter which is supposed to be a low pass filter centered around DC? how does that one look like?

A matched filter is matched to the received pulse shape. If the matched filter is implemented in the baseband, i.e., after down-conversion, the received pulse shape can be complex-valued, e.g., if the channel frequency response is not symmetric with respect to the carrier frequency. In that case the matched filter is a complex-valued low-pass filter matched to the received (complex low-pass) pulse shape. The frequency response of such a filter is centered around DC but it is not symmetric with respect to DC.


FIR filter design with frequency sampling method - setting proper phase response for IDFT I'd like to implement freuqency sampling method for linear phase FIR filter design using IDFT transform.

My procedure goes like this :

determine desired magnitude response value in frequency points
I add linear phase response function with group delay (N-1)/2 to get complex frequency response values and equidistant points in frequency
use IDFT to calculate system's impulse response.

This approach works for odd-length FIR filter but I have problems with even-length FIR filters - I just don't  get proper result.

I guess I'm doing something wrong maybe with phase function or something else - but what ?

Thanks in advance,
regards,
Bul.

Your procedure is of course correct. It is in the details where things usually go wrong. One important thing is how to extend the desired frequency response beyond Nyquist taking the required symmetry into account. In order for the filter coefficients to be real-valued, the desired frequency response must satisfy

$$D[k]=D^*[N-k],\quad k=0,1,\ldots,N-1\tag{1}$$

where $*$ denotes complex conjugation, and $N$ is the FFT length (= filter length). From (1) it follows that for even $N$, the elements of the complex desired frequency response are

$$D_0,D_1,\ldots,D_{N/2},D^*_{N/2-1},\ldots,D^*_1$$

For odd $N$ you get

$$D_0,D_1,\ldots,D_{(N-1)/2},D^*_{(N-1)/2},\ldots,D^*_1$$
Note that in order for (1) to be satisfied, $D_0$ and, for even $N$, $D_{N/2}$ must be real-valued.

This little Matlab/Octave code works for even and odd filter lengths.

% frequency sampling design of linear phase FIR filter

N = 64;                         % FFT length = filter length
np = floor(N/2) + 1;            % number of independent frequency points
n = 0:np-1;
w = n*2*pi/N;                   % frequency vector
M = sin(n*pi/(np-1));           % some desired magnitude response
D = M.*exp(-1i*(N-1)/2*w);      % desired complex frequency response (linear phase)
D = [D,conj(D(N-np+1:-1:2))];   % append redundant points for IFFT
h = ifft(D);                    % compute impulse response
max(abs(imag(h)))               % should be very close to 0
h = real(h);                    % remove numerical inaccuracies

% check result
[H,w2] = freqz(h,1,4*N);
plot(w/2/pi,abs(D(1:np)),'.',w2/2/pi,abs(H))


How to find the difference equation directly from Direct Form II signal flow graph I am trying to solve for the difference equation of the following signal flow graph:

I am aware that Direct Form II can be converted to Direct Form I, which finding the difference equation directly is much easier.
What I have come up with is defining an intermediate signal $v[n]$ as the top middle node. I get the following result,

$$\begin{align}
v[n] &= x[n] + \tfrac{3}{5}v[n-1] - \tfrac{38}{75}v[n-2] - \tfrac{2}{15}v[n-3] \\
y[n] & =v[n]-\tfrac{3}{10}v[n-1]+\tfrac{1}{3}v[n-2] \\
\end{align}$$

As you can see, $v[n]$ is recursively dependent on itself. I know I could transform the two difference equations into their system functions and multiply together, find the total system function, and then transform back to their total difference equation. But for practice, I would like to know if their is another way to do this. I would like to somehow find the total difference equation without using their system functions.

Any advice would be appreciated! Thank you!

Darkling, these things are quite clearly explained in standard signals & systems textbooks. But I assume you have little time to read (as most undergrad courseware are full of homeworks, preliminaries, projects that most of which I believe become harmful as they consume (more than necessary amount of) a student's learning time and just completing a homework or project without full understanding makes what sense..? anyway)

Consider the LCCDE associated with a causal LTI system (such a system is said to be at initial rest which assures that all the initial conditions ard zero; i.e., no output occurs before an input is applied) :

$$ \sum_{k=0}^{N} a_k ~y[n-k] = \sum_{k=0}^{M} b_k ~ x[n-k] $$ The order is $\max\{M,N\}$.

Now the LCCDE can be considered to be a cascade of two LTI systems; the RHS represents an IIR-LTI system represented by $A(z)$ and LHS represents an FIR-LTI system represented by  $B(z)$. So the system can be represented with the following block diagrams:

$$ x[n] \rightarrow \boxed{ B(z) } \rightarrow v[n] \rightarrow \boxed{ \frac{1}{A(z)} }\rightarrow y[n] \tag{1}$$

or equivalently
$$ x[n] \rightarrow \boxed{ \frac{1}{A(z)} } \rightarrow v[n] \rightarrow \boxed{ B(z) } \rightarrow y[n]  \tag{2}$$

This follows from the fact that both cascaded systems are LTI and order of LTI systems can be interchanged which follows from the fact that convolution operator is commutative and associative; i.e., $$h_1[n] \star h_2[n] \star x[n] = h_2[n] \star h_1[n] \star x[n] $$

Now, the block diagram of (1) can be implemented with the following signal flow graph which yields a Direct-Form-I structure :

And, the block diagram of (2) can be implemented with the following signal flow graph which yields a Direct-Form-II structure:

The intermediate memory (delays) are merged and the signal flow graph is simplified into the canonical direct-form II form finally as (assuming $N > M$)

From this last figure, it can be deduced that your LCCDE relating the input $x[n]$to the output $y[n]$ is :

$$ y[n] - \frac{3}{5} y[n-1] + \frac{38}{75} y[n-2] + \frac{2}{15} y[n-3] = x[n] - \frac{3}{10} x[n-1] + \frac{1}{3} x[n-2] $$


Does Reed Solomon successful decoding guarantee valid payload? I've got a CCSDS message. It was decoded by viterbi, de-randomized and processed by Reed Solomon decoder.

Decoder said that 14 bits were corrected. However when I decode payload according to the protocol, I'm getting invalid data.

So it looks like I've got random noise, which was accidentally decoded by Reed Solomon? Can it be possible? Or am I missing something?

No, you're missing nothing.

Quite obviously, any valid code word (sequence) should decode correctly. And since noise can lead to any potential code word (sequence) at the receiver, these "false corrects" are as likely as any other reception, if the reception was uncorrelated to the transmission.

It's just that codes are designed so that it's unlikely to happen that by pure (bad) luck you end up decoding the wrong code word when a code word was sent. But yes, it's possible, and yes, it does happen.


The necessity of null subcarrier in OFDM? In looking to LTE specifications, with the subcarrier spacing $\Delta f = 15 \textrm{kHz}$,

for bandwidth $10 \textrm{MHz}$, fft size $N_{fft} = 1024$ that needs a sampling rate $F_s \geq 1024 \times 15\textrm{kHz} = 15.36 \textrm{MHz}$. The useful subcarriers are the 600 centered ones (excluding the DC). There are $1024-600=424$ null subcarriers.
for bandwidth $20 \textrm{MHz}$, fft size $N_{fft} = 2048$ that needs a sampling rate $F_s \geq 2048 \times 15\textrm{kHz} = 30.72 \textrm{MHz}$. The useful subcarriers are the 1200 centered ones (excluding the DC). There are $2048-1200=848$ null subcarriers.

I want to know what is the objective of band-edge null subcarriers in LTE, or in OFDM ? I have found the photo above in Internet, it seems that these null subcarriers are guard band.

Question 1 : are these null subcarriers really sent over propagation channel ? I mean the occupied bandwidth of LTE in "$10 \textrm{MHz}$" case is really $10 \textrm{MHz}$ or $15 \textrm{MHz}$ ?

Related question, BW of OFDM question, but the answer is still unclear about what is the occupied bandwidth in the shared wireless medium.

Question 2 : if the occupied bandwidth is only $10 \textrm{MHz}$, can other systems use this bandwidth (from subcarrier 300 to 511 and from -300 to -511) ? Is there degradation of its own performance ? Does it cause inteference to LTE data subcarriers ?

Question 3 : if the occupied bandwidth is $15 \textrm{MHz}$, the factor $424/1024 = 41.4 \%$ seems a lot. How is it dimensionned ? Are they "guard band" to allow the side lobe of sinc pulse shaping (from FFT) to vanish ?

For the question 3, I have found this question in stackexchange a related question with the answer

  Specifications will often call for some null carriers around the edge
  of the band (to allow for rolloff in the antialiasing filters in the
  front end)

Hm, $41.4\%$ seems too much for rolling off, doesn't it ?

  Update : my questions are about the band-edge null subcarriers, but one may find interesting information about the DC null carrier in
  the answer. Note that the DC subcarrier is skipped at OFDM baseband
  signal generation, i.e. not be used to transmit data. But in wireless
  environment, there may exist DC carrier to help mobile devices to find
  the carrier frequency, at least it is true for Downlink LTE (see LTE
  UE radio TX and RX spec, figure 5.4.2-1).

The only portion of the spectrum that is "occupied" or "transmitted" are the subcarriers that are actually used (i.e. 10 MHz instead of 15 MHz). The null subcarriers referenced in your question are used for two purposes:

The null subcarriers along the edges are used to allow for easily-realizable anti-aliasing filters. Given a sample rate of 15.36 MHz (complex), and a passband of 10 MHz, it's relatively easy to design a good, simple lowpass filter that does the job. Simple is the key here; LTE receivers are made by the millions, and any cost savings that can be achieved by design of the waveform are attractive to hardware manufacturers.
The null subcarrier at DC is used to allow the use of simple (again, cheaper) direct-conversion (zero-IF) RF receivers. One disadvantage of zero-IF receivers is that they typically encounter a strong interfering signal at DC. This can be caused by direct leakage of the LO signal into the output as well as a nonzero DC offset in the receiver's analog-to-digital converter. This is sidestepped in LTE by just not transmitting a subcarrier at the center of the band.


Fourier transform of Autocorrelation function, what am I missing here? (Power spectrum) Given this problem

I know that the autcorrelation and power spectrum are fourier pairs, so when taking taking fourier transform of Rxx, one should end up with Sxx

However, when I take the transform of $R_{xx}(m)$, I end up with $S_{xx}(\omega)=\frac{2\alpha }{\alpha ^{2} + \omega ^{2}}$, which also matches up with fourier transform pairs in the book. Here's by hand:

$S_{xx}(\omega )= \int_{-\infty}^{\infty} R_{xx}(m)e^{-i\omega m} dm=\widehat{R}_{xx}(\omega)$

$=\int_{-\infty}^{\infty} R_{xx}(m)e^{-i\omega m} dm=\int_{-\infty}^{\infty} e^{-\alpha |m|}e^{-i\omega m} dm=$
$=\int_{-\infty}^{0} e^{\alpha m}e^{-i\omega m} dm + \int_{0}^{\infty} e^{-\alpha m}e^{-i\omega m} dm=$

$ =\int_{-\infty}^{0} e^{m(\alpha -i\omega)} dm + \int_{0}^{\infty} e^{-m(\alpha +i\omega )} dm=$

$=\frac{1}{(\alpha -i\omega)} - \frac{0}{(\alpha -i\omega)}+\frac{0}{(\alpha +i\omega)}-\frac{1}{(\alpha +i\omega)}$

$=\frac{1}{(\alpha -i\omega)}+\frac{1}{(\alpha +i\omega)}=\frac{\alpha +i\omega+\alpha -i\omega}{(\alpha -i\omega)(\alpha +i\omega)}=\frac{2\alpha}{\alpha^{2} +\omega^{2}} $

Which doesn't look like equation given in the problem, unless I'm missing some identities or something here.

Also, as for problem h), that should be pretty straight forward. Avg. power of $X=R_{xx}(0)$ right?

That question was posted quite some time ago! Have you tried the discrete time Fourier transform of $R_{XX}[m]$, instead of the Fourier transform? I believe this is what was meant, there (the text talks about a "sequence").

An "easy" way to solve this, according to me, is to compute the Z-transform, then identify $S_{XX}(\omega) = Z[R_{XX}](z=e^{j\omega})$. The Z-transform is given as:
\begin{align}
Z[R_{XX}](z) &= \sum_{m=-\infty}^{\infty} R_{XX}[m] z^{-m} \\
   &= \sum_m e^{-\alpha |m|} z^{-m} \\
   &= \frac{1 - e^{-2\alpha}}{(1-e^{-\alpha}z^{-1})(1 - e^{-\alpha}z)}
\end{align}

Hope that still helps!


Maximum NCO frequency and filter design for down-conversion of complex input signal I am trying to perform digital down conversion by using the DDC from Matlab dsp toolbox. My question is that the bandwidth property is limited to Fs/2 and the filters created are real filters. If the input is comlex then the frequency range will be from 0 to Fs and the decimating filters need to be complex but the toolbox does not allow that. Is this some issue with the Matlab toolbox or is this a theoretical limit that DDC NCO frequency cannot exceed Fs/2 even if the input signal is complex?

Given the input signal is complex, then this can possibly be down-converted with a single real NCO as in the block diagram as follows as long as the frequency distance from DC or the sampling rate to the digital IF carrier is greater than the half the signal bandwidth (meaning none of the spectrum crosses those boundaries). By doing so increases the challenge of the low pass filter for isolating the images, especially when the signal is close to the DC or sampling rate boundary.

This is intuitively explained through observation of the "unrolled digital spectrum", such that we extend the frequency spectrum to plus and minus infinity and know through sampling that the spectrum from DC to the sampling rate will repeat in every additional integer multiple of the sampling rate. And further for real signals, we known that the positive and negative frequencies must be complex conjugate symmetric. We see all this in the diagrams below for an arbitrary complex IF waveform that is above half the sampling rate. These diagrams show the frequency domain spectrums for the down-conversion process with complex NCO's as typically done (and what happens if we rotate the wrong way), and finally what the above block diagram would result for a digital baseband spectrum.
The top diagram would be the most typical approach, using a complex NCO (not as in the block diagram above, but a full complex multiplier requiring four real multipliers and two adds, together with the complex NCO). The second diagram is what would occur if the NCO was wired up incorrectly (Sin/Cos swapped), and the bottom diagram is what we could do with the simpler diagram above, and shows the greater challenge in filtering out the image (but otherwise quite feasible). Note that because the real NCO signal actually creates two tones in the unique frequency span of plus or minus half the sampling rate (or equivalently  DC to the sampling rate), we can indeed use the real NCO to down-convert any signal in the span from DC to the sampling rate (as long as the signal's bandwidth doesn't cross those boundaries as already noted). We see this directly with Euler's formula: $2\cos(\omega t) = e^{j\omega t} + e^{-j\omega t}$; a single sinusoid consists of two "tones", defining a "tone" as an impulse on the frequency domain.

The Complex IF is multiplied in time by the NCO, so in the frequency domains these spectrums would convolve resulting in the shifted spectrums at complex baseband. I find some grasp this quicker by knowing first that each impulse shown in the frequency domain represents a single spinning phasor in the time domain (given by $e^{j \omega_o n}$), so if you imagined an impulse at the center carrier for the complex IF given by $e^{j \omega_o n}$, you can see that it would be the one impulse out of many that is located at $e^{-j \omega_o n}$ that would move our signal to baseband, since in the time domain the product would be  $e^{j \omega_o n}e^{-j \omega_o n} = e^{j0}$, (and similarly you will see there are infinite combinations that do the same since the spectrums as shown repeat, but ultimately we know that everything from DC to fs repeats, so once we solve for all the signals in DC to Fs, we can simply populate the rest in this graphical approach). I highlight one example in the top diagram with the carrier at $f_c$ and the specific component in the complex NCO spectrum that is created at $-f_c$ that moves the complex IF to baseband through the time domain multiplication process. Through further observation we also see the same for the component in the fourth Nyquist zone just below $2f_s$, as we can find the component in the complex NCO spectrum that is the negative of this that would do the same (again this is just a visualization of the one digital spectrum that is unique over any frequency span of the sampling rate, such as $-f_s/2$ to $+f_s/2$ of $0$ to $f_s$, another approach is to view it as a cylinder but I find this clearer as can be represented easily on 2-dimensional plots). The math is very straight forward but given alone may mask this intuitive insight which helps in explaining many multi-rate and mixed-signal (analog/digital) principles.
A complex NCO will create a unique tone given by samples of $e^{j 2\pi f t}$ for any frequency from DC to the sampling rate (the sampling rate is an alias of DC, so to be accurate it would be unique from DC to one frequency spacing (given by $f_{\Delta} = f_{clk}/2^{acc}$ where $f_{clk}$ is the NCO clock rate and $acc$ is the accumulator bit size.). This would require an NCO with two outputs representing the real and imaginary components. For implementation this can be done with two real NCO's synchronized in quadrature given the relationship:
$$e^{j 2\pi f t} = \cos(2\pi f t) + j\sin(2\pi f t)$$
Implemented with "I" and "Q" where the $j$ is implied in the "Q" path. Ultimately the sinusoidal look-up-table would be shared with a simple offset in the phase accumulator of a quarter cycle (and only a quarter cycle needs to be stored to implement a full cycle with efficient up/down counting and sign change).
A complex NCO frequency translating a complex input is shown in the block diagram in this post:
Recovering signal for psk
Often with a DDC implementation, the input is real and the NCO is complex with a product of real and complex with complex output. The positive and negative half spectrums both get translated in the same direction (either left or right) such that one will end up at baseband (if shifted to the right, meaning the negative half spectrum, a spectral inversion will result which is easy to correct). The remaining higher frequency component needs to be filtered out after the multiplication.


Why does this manual bilinear transformation yield different results from Matlab's? I've got a first-order Butterworth filter with the cutoff frequency $\omega_c$. Its transfer function is then

$$H(s) = \frac{\omega_c}{s+\omega_c}$$

Using the bilinear transform to find an $H(z)$ (what is that function called?), I get

$$H(z)=\frac{\omega_c}{\frac{2}{T}\frac{z-1}{z+1} + \omega_c}
= \frac{\omega_c z + \omega_c}{\left(\frac{2}{T}+\omega_c\right)z + \omega_c-\frac{2}{T}}$$

However, I can't reconcile this result with what Matlab is doing. It seems wrong, no matter what value of $T$. I assume that B and A below are the coefficients of $H(z)$.

>> [B,A] = butter(1,0.5)
B = 0.5000    0.5000
A = 1.0000   -0.0000
>> [B,A] = butter(1,0.6)
B = 0.5792    0.5792
A = 1.0000    0.1584
>> [B,A] = butter(1,0.7)
B = 0.6625    0.6625
A = 1.0000    0.3249
>> [B,A] = butter(1,0.8)
B = 0.7548    0.7548
A = 1.0000    0.5095

What am I misunderstanding?

A couple things:

Before making the substitution $s = \frac{2}{T} \frac{z-1}{z+1}$, you need to prewarp the cutoff frequency by making the substitution:

$$
\omega_{c,w} = \frac{2}{T} \tan (\omega_c\frac{T}{2})
$$

where $\omega_{c,w}$ is the warped cutoff frequency. This is necessary because the bilinear transform maps the left-half plane in the Laplace domain (used in analog filter design) to the unit circle in the $z$-domain in a nonlinear fashion. Therefore, as you approach the Nyquist rate (digital frequencies of $\pm \pi$), the approximation to the analog filter prototype becomes inaccurate.

Also, the second parameter that you are passing to the butter function is the normalized cutoff frequency, not the sampling interval $T$. The normalized frequency used by that function is in the interval $(0,1)$ and is equal to the ratio of the desired cutoff frequency to the Nyquist rate:

$$
\omega_n = \frac{\omega_c}{2 \pi \frac{f_s}{2}}
$$

$$
\omega_n = \frac{\omega_c}{\pi f_s}
$$

$$
\omega_n = \frac{\omega_c T}{\pi}
$$


Convolution of Signal with a Gaussian Filter / Kernel Following up on Analytical Solution for the Convolution of Signal with a Box Filter, I am now trying to convolve a Gaussian filter with the sine signal by hand.

My method is to use the definition of convolution and attempt to integrate:
\begin{align*}
  \overline{\phi} &= \phi(x) * h(x) = \int_{-\infty}^{+\infty} \! \phi(x') h(x - x') \, \mathrm{d} x' \\
  \Rightarrow \overline{\phi} &= \int_{-\infty}^{+\infty} \! \sin(x') \left( \frac{6}{\pi \Delta^2} \right)^{1/2} \exp(-\frac{6 (x - x')^2}{\Delta^2}) \, \mathrm{d} x' \\
  \overline{\phi} &= \left( \frac{6}{\pi \Delta^2} \right)^{1/2} \int_{-\infty}^{+\infty} \! \sin(x') \exp(-\frac{6 (x - x')^2}{\Delta^2}) \, \mathrm{d} x' \\
\end{align*}
At this stage, the integration seems impossible by hand (e.g., integration by parts).

I also tried moving to wave space by taking Fourier transforms and instead using multiplication:

\begin{align*}
  \mathscr{F}\{\phi(x)*h(x)\} &= \mathscr{F}\{\phi(x)\} \cdot \mathscr{F}\{h(x)\} \\
\end{align*}

I get the Fourier transform of sine as $\sqrt{\pi/2} i (\delta(k - 1) - \delta(k + 1)$, but am not sure how to multiply this with the exponential and how the inverse Fourier transform will work.
Am I missing something clever and/or elegant? Is such a solution possible? I don't have a good understanding of the Dirac delta and how to multiply it with an exponential in wave space.

If you have an implementation you want to validate the best way to do it is by validating it versus a reference implementation.

For instance, the above, in MATLAB, would be something like:

% Setting the Grid Parameters
leftBound = -10;
rightBound = 10;
numSamples = 10000;

% Signals Parameters
gaussianKernelStd   = 2;
gaussianKernelMean  = -1; %<! Basically its shift on the grid
sineAmp             = 0.1;
sineFreq            = 1;
sinePhase           = 0;

% Generating the grid
vX = linspace(leftBound, rightBound, numSamples + 1);
vX = vX(:);

samplingInterval    = mean(diff(vX));
samplingFrequency   = 1 / samplingInterval;

% Generating the Sine Signal
vS = sineAmp * sin((2 * pi * sineFreq * vX) + sinePhase);
% Generating the Gaussian Kernel
vK = exp(-((vX - gaussianKernelMean) .^ 2) / (2 * gaussianKernelStd * gaussianKernelStd));
vK = vK / sum(vK * samplingInterval); %<! Normalization

vY = conv(vS, vK, 'same');

figure();
% plot(vX, vK);
% figure();
% plot(vX, vS);
plot(vX, [vS, vK, vY], 'LineWidth', 2);
legend({['Sine Signal'], ['Gaussian Kernel'], ['Convolution Result']});

The result is given by:


How to apply a 1D filter (CSF) to an image? I am facing the following problem:
I have an image and a contrast sentivity function (CSF). The CSF is a function defined in the frequency domain and it is just a band-pass filter. The problem is that it is a 1D filter so I do not know how to apply it to the Fourier transform of the image.

I thought about creating a 2D filter based on the revolution around the zero frequency of the 1D filter but I do not know if it is mathematically correct or if it is the common procedure.

Thanks in advance!

I believe what you want to do is take your 1D filter column vector, call it $x$, and create a 2D filter with $xx^T$.  The result of this outer product (column vector times row vector) is a matrix that you can use to do element-wise multiplication with the 2D Fourier transform of the image.


interpret wavelet scalogram My knowledge of wavelets is less than epsilon. Bear with me.
If I have a signal of two well separated sinusoids (15 and 48 Hz) plus some random noise, I can clearly make out the two in a spectrogram (the two stripes in my picture);

t=0:0.002:1;   % fs = 500 Hz
x=4*sin(2*pi*15*t)+2*cos(2*pi*48*t);
xn = x + randn(size(x));
figure(1);
plot(xn);
figure(2);
spectrogram(xn, 64, 60, [], 500);

but, using 'wscalogram', I cant say I can tell the components apart:

    coefs = cwt(xn,1:64,'db8','scalCNT');
    wscalogram('image',coefs,'scales',1:64,'ydata',xn);

Is there a way to read from the scalogram that there are 2 distinct sinusoids and if so, how can I separate them using wavelet decomposition and filtering? I thought that perhaps my choice of wavelet ('db8') is not optimal but I can't see much difference using other types ad larger (or smaller) number of scales. Clearly I am missing something about where and when and on what to apply wavelets.

Thanks

Wavelets are ideal for localized events. The Fourier Transform represents a function as a sum of sines and cosines, neither of which are localized. The spectrogram does keep some time information, at the expense of frequency resolution

In your case, the signal is not localized at all.  The spectrogram smears your 15 Hz band over several Hz, as it captures some time information, and the scalogram flat out fails.


LTI and Deterministic Channel Can someone explain me the statistical relationship of random signals with a linear time-invariant and deterministic channel to the output of the channel when having random signal as an input?

If your input "random signal" is modeled as a zero-mean wide-sense-stationary random process $\{X(t)\}$ with autocorrelation function $R_X(t)$, then the output of the channel (which is basically an LTI system with impulse response $h(t)$) is also a zero-mean wide-sense-stationary process $\{Y(t)\}$ with autocorrelation function $R_Y(t)$ given by
$$R_Y = R_h \star R_X$$
where $R_h$ is the autocorrelation function of the impulse response $h(t)$ $\big($that is, $R_h(t) = \int_{-\infty}^\infty h(t)h(t+\tau) \, \mathrm d\tau$ or $\int_{-\infty}^\infty h(t)h(t-\tau) \, \mathrm d\tau$ for left-handed folks$\big)$ and $\star$ notes convolution. Indeed, since both $R_h$ and $R_X$ are even functions, that convolution can be written as a cross-correlation if you like, which makes for a nice mantra to murmur to impress your boss during your presentation: "The output autocorrelation function is found by cross-correlating the input autocorrelation with the channel autocorrelation".

If you prefer power spectral densities, then
$$S_Y(f) = |H(f)|^2S_X(f)$$
where $H(f)$ is the channel transfer function.


Eroding noisy borders of polygons in an image I am new to image processing, and this is the task at hand.
I have a mask image of shape 500x500x1 that contains polygons.
These polygons (=vector data: common term in remote sensing and GIS) are labeled, and will be used as my target later on to classify crops. This vector data has been rasterized (another remote sensing and GIS term) and resulted in the 500x500x1 mask. As you can see the borders of all polygons are noisy (contains pixels of different colors==> different classes) although each polygon should correspond to one class (one color). I am guessing that something went wrong during the rasterization process;  I tried to erode the borders using simple morphological operations:
kernel  = cv2.getStructuringElement(cv2.MORPH_RECT, (2, 2))
erosion = cv2.erode(lulc_mask, kernel, iterations = 3)

However, this resulted in reducing the area inside the polygons and not the border.
Each time I increase the size of the kernel I basically erode the data from within the polygon, but the border remains the same.
Any idea on how to solve this?

I have a mask image of shape 500x500x1 that contains polygons.

. . .

These polygons [...] are labeled

These two statements are incompatible. Assuming that the third number is the color depth, then 1-bit implies a binary mask, but each class being assigned to a different level implies a color depth >1.
With a 1-bit mask, binary erosion is applied. Otherwise grayscale erosion is applied, which is why you see those "out of class" values at the borders of the mask.
The simplest way to deal with this is to apply the erosion at each class within the image iteratively. So, collect the unique values from the mask in unique_classes. Then, for every level in unique classes, threshold the mask image (this returns a binary mask) apply the erosion and then bring it back to the level it is supposed to be.
Alternatively, duplicate the mask into two images, hold one back to be the one that describes the classes and threshold the second one so that it becomes a binary mask. In this setting, you apply erosion once, to the second mask, but use the first mask to decide which class a given pixel is a member of.
Hope this helps.


Applying a window function to a speech signal I'm currently developing a speaker recognition system. I have two main problems:

What is the most suitable windowing function to apply to a speech signal prior to the STFT?
Is it necessary to apply more than one windowing function?

I don't really understand what "apply a suitable windowing function" "prior to the STFT" means. Most STFT implementation incorporate apodizing windows. And the best is always related of some cost function, objective or processing.

Note that there are people called Mr Hanning, but they are not (AFAIK) at the origin of the raised-cosine or von Hann window (Julius von Hann, Austrian meteorologist). You can find a mention at page 97 from the book  The measurement of power spectra from the point of view of communications engineering (archive.org), 1958, by Blackman and Tuckey. Noteworthy, they spell hamming and  hanning in small letters.

Hamming or Hann windows are very standard. $50\%$ or $75\%$ overlap are common. In  the recent works I know of, especially for blind source separation,  people often use at least two window sizes: a shorter one for transients, and a longer one.

This is a first step toward using different windows depending on the scale of observation. Some have used a multiscale context, such as in the
ERBlet non-stationary Gabor filterbank, ERB for Equivalent Rectanguar Bandwidth, with a constant-Q transform:

Some have been using asymmetric windows, better suited to the skewed shape of speech onsets (see eg Using asymmetric windows in automatic speech recognition).

The choice of the best windows, to me, even for speech, is not  clearly settled.


How to make a Power Spectral Density Plot in R I have a time series point process representing neuron spikes. I have computed and plotted autocovariance using acf but now I need to plot the Power Spectral Density.

Power Spectral Density is defined as the Fourier Transform of the autocovariance, so I have calculated this from my data, but I do not understand how to turn it into a frequency vs amplitude plot.

I have used the following code

# X is some set of Wait times between spikes, below is just an example
X <- c(56, 3, 4, 119, 3, 4, 121, 3, 3, 121, 3, 4, 120, 3, 4, 4, 115)
acf <- acf(X,type="covariance")
psd <- fft(acf$acf)

Now psd is a complex valued array across the default 24 lags from the acf function.

How do I turn this array into a PSD plot?

From further research I've discovered that the frequency is given by the index of the FFT multiplied by the sampling rate and divided by the size of the array. And the amplitude is the magnitude of the complex number.

So the full code for such a plot would be as follows

# Load ggplot library
library(ggplot)
# X is some set of Wait times between spikes, below is just an example
X <- c(56, 3, 4, 119, 3, 4, 121, 3, 3, 121, 3, 4, 120, 3, 4, 4, 115)
acf <- acf(X,type="covariance")
ft <- fft(acf$acf)
freq <- (1:nrow(ft))*1000/nrow(ft) #In my case we sample by ms so 1000 hz
A <- (Re(ft)^2 + Im(ft)^2)^.5 #amplitude is magnitude
PSD <- cbind.data.frame(freq,A)
ggplot(PSD, aes(x=freq,y=A)) + geom_line()


I have run an FFT now want to get the mean value in the 8hz - 13hz range? I have a  matrix [66 x 100,000] (it's some EEG data, across 66 electrodes) and want to find the mean power in the 8-13hz range. How can I do this?

I have so far: x = fft(data) What do I have to do next?

Furthermore, when I plot this I get the following::

I was hoping for something more like this:

Did you ever tried to follow example from MATLAB help? Result of the FFT is complex, containing information both for magnitude and phase. Thus in order to extract the magnitude you must take the absolute value of your result - abs function. What you are doing currently is plotting complex values.

Please refer to basics of DFT if you have any doubts (i.e. to understand why DFT spectrum is mirrored) and mentioned example for MATLAB fft. Here is one of ways to create frequency vector for magnitude response.

% fs --> sampling frequency to be defined.
N = length(data);
X = abs(fft(data))/N;
f = (0:N-1).*(fs/N);

If you wish to show your results in decibel scale, then convert the values with usage of 20*log10(X/ref); where ref is your reference value - can be max(abs(X)).

I am sure that having the frequency vector and magnitude values you can take it from here in order to calculate the mean value in appropriate range.

Extra things you might want to consider:

Apply the windowing to your signal if it is necessary (probably it is). When windowing is applied, scale not by N, but use the sum of window samples.
Use the Power Spectral Density estimation instead of FFT. Second figure you provided is the example of PSD - in fact this is copied from MATLAB help regarding PSD estimates using FFT. I believe you have all necessary info in there regarding periodogram.


What is the equivalent smoothing function to running the same Gaussian 8 times? Suppose there is an image which is to be smoothed by convolving it with a Gaussian kernel with standard deviation $\sigma$. If the image is then smoothed with this kernel 8 consecutive times, is the total smoothing still Gaussian with some calculable standard deviation $\hat{\sigma}$? Or is it some other distribution instead?

From Wikipedia:

  Applying multiple, successive gaussian blurs to an image has the same effect as applying a single, larger gaussian blur, whose radius is the square root of the sum of the squares of the blur radii that were actually applied.

So in your case, that would be equivalent to blurring the image with a single Gaussian with $\hat{\sigma}$:
$$
\hat{\sigma} = \sqrt{8\sigma^2} = 2\sqrt{2} \sigma
$$


Transfer Function from input to node From the IIR filter flow graph below i don't understand how the transfer function is calculated in every node:


Zero Forcing Equalizer in ISI channel I'm studying on designing 9-tap zero forcing equalizer to remove ISI in the channel. I assumed that the channel response is $h[n] = [0.41, 0.815, 0.41]$ and I tried to calculate bit error rate. When I choose $h[n]=[0.21, 0.815, 0.21]$ for example, I get good results. BER goes to 1e-4 for high SNR values. But the code below gives bit error rate of 0.5 when the channel response is $[0.41, 0.815, 0.41]$. Can anybody help me about this problem? I did something wrong but I couldn't find it.
N  = 10^4; % number of bits or symbols
Eb_N0_dB = [0:2:10]; % multiple Eb/N0 values

for ii = 1:length(Eb_N0_dB)

    % Transmitter
    ip = rand(1,N)>0.5; % generating 0,1 with equal probability
    s = 2*ip-1; % BPSK modulation 0 -> -1; 1 -> 1

    % Channel model, multipath channel
    ht = [0.41 0.815 0.41]

    chanOut = conv(s,ht);
    n = 1/sqrt(2)*[randn(1,N+length(ht)-1) + j*randn(1,N+length(ht)-1)]; % white gaussian noise, 0dB variance

    % Noise addition
    y = chanOut + 10^(-Eb_N0_dB(ii)/20)*n; % additive white gaussian noise

    kk=4; % 9-tap filter
    L  = length(ht);
    hM = toeplitz([ht([2:end]) zeros(1,2*kk+1-L+1)], [ ht([2:-1:1]) zeros(1,2*kk+1-L+1) ]);
    d  = zeros(1,2*kk+1);
    d(kk+1) = 1;
    c  = [inv(hM)*d.'].';

    % matched filter
    yFilt = conv(y,c);
    yFilt = yFilt(kk+2:end);
    yFilt = conv(yFilt,ones(1,1)); % convolution
    ySamp = yFilt(1:1:N);  % sampling at time T

    % receiver - hard decision decoding
    ipHat = real(ySamp)>0;

    %     % counting the errors
    nErr(1,ii) = size(find([ip-ipHat]),2);
end

simBer = nErr/N; % simulated ber

Your channel has 3 taps, your EQ has 9 taps which means that your cascaded response has 11 taps. You only force the middle 9 taps of the response to be a unit impulse. If you convolve the channel with your EQ, you see that the middle is indeed a unit impulse but the first and the last sample are very non-zero.

The actual inverse your channel would be an IIR filter. To approximate this with an FIR filter, you need to have a fairly well behaved channel and spent a fair amount of FIR taps. Your channel is NOT well behaved since it got zeros on the unit circle, which means it's not invertible.


Is SSIM a measure of mutual information? As i have seen contradictory information online: is the quality measure SSIM between two images (one base-line, one distorted) directly correlated to the mutual information between the two. One one hand, MSSIM can be partially interpreted as normalized mean square errors between the two on the other hand there are papers online that explicitely extend MSSIM by a mutual information term, suggesting they are not the same. I know this is not very in-depth, but for now i'm look for a simple yes/no -answer so I can decide in which direction to extend my further research.
Thanks!

...is the quality measure SSIM between two images (one base-line, one distorted) directly correlated to the mutual information between the two [?]

The short answer to this is "yes" but it tells you nothing. Because the question is a little bit ill posed.

You have two signals, let them be signals or images, it doesn't matter, you have two observations. And you want to assess how much they are alike.

What does alike mean? You might have an idea about this already, if it is sound do they sound similar? If it is images, do they look similar?

So, things look alike when they co-vary...Sound pressure over time covaries between two sound clips, brightness (reflected light flux) over space covaries between two images.

One of the first metrics of assessing how much alike two signals (images) are that someone might encounter is cross correlation .

If you employ (normalised) cross correlation, on your baseline image (i.e calculate its autocorrelation), you will observe the highest value at the centre of the spatial cross correlation. A value of 1.0.

Now, notice here, we decoupled the concept of co-variance from the way by which we assess it. Co-variance is this notion that things might be going up and down at the same time/place. Cross correlation is a quantification of that notion.

And, there are others. Cross correlation, mutual information, transfer entropy, granger causality, partial correlation, ... and many others.

ALL of these are correlated with each other. Their ranges, sign and even interpretation might be different, but if you adjust for all of these, in general, all of these metrics go in one direction when the observations co-vary and in another when they don't co-vary.

Let's turn to the SSIM, the SSIM is the product of three factors: luminance, contrast and structure with the following expressions:

(Image captured on wikipedia)

Where the $\mu$'s denote means and the $\sigma$'s denote variances with the $\sigma_{xy}$ in particular denoting covariance. By the way, you might want to notice the similarity of the $l,c,s$ components with the covariance and even the correlation coefficient.

Luminance doesn't even count, if I remove the mean from two images, the luminance component goes to zero. It is a little bit more difficult to zero the Structure but I can do it, if the two images do not covary with each other, a condition that is met if the two images are pure white noise. This will leave the contrast with some value that would have to be equal to the separate variances.

Long story short: SSIM contains a covariance term, covariance and mutual information are quantifications of the same concept, therefore, SSIM and mutual information are correlated.

BUT, this tells you nothing. The point is how are these two correlated and the interpretation. So, there is a different mapping between cross correlation and mutual information because mutual information has that $\log$ factor in it and, if it is mutual information then would it assume memory of some length?

So now we are on the next level, WHAT are you trying to achieve? To then move on to WHICH aspect of your signal do you need to quantify (and how) in order to get there?

Hope this helps.


Matlab code for Simulating of Phase Noise for OFDM systems I am trying to simulate phase noise for OFDM systems.
I have been reading about phase noise, some literature suggests that phase noise is cumulative as follows:

$\phi(n)=\phi(n-1)+\omega(n)$

where $\omega(n)$ is a Gaussian random variable.
The received OFDM signal is expressed as:

$r(n)=x(n)e^{j\phi(n)}$
where $x(n)$ is the transmitted OFDM signal.

I am looking for Matlab code for phase noise simulation.
thanks

Below is a short code to do this for 1024 samples. "phi" array below is the required phase for each samples

N=1024; %number of samples

phi=zeros(1,N);
phi(1,1) = pi/4;. %initial phase, could be set to zero if the phase is set to zero by the front end processing

Phase_snr=10; %snr in phase estimation,equivalent to 1db

for n=1:N-1
phi(1,n+1)=phi(1,n) + (1/sqrt(phase_snr))*randn(1);
end


Minimize the Cost Function of Values of Vectors Based on Their Amplitude I have  two vectors $X = [x_1,x_2,x_3,x_4]$;  and $Y = [y_1,y_2,y_3,y_4]$; I know that $|x_1|$ = $|y_1|$, and $|x_2|$ = $|y_2|$,... so on. it means the difference is only in the sign. it might be either similar to same sign.
Those two vectors have been affected by such noise, and I need to create a cost function to recover them. Is it possible to perform that using such algorithm of optimization ?

Since there is no prior at the Vector level this is basically element wise problem.
Moreover, if we assume the noise to be White Noise with zero mean then the answer can be very simple.
Since the phase difference is always a multiply of 180 [Deg] we can, without loss of generality, assume they are on the real axis.
So what we have can be modeled as:
$$ {z}_{1} = {w}_{1} + {n}_{1}, \; {z}_{2} = {w}_{2} + {n}_{2}, \; \left| {w}_{1} \right| = \left| {w}_{2} \right| $$
So the answer will be $ \hat{w}_{1} = \operatorname{sign} \left( {z}_{1} \right) \frac{\left| {z}_{1} \right| + \left| {z}_{2} \right|}{2} $ and $ \hat{w}_{2} = \operatorname{sign} \left( {z}_{2} \right) \frac{\left| {z}_{1} \right| + \left| {z}_{2} \right|}{2} $.
This comes from a simple intuition. Let's assume bot are positive and $ {w}_{1} = {w}_{2} $. Then of course the answer will be, for any white noise, to average the results. So we are doing the same, averaging them, just taking care of the case they have opposite phases (In the complex plane they are on the same circle just could be on opposite directions).
If you assume AWGN and build the Maximum Likelihood Estimator you'll get the same result.


How does OpenCV find chessboard corners? What are the internal workings of the OpenCV function findChessboardCorners?

The source code of OpenCV is available, so I'd recommend just taking some time going through the code. The relevant file for this particular function is:

calib3d/src/calibinit.cpp

I've not looked into it in detail, but it looks like

CV_IMPL
int cvFindChessboardCorners( const void* arr, CvSize pattern_size,
                             CvPoint2D32f* out_corners, int* out_corner_count,
                             int flags )

is the main implementation of this method. In here they

Use cvCheckChessboard to determine if a chessboard is in the image
Convert to binary (B&W) and dilate to split the corners apart
Use icvGenerateQuads to find the squares.

The code then seems to go though a set of checks to condense these quads to chessboard corners, including icvFindConnectedQuads, icvCleanFoundConnectedQuads to remove extra corners, icvCheckQuadGroup, and icvCheckBoardMonotony.

All of these functions are implemented in the same file, apart from cvCheckChessboard which is in calib3d/src/checkchessboard.cpp. Depending how well you wanted to understand the code, there seem to be a number of debugging lines, which can be included if you #define DEBUG_CHESSBOARD, that may help you see what's going on.


Reflecting negative frequency to positive frequency I am trying to synthesize sinusoids by using window function in the frequency domain.

It involves:

In frequency domain, shift the window to center around the peak frequency
To generate a DFT frame, sample a few values of the window around the peak as the spectral motif
Inverse-Fourier transform the spectral motif, to generate the sinusoid in the desired frequency

This approach works great in general, except for synthesising low frequencies. Because when shifting the window to the low frequency, the left side of the window will sit in the negative frequency domain.

The figure in the middle demonstrates the issue (T(k) sits in the negative domain)

I found a solution here, it suggested to add the complex conjugate value of the left tail of the window (in the negative domain) to the DFT bin that’s on the right tail of the window (in the positive domain). Which I can’t make sense of, and following this solution creates even more distortion. So I wonder if anyone knows how to do it properly. Any suggestions would be very appreciated!

Some excerpt from the aforementioned solution:

  The reflection about the k=0 axis is due to the specific embodiment
  described herein for synthesizing a sinusoid. For each real sinusoid,
  one peak exists in the positive frequency bins and another peak exists
  in the negative frequency bins. In the embodiment wherein only the
  peak in the positive frequency bins is synthesized, a peak centered
  about a low positive frequency bin spills into the negative
  frequencies (as shown by the plot for Ht(k−bc) in FIG. 3). Similarly,
  a peak centered about a low negative frequency bin spills into the
  positive frequencies. The portion of Ht(k−bc) in the negative
  frequencies that is reflected, or T*(−k), represents the portion of
  the peak centered about the negative frequency bin that spills into
  the positive frequencies.

PS. Some time ago I've raised this question on DSP-related forum. I've got very detailed suggestion from Robert, to ignore the bins lying on negative domain to specify the whole positive frequency, and then complex conjugate to reflect it to the negative frequency, which has improved the problem but it still can’t go down below 80 Hz. So I thought I’ll post again here.

You could use another approach: allow negative frequencies, do a complex IFFT, discard the imaginary parts of the time domain samples, and multiply the result by 2.
Let's try it in Octave (MATLAB clone) with a motif $[1 + 2i, 3 + 4i, 5 + 6i, 7 + 8i]$ shifted so that its leftmost bin ($1 + 2i$) lands on a negative frequency. FFT length is 8. (I rewrote the results for readability.)

> g = [3 + 4i, 5 + 6i, 7 + 8i, 0, 0, 0, 0, 1 + 2i]
g =

   3 + 4i   5 + 6i   7 + 8i   0 + 0i   0 + 0i   0 + 0i   0 + 0i   1 + 2i

> y = real(ifft(g))*2
y =

   4.00000  -0.89645  -2.00000   0.98223   1.00000  -1.60355   0.00000   4.51777

> fft(y)
ans =

   6 + 0i   6 + 4i   7 + 8i   0 + 0i   0 + 0i   0 + 0i   7 - 8i   6 - 4i

The same time-domain signal $y$ can be generated by synthesizing the first 5 bins of $\text{fft}(y)$ and feeding them to real IFFT (which I assume you use). The 3 remaining negative-frequency bins the IFFT kind of generates internally. You can follow these steps for the synthesis:

1) For each bin of the shifted motif that lands on a positive-frequency bin, copy the value verbatim. (You are already doing this)

2) For a bin of the shifted motif that lands on the 0 Hz bin, double its real part and write that to the 0 Hz bin, discarding the imaginary part, like $3 + 4i \rightarrow 6 + 0i$ in the example. (You are already doing this)

3) For each bin of the shifted motif that lands on a negative-frequency bin, add its complex conjugate (same real part, sign of imaginary part flipped) to the positive-frequency bin that is at the same distance from the 0 Hz bin. This will be one of the bins already written to in step 1. In the second bin of the example $5+6i$ and $\text{conj}(1 + 2i) = 1 - 2i$ were summed to get $6 + 4i$. Don't expect the imaginary parts to cancel, because even if the motif might be symmetrical, those two motif bins are not located symmetrically on the opposite sides of the motif, unless you are synthesizing a 0 Hz sinusoid. It is perfectly fine if the imaginary part remains non-zero; frequency domain bins can have complex values and still yield real-valued time-domain data.

Oh, and if the shifted motif reaches the Nyquist bin (corresponding to half the sampling frequency) or beyond, do exactly the same things there, mirroring around the Nyquist bin.


Magic of twiddle factor in DFT In DFT calculation following formula is used
$$X[k] = \sum_{n=0}^{N-1} x[n] e^{-j2\pi\frac{kn}{N}}$$
where the Twiddle factor is known as complex root of unity, that is a complex number $W_N$ such that $(W_N)^N = 1$. There are $N$, $N$-th roots of unity.
$$\big( W_N \big)^{k} = e^{j2\pi\frac{k}{N}} = (\omega_N)^k$$

Suppose I am doing 8 point DFT. Then I am able to get the 8 frequency bins output of this DFT formula : $X[0], \ldots, X[7]$.

Please suggest what magic this twiddle factor does that when each twiddle factor is multiplied with samples & summation is done then we get the respective frequency bin?

Please see Different way to separate a particular frequency from a signal where I explain the difference and similarity between homodyning the signal and homodyning the coefficients. The DFT is an example of homodyning the coefficients, such that multiplying the coefficients of what would otherwise be a moving average low pass filter by $e^{j2\pi f_c t}$ creates a bandpass response around frequency $f_c$. Twiddle factors are simply the discrete samples of $e^{j2\pi f_c t}$!

Observe specifically in the DFT matrix that the second row steps through all the samples in the roots of unity, so goes around the unit circle once (representing the normalized frequency of one sample per cycle).  The third row steps through every other sample in the roots of unity and in the process goes around the unit circle twice, the fourth row goes around the unit circle
three times, etc. (The first row is all ones representing no rotation, DC, or the base moving average filter). Note that going past $f_s/2$, where $f_s$ is the sampling rate, that the higher positive rotations alias to negative rotations. For example in the 4 point DFT as illustrated below, the normalized frequencies of 0, 1, 2 and 3 are also 0, 1, -2, -1; rotating twice around the unit circle counter-clockwise is the same as rotating twice around the unit circle clockwise when the twiddle factors are just 0, j, -1 and -j: (0, -1, 0, -1). Rotating three times counter-clockwise is the same as rotating once clockwise: (0, -j, -1, j).

The DFT as a filter bank is much clearer if you consider the case of a "streaming DFT" where we shift through a much longer sequence and compute the DFT once after each shift. In this case every point of the DFT would be identical to the output of a FIR filter with the coefficients as given by the DFT row associates with that output.


bandwidth between the signal and the channel I have a signal with a bandwidth $B=4$ kHz and the bandwitdh is between $6$ kHz and $10$ kHz. I want to transmit this signal in a channel with bandwidth $B'=4$ kHz but for the channel the bandwidth is between $1$ kHz and $5$ kHz.

This is my question:
Can we pass a signal in this channel even if the bandwidth support is not the same?

If I understand correctly you would like to pass a bandpass signal of 4Khz bandwidth in the frequency range 6 to 10 Khz, but you would like to transmit this directly over a channel which has good channel gain in 1 to 5Khz. Clearly, the signal when passed over a channel with this frequency response will suffer complete degradation.

Hence, we cannot transmit this signal directly over this channel. You could first downcovert the signal by 5 Khz and then transmit.


Hermitian FFT vs discrete Fourier transform. Python's $\tt hfft$ vs $\tt fft$ What is the difference between the Hermitian FFT and DFT? Particularly in Python, there are two functions fft and hfft.

numpy.fft.hfft(signal) vs numpy.fft.fft(signal)

What I simply could find out is:
The Hermitian has to do something with symmetry and needs 50 times longer to calculate, while producing a 'slightly' different result than the 'discrete' FFT. (tested on an audio file of machinery sounds and length of 1.5 sec).

You have to understand that fft is the general function, which always works. hfft, fftr and their derivatives are optimized for special signal constellations and should either be faster or more accurate (if this is at all necessary with modern floating point accuracy).

The hfft function is used, when you expect the spectrum to be real-valued. This happens, if your input signal has even symmetry. I doubt that your audio signal does have this property. Have you been looking for fftr (i.e. fft for a real-valued signal)?

In case your input signal has the even symmetry property, the output of fft and hfft should be exactly the same (up to numerical residual errors). Also, hfft should be a bit faster in this case (otherwise, it would make no sense to use it).


Estimation of Instantaneous Amplitude I'm reading a paper on EMG analysis. The formulas are all clear to me, but the paper refers to the signal amplitude as "instantaneous" amplitude.

I know what instantaneous mean, but what does it mean in the context of signal processing?

Instantaneous amplitude is usually done in conjunction with analytic signals and the Hilbert transform. See this answer on this site.

  For signal $s(t)$, given its Hilbert Transform $\hat{s}(t)$ it is defined as a composition:

  $$s_A(t)=s(t)+j\hat{s}(t) $$

  The Analytic Signal that we obtain is complex valued, therefore we can express it in exponential notation:

  $$s_A(t)=A(t)e^{j\psi(t)}$$

  where:

  $A(t)$ is the instantaneous amplitude (envelope)

  $\psi(t)$ is the instantaneous phase.


Advantage of using nonideal filter in purely digital application Say I'd like to bandpass filter a signal with corner frequencies 1 kHz and 2 kHz (i.e., passes signals between 1 kHz and 2 kHz and suppresses other frequencies) using software like MATLAB. Intuition tells me there's two ways to do this:

Use filter design tools to create, say, a 10th order Butterworth filter. Use the filter command with the filter coefficients and input signal. This works, but the filter is nonideal (specifically, there will be some passband ripple, and the frequency cutoffs won't be immediate -- a signal at say, 2.1 kHz, will still be partially passed).
Perform an FFT on the input signal. Set any FFT element corresponding to a frequency outside the 1 kHz-2 kHz range to zero. Inverse FFT the input signal. This will essentially be an ideal filter -- anything between 1kHz-2 kHz is passed with unity gain (no passband ripple), and anything outside 1  kHz-2 KHz is exactly zero (immediate cutoff).

Intuition tells me that #2 is better, because it's an ideal filter. I've learned that an ideal filter can't be built using analog components, but if a certain application is purely digital, why not use the ideal filter?

That being said, I've seen lots of code that uses tools like a nonideal 10th order Butterworth filter -- so what am I overlooking? What's the advantage of using a nonideal filter?

Your assumption #2 is false.  Zeroing bins does not produce an ideal filter with unity gain in the passband.  Far from it if you look closely (between FFT bins).


Calculate the signal's average power We have this signal:

$$x(t)= 2 + 2\cos(2\pi f_0 t) + \frac{2}{T_0} \operatorname{sinc}\left(\frac{2 t}{T_0}\right)e^{j2 \pi 4/ T_0} + \operatorname{sinc}\left(\frac{2 t}{T_0}\right)e^{-j2 \pi 4/ T_0}.$$

I have calculate Fourier trasformed and have found:
$$
X(f)=\delta(f)+ \delta(f-f_0) + \delta(f+f_0) + \operatorname{rect}\left(\frac{f-\frac 4{T_0}}{\frac 2{T_0}}\right) + \operatorname{rect}\left(\frac{f+\frac 4{T_0}}{\frac 2{T_0}})\right)
$$

I must find the average power. If use Parseval's have:
\begin{align}
P_x&=\frac{1}{T_0} \int_{-\infty}^{\infty}\lvert X(f)\rvert^2df \\
&= \frac 1{T_0}\left(\int\lvert\delta(f)+\delta(f+f_0)+\delta(f−f_0)+\operatorname{rect}+\operatorname{rect}\rvert^2df\right).
\end{align}
How can I solve this?

A signal either has finite energy, finite power or even infinite power. If it has finite energy, it will have zero average power, according to your definition
$$P_x=\lim_{T_0\rightarrow\infty}\frac{1}{T_0}\int_{-\frac{T_0}{2}}^{\frac{T_0}{2}}|x(t)|^2dt.$$

Knowing that the sincs are orthogonal to each other, as well as the cos function, the average power is given by

$$P_x=P(2)+P(2\cos(2\pi f_0t))+P(\operatorname{sinc})$$

where $P(...)$ is the power of the components in brackets. The sinc function has finite energy, hence it has zero power on average. Hence, we have

$$P_x=P(2)+P(2\cos(2\pi f_0 t)$$
with
$$P(2)=4$$
and $$P(2\cos(2\pi f_0 t)=2^2\cdot\frac{1}{2}=2.$$

Therefore, the overall power is given by

$$P_x=6.$$

For the calculation of the $\cos$, consider
$$P_{\cos}=\frac{1}{T}\int_{0}^{T}\cos^2(2\pi t/T)dt=\frac{1}{2}$$


wiener filter (non-causal) calculating covariance matrices I'm learning about the Wiener filter, and I'm working on my own implementation. I'm starting out with the non-causal filter, and I need to calculate some covariance matrices. If the signal is assumed to be a linear combination of a target signal and noise, then the signal can be modeled as: $$x[n] = y[n] + w[n]$$ where x is the observed signal, y is the target and w is noise. If I understand correctly, the weiner filter is a linear-minimum-mean-square-error estimation process. What I'm getting hung up on is that we need to compute a matrix of expected values based on the input signal x, namely $$E[xx^{T}]$$ and I know that to do this we need to apply the expectation operator to each element in the matrix $$xx^{T}$$ but I don't know how to do that computation, or what it even means.

The expectation operator can not be applied in real.In Weiner filter you will require to find the value of R i.e. correlation matrix of input signal which is ideally E[X(n)X'(n)] and also P i.e  cross correlation matrix of input signal and desired signal E[X(n)d(n)]. Actual meaning of this expectation is that the correlation is obtained by using a large number of realisation of the input signal called ensemble .But in reality we have only one realisation and that is the input signal which you are going to apply to the weiner filter to calculate output. hence its exact value cant be find out. and therefore to use weiner filter in reality we require to apply various approximations to calculations of covariance matrix which finally result in various adaptive algorithms .As an example ,LMS algorithm which uses the following approximations for above matrices,
               R=X(n)*X'(n)   and P= X(n)d(n)
Now these values can be calculated using the only available realisation of signal.Although these Algorithms do not provide full convergence to weiner solution but approximates weights close to weiner solution. This error in solution is the result of approximation used. Weiner solution is actually the ideal solution of a adaptive filter hence it can not be attained because the Expectation operation can be applied in theory but practical implementation is impossible.
LMS,NLMS,BLMS,RLS these are some other algorithms that can also be used .each have its own good side and a bad side and have to be applied according to necessity.


Overlapping with flat-top window MATLAB I am applying a flat-top window to signal because I am very interested in its actual amplitude in the time-domain. I have read that overlapping will improve results. I am using MATLAB and the following is the code that I am using to compute FFT with flat-top window:

time = 0:0.1:500;
data = 1.8*sin(2*pi*0.2003*time);Fs = 10;
N = length(data);
w = flattopwin(N);
wdata = data(:).*w;
ws = sum(w);
Y = fft(wdata)/ws;
freq = 0:Fs/N:Fs/2-Fs/N;
freq = freq';
amplitude = 2*abs(Y(1:floor(N/2)));
plot(freq,amplitude);

Can someone please help me to modify the code so that I can include overlapping? I am relatively new in this field and am still confused. Thanks.

Higher size of your FFT best the resolution of frequencies, but of course it has a cost resulting in a longer to process !

To apply overlap you need define one window size to analyse and one hop/step to walk around your data!

This picture illustrates how this works:

Here a simple code using Window size = 2048pts and 50% overlap step!

n_samples = length(yourdata);
winSize = 2048;

%50% overlap or 0 to not use overlap
OverlapStep = 50;

if OverlapStep > 0

    Overlap = floor((OverlapStep*winSize) / 100);
    nFrames=floor(n_samples/Overlap)-1;
else
    Overlap= winSize;
    nFrames=floor(n_samples/Overlap)-1;
end

k=1;
while ( (k+winSize-1) <= n_samples )

    FrameSignal = yourdata(k:k+winSize-1);

    %%%%%%%%%% Start process FFT here  %%%%%%%%%%%%%%%%%
    w = flattopwin(winSize);
    wdata = FrameSignal(:).*w;
    Y = fft(wdata)
    .
    .
    .

   k=k+Overlap;
end


Non-zero mean of Kalman innovation I am using a Kalman filter to fuse gyro and inclinometer data. The prediction step is given by:

$\hat{\alpha_{i}}^- = \hat{\alpha}^+_{i-1} + \omega_{i}$

Where $\hat{\alpha_{i}}^-$ is the prediction for the angle, $\hat{\alpha}^+_{i-1}$ is the previous a-posteriori state estimate and $\omega_i$ is the change in angle measured by the gyro.
The innovation of the filter is given by:

$\alpha_i - \hat{\alpha_{i}}^- $

where $\alpha_i$ is the current measurement of the inclinometer. As far as my understanding of the Kalman filter goes, the innovation is expected to have a zero mean. However, the following Figure shows the histogramm for some 25000 datapoints of the innovation:

The mean is obviously non-zero. Would it be correct to claim that this means that some error source is present, for example an uncorrected gyro offset? Or is the non-zero mean more an "academic" property of the filter, can't be expected to be true in a real-world implementation, and thus a non-zero mean should not be overrated?

Mean of innovation sequence is zero if the Kalman filter is optimal. A non zero mean implies that there are certain errors present in your assumptions of covariance of noise or even your system matrices(if you are using a state space model).


Image averaging for impulse noise I understand the mathematics behind the use of image averaging for gaussian noise. Why can't it be used for impulse noise ? Would it be because there is no way to prove that impulse noise is additive ? And that it has no mean ?

Also, am I right to say that impulse noise does not affect every pixel in an image, as opposed to gaussian noise ?

Thanks !

I'm not in the image processing filed, but I would say:
Yes, Gaussian noise means that each of your pixel is assumed to be produced by a stochastic process. Hence you can for instance think of it as a "real" intensity corrupted by an additive normally distributed white noise. Normally distributed means also that if you take $n$ observations of the same pixel you most likely will end up with $n$ wrong measurements. I think this is actually the worst kind of noise since given any pixel, then any other pixel inside a given neighbourhood will be corrupted, hence telling which observation is to trust more is difficult. The way I would go is probably (low pass) filtering.

You can think to the impulsive noise again as an additive noise potentially present in all the pixels, but potentially also localized in only a subset of pixels or observations (either it is or it is not). I personally would not go for averaging or low pass filtering since it will never erase completely the noise.
Assume indeed to have $n$ images of the same scene. Assume that for a given pixel $p$ whose "actual" value is $p=0$ you measure it as $255$ in one image and as $0$ in $n$ images. If you take the mean then you will chose $p = 255/n$ as the actual value. Infact a low pass filter is a process with infinite impulse response and it tends to zero only asymptotically.
In such case for me it looks more effective to take the median, over the $n$ samples. With the median you obtain that you estimate of the actual values of $p$ is exactly $0$.
You can also take the most-frequent value (the mode) and it will perform better than the mean. In general you can also use a median or a mode inside the same image by taking them in a neighbourhood of a given pixel. Infact for example, take a $10\times 10$ sub-image and assume it corresponds to a white background. Assume that 30 over the 100 pixels are black due to the impulsive noise, then you can substitute to each pixel the median over the entire sub-image obtaining a completely white patch. A mean will lead to a 30% grey (76 / 255)


how twiddles are calculated in the code snippet below from ti.com? I am trying to port a C code for fft/ifft to a DSP (TMDSEVM6678LE) from ti.com. In one part of the code it tries to calculate the twiddle coefficients as you can see below; however, I can't understand why the matlab code is so different ...

I would appreciate any comments.

void tw_gen_cn (double *w, int n)
{
    int i, j, k;
    const double PI = 3.141592654;

    for (j = 1, k = 0; j <= n >> 2; j = j << 2)
    {
        for (i = 0; i < n >> 2; i += j)
        {
            w[k]     = cos (2 * PI * i / n);
            w[k + 1] = sin (2 * PI * i / n);
            w[k + 2] = cos (4 * PI * i / n);
            w[k + 3] = sin (4 * PI * i / n);
            w[k + 4] = cos (6 * PI * i / n);
            w[k + 5] = sin (6 * PI * i / n);
            k += 6;
        }
    }
}

The matlab code for twiddle coeffs. are as below:

% Program for generating n-length FFT's twiddle factor
% By: Denny Hermawanto
% Puslit Metrologi LIPI, INDONESIA
% Copyright 2015
function [real_twiddle, im_twiddle] = tw_gen_cn(fft_length)
    if nargin < 1
        [real_twiddle, im_twiddle] = tw_gen_cn(8);
        return;
    end
% fft_length = input('Enter FFT length:');
for mm = 0:1:(fft_length-1)
    theta = (-2*pi*mm*1/fft_length);
%   Twiddle factor equation
%   twiddle = exp(1i*theta);
%   Euler equation for complex exponential
%   e^(j*theta) = cos(theta) + j(sin(theta))

    twiddle(mm+1) = cos(theta) + (1i*(sin(theta)));
    real_twiddle = real(twiddle);
    real_twiddle = real_twiddle';
    im_twiddle = imag(twiddle);
    im_twiddle = im_twiddle';

end;
end

These are similar and provide the same result. In the c code every other sample:  w[k+1], w[k+3], etc, are the imaginary twiddle factors in the MATLAB code:

For instance the lines in C:

w[k]     = cos (2 * PI * i / n);
w[k + 1] = sin (2 * PI * i / n);

Corresponds to this line in MATLAB :

twiddle(mm+1) = cos(theta) + (1i*(sin(theta)));

Where "theta" in MATLAB is "2*PI*i/n" in c.


Can CV or DS be seriously considered DSP roles? This is not a technical question but can have a yes/no answer, and a why statement. My understanding about digital signal processing (DSP) is that it allows to take real-world analog signals (1D, 2D, etc...), sample them to create digital representations on the computer, and either way: extract useful information from filters or modify the signal to take it back to the analog world (using a digital-to-analog converter if you like). There is an overlap with other disciplines like pattern classification or statistics as DSP implements dimensionality reduction techniques that make easier to classify samples; or with embedded systems, that have DSP ICs. I am not considering analog signal processing which is another topic for a question.

There a talk explaining that there is an "identity crisis" in DSP nowadays with all the machines learning (ML) methods because they can process raw data without having to apply the carefully crafted DSP knowledge (let's say, using a CNN instead of extracting feature vectors). Of course, sacrificing computational power, time tunning hyperparameters and or the effort to create enough training data. But people don't care as long as they get the processing they want and don't have CPU, power, or memory constraints. This also reminds me this article pointing that DSP roles are not longer the same they used to be 25 years ago.

I have observed certain similarity between computer vision (CV) and data science roles with DSP, but I am still trying to convince my self that someone can really make some DSP working in them. On a CV role you are in charge of dealing with 2D+ signals and perform tasks such as face recognition or pose estimation (which has historically been a DSP task) but the problem is that the majority of the work is to reuse components-out-of-the-shelf (COTS), like openCV, to specific areas, without having to demonstrate full understanding of the yielding techniques (of course, I am talking about someone that hasn't to work in a lab where some research work needs to be done). And about data science, I know that they use a lot of ML, data visualization, statistics, business logic, A/B testing; and have to use the data gathered by data engineers, but they do not use DSP-core techniques like sampling and interpolation or FIR/IIR/adaptive filtering. A former data scientist said that you can even not use ML at all if you want.

Can someone really do DSP with these positions?

P.S. I know there are a lot of roles that have a better overlap with DSP (e.g. EEG, acoustics, radar, music information retrieval, ASR, media codecs) but I am focusing on the ones I mentioned, mostly because they seem to be doing a big impact today. And I believe someone would answer something like "it depends on how good you want to be at your job" but I'd like to consider roles where you are naturally forced to build DSP ad-hoc solutions.

in my opinion

your question seems based on a combination of thoughtful reflection and possibly, pre Alzheimer’s cognitive biases. ;)

There was an article in the WSJ not too long ago , I’m outside the paywall so I can’t provide a link, that looked at IBM’s efforts using Watson at major Cancer research centers.

The conclusion was that Watson wasn’t offering much in terms of originality. It wasn’t revealing those deep unifying subtle links in the masses of data it was being feed. It tended to find what was already known.

hard problems are still hard particularly if the way the data is represented ( labeled) isn’t relevant to what might be hidden. The threshold of what might be considered for brute force search might be advancing but NP hard hasn’t gone away.

I don’t believe I answered your question but I don’t think your question has a clear answer.

When I was an undergrad, the EE department had just the year before my junior semester had dropped land surveying as a core requirement. skill mixes change.  not many years ago i was calculating viewsheds to estimate RF coverage. Land Surveying would have probably had some use.

Are patterns truth?


Bandpass filter that automatically adapts its bandwidth when a transient is detected (to avoid to smoothen the transient) Let's say we want to isolate a band 1000 hz +/- 50 hz.

Obviously, limiting the bandwidth by applying a passband filter will always destroy a bit the sharp transients (a Dirac or a rectangular envelope / Heavyside step function requires all frequencies, so if limit the bandwidth, we lose a part of it, and it becomes smoother).

Question: are there some adaptative band-pass filters that would auto-extend their bandwidth for a short time when a transient is detected, in order to not lose the sharp transients?

In this example (1000 hz sinusoid modulated by a rectangle envelope, input in blue):

the filter would of course still focus on the 1000hz +/- 50hz band, but it would extend its bandwidth near the transient so that the transient is not smoothened like with a normal filter (signal in red).

Does such an adaptative bandpass exist, and is it available easily in most languages (Matlab, Python, etc.)?

NB: on this graph there is nothing else except the 1000hz sinusoid, so you may wonder "why bandpass filtering?", but it's just an example, in the general case, it would be a broadband signal.

Ok, I think I got your point. You want a BPF, $H(z)$, that auto extends its bandwidth accordingly to the energy distribution in the magnitude spectrum. If you have a pure 1k Hz sinusoidal tone (that corresponds, in the frequency domain, to a dirac delta located at $\omega_0=\pm2\pi 1$k rad/s), you want to pass only frequencies in the 1k$\pm 50$ Hz range, and if you have a transient event with a white noise-like distribution, you want an all-pass filter to preserve the sharp attack.

What you need is a resonator filter [1]:
$$
H(z)=\frac{(1-\lambda)\sqrt{1+\lambda^2-2\lambda\cos(2\omega_0)}}{1-(2\lambda\cos(\omega_0))z^{-1}+\lambda^2 z^{-2}},
$$

its behavior for different values of $\lambda \in [0,1]$ is like so:

so for $\lambda\to 0$ you will get a flat response to catch transient events, and for $\lambda\to 1$ you will have a localized filter at the desired frequency. Here, for illustration purposes, I set $w_0=\pi/2$, but you can change the desired frequency using the formula $w_0=2\pi F_0/F_s$.

For setting $\lambda$ automatically you can use the spectral flatness estimator [2]:

$$
f = \frac{\left(\prod_{n=0}^{N-1}{x[n]}\right)^{1/N}}{\frac{1}{N}\sum_{n=0}^{N-1}{x[n]}},
$$

which is $f=1$, when the magnitude spectrum is completely flat, and $f=0$, when the magnitude spectrum is completely localized. Therefore, you can make $\lambda=1-f$. I wrote the following code to exemplify how you can apply this control:

Fs=16e3;
F0=1e3;
w0 = 2*pi*F0/Fs;
x1 = [zeros(1,50),2*rand(1,50)-1];
x2 = 0.7*sin(w0.*[1:100])+0.3*rand(1,100);
x3 = 0.7*sin(3.5*w0.*[1:100])+0.3*rand(1,100);
y = [adaptiveResonatorFilter(x1,w0), adaptiveResonatorFilter(x2,w0), adaptiveResonatorFilter(x3,w0)];
plot([x1,x2,x3],'linewidth',2)
hold on
plot(y,'linewidth',2)
xlabel('Samples')
ylabel('Amplitude')
legend('Original','Filtered')

function y = adaptiveResonatorFilter(x,w0)
    X = fft(x);
    mX = abs(X);
    mX = mX/max(mX);
    sf = mean(mX,'g')/mean(mX,'a')
    lambda = ifelse(0.5<1-sf, 0.99, 0.0)
    B = (1-lambda)*sqrt(1+lambda^2-2*lambda*cos(2*w0));
    A = [1,-2*lambda*cos(w0), lambda^2];
    [H,W] = freqz(B,A,linspace(-pi,pi,length(mX)));
    Y = X .* fftshift(H);
    y = real(ifft(Y));
end

which gives the following output:

where you can see that transient part is kept untouched, the 1k Hz pure tone contaminated with noise has been cleared and the 3.5k Hz pure tone has been attenuated, as you wanted.

Note: I am taking this as the definition of "transient attack". Please correct me if I misunderstood.

M. Vetterli, P. Prandoni. Signal Processing for Communications. EPFL press.
https://en.wikipedia.org/wiki/Spectral_flatness


Impulse response with sweep from measurement Does someone have a matlab script to get the impulse response from the recorded measurements using log sweep? I have the input file (sweep) and the output and I need to get the impulse response and then calculate reverberation times for octave bands.

referenceFilename = 'reference.wav';
responseFilename = 'response.wav';

[reference, referenceFs] = audioread(referenceFilename);
[nReferenceSamples, nReferenceChannels] = size(reference);

[response, responseFs] = audioread(responseFilename);
[nResponseSamples, nResponseChannels] = size(response);

if referenceFs ~= responseFs
  error('reference audio and response audio files must have same sample rate');
end

nFFT = 2^( ceil( log2(nReferenceSamples+nResponseSamples) ) );

x = zeros(nFFT, 1);
y = zeros(nFFT, 1);

x(1:nReferenceSamples) = reference(:,1);    % use only left channel
y(1:nResponseSamples) = response(:,1);      % use only left channel

h = ifft( fft(y) ./ fft(x) );

h = h(1:nFFT/2);                            % truncate latter half

plot(h);

sound(h, responseFs);


Comparison Between Average Kernel (Box Kernel) and Gaussian Kernel In image processing, we have two kinds of major kernels that are average kernel and Gaussian kernel. For image segmentation, which is difference between average kernel and Gaussian kernel? I found some paper said that they are similar, and average kernel implement faster than Gaussian kernel, right?When we use average kernel

If you pre calculate the filter coefficients the complexity of the convolution is set by its radius only.

Yet, if all coefficients are the same, it could be reduced into summation and one normalization multiplication.

Yet, I don't think it will have major effect on modern CPU's.
Certainly not on those who supports fused multiply and addition.


Distinguish between sinusoidal oscillations and rythmic pulses I want to distinguish between a sinusoidal oscillations and rythmic pulses.
Which method can i use to do that?

I want to implement a detector that can distinguish between this two types of signals.

i have calculated the fast fourier transform of this signals. however to determine if it is a rhythmic pulses or sinusoidal oscillations, i have to fixe a threshold manually. It seems not good.

You are going to want to take your DFT (Discrete Fourier Transform) on a whole number of repeat cycles.

Here is a good way to identify boundaries for your repeat pattern:  Use the difference method I describe in my blog article "Exponential Smoothing with a Wrinkle".  Smooth your signal heavily and locate your boundaries by finding positive to negative zero crossings.  Select a duration of at least three cycles, but keep the number limited so you know your cycles have a similar shape within your duration.

Take a DFT of your original signal on your duration.  A FFT is a fast version of this.  For real valued signals, only the first half of the return values are of interest.  The second half are a mirrored complex conjugate of the first and don't add any new information.

For a pure tone, only one bin will have a value in it.  The rest will be zero or near zero.  For a more complicated waveform, the bin values at multiples of the number of cycles you included will give the data to calculate Fourier coefficients for your waveform.

By a "not pure sinusoidal oscillations", I assume you mean a pure tone with a little harmonic distortion.  In this case, the magnitudes of the bins of the harmonic tones will be smaller compared to a sharp "rhythmic pulse" like the ones in your graph.

If you are using MATLAB, be aware that it is one based, and DFTs are zero based so you will have to subtract one from the MATLAB index to get the DFT bin number.

Hope this helps.

Ced

Followup

If your rhythmic signals are like the one you showed, you can simply count zero crossings of your signal within the bounded intervals defined above.


RC high pass filter differential equation I have the question to derive:
$$RC \frac{dy(t)}{dt} + y(t) = RC \frac{dx(t)}{dt}$$
which describes a high pass filter or $CR$ circuit with $x(t)$ the input and $y(t)$ the output.
I am really struggling to derive this
I start with Kirchhoff's voltage law
$$x(t) = Ri(t)+y(t)$$
Use the differential equation for capacitor current. (in series so resistor and capacitor current are the same)
$$i(t)=C \frac{dy(t)}{dt}$$
I then substitute this in to Kirchhoff's voltage law
$$x(t)=RC \frac{dy(t)}{dt} +y(t)$$
This is the same workings out as for a low pass filter (RC circuit). Is this correct?
If so how do I get to the initial equation?

In an RC high-pass filter, the output is taken across the resistor instead of the capacitor. Otherwise, if the output was taken across the capacitor, the RC high-pass filter would start to work as an RC low-pass filter.
An example RC high-pass filter circuit is given below:

The basic idea of using an RC high-pass filter is to filter out the input signal’s low frequency components and let the high frequency components to pass through the succeeding stages of a network.
Actually, this filtering operation depends heavily on the capacitor’s properties for varying frequencies. The formula of a capacitor’s capacitive reactance which represents its resistive nature for varying signal frequencies is given below:
$$X_C{(j\omega)} = \frac{1}{j\omega{C}} \Omega$$
As it can be seen, the formula results in a complex value. According to the formula above; if the frequency $f$ of the signal which is passing through the capacitor decreases, so do the angular frequency $\omega$ because of the formula $\omega = 2\pi{f}$, the magnitude of the capacitive reactance increases. Therefore, basically, the signal which is passing through the capacitor is attenuated.
So, in the light of these knowledge, let’s analyse the given circuit.
If the Kirchhoff’s voltage law is applied on the circuit, we get;
$$-V(t) + V_C{(t)} + V_R{(t)} = 0$$
$$-V(t) + {\frac{1}{C}\int_{t_0}^{t} i(t) \, dt} + V_C{(t_0)} + V_O{(t)} = 0$$
If the assignments $V(t) = x(t)$ and $V_O{(t)} = y(t)$ are done in order to involve classical input and output notations as in the systems, the equation becomes;
$$-x(t) + {\frac{1}{RC}\int_{t_0}^{t} y(t) \, dt} + y(t) = 0$$
$${\frac{1}{RC}\int_{t_0}^{t} y(t) \, dt} + y(t) = x(t)$$
$${\int_{t_0}^{t} y(t) \, dt} + RCy(t) = RCx(t)$$
If the time derivatives of the both sides of the last equation is taken, the resultant equation becomes;
$$y(t) + RC\frac{dy(t)}{dt} = RC\frac{dx(t)}{dt}$$
According to this differential equation which describes the input-output relationship of the given RC high-pass filter, the frequency domain analyses can be done.
The same process can be followed when analysing an RC low-pass filter too.


Frequencies of complex exponentials in Discrete Fourier Transform What I understand is that using DFT, we are representing a given discrete signal using a basis of complex exponentials of different harmonic frequencies. If I am taking 16 point DFT of a signal sampled at 1 KHz, k=0 corresponds to 0 Hz, k=1 corresponds to 62.5Hz etc. Is this correct?
In the figure, a cosine of 187.5Hz is shown and it is sampled at 1KHz. The sampled signal is actually the cosine part of the complex exponential of k=3 in a 16 point DFT. What confuses me is that, although the continuous time signal is periodic with frequency 187.5Hz, the sampled signal is not having the same period.
If the sampled signal is not of frequency 187.5Hz, how can we say that basis vectors of k=3 picks up the signal component of 187.5Hz?

It's not necessary for a discrete-time signal to have the same period as the continuous-time signal that it represents. As long as the continuous-time signal is band-limited and as long as the sampling frequency is greater than twice the bandwidth of the continuous-time signal, the discrete-time signal can perfectly represent the continuous-time signal.
In your example, the continuous-time signal is
$$x(t)=\cos(2\pi f_0t)\tag{1}$$
and the discrete-time signal is
$$x_d[n]=x(nT)=\cos(2\pi f_0Tn)\tag{2}$$
$x_d[n]$ perfectly represents $x(t)$ if $f_0T<\frac12$ is satisfied. In general, $x_d[n]$ is not even a periodic sequence. It is only periodic if $f_0T$ is rational.
With the choice $f_0T=3/16$ we make sure that $x_d[n]$ perfectly represents $x(t)$, and we also make sure that $x_d[n]$ is periodic with period $16$. Note that there can't be $3$ periods of the discrete-time signal in the DFT frame, simply because $16/3$ is not an integer.
Since $x_d[n]$ has period $16$, taking a length $16$ DFT results in a single frequency component. With $f_0T=3/16$, that component occurs at frequency index $k=\pm 3$, which exactly represents a continuous-time sinusoid at frequency $3/16\cdot f_s$ (with $f_s=1/T)$.


Fourier Transforms, symmetry, real/imaginary I was hoping to clarify if the following was correct:

A real function (neither even nor odd) in time exhibits conjugate symmetry in frequency, so the real part of the frequency response is even, and the imaginary portion is odd
A real, even function has a frequency response that is strictly real and even
An imaginary, odd function has a frequency response that is strictly imaginary and odd

-also: Can anything be said if the function were imaginary in time (neither odd nor even)? Would it be anti-conjugate symmetric, so the frequency response would be such that the real part is odd and the imaginary portion is even?

Your suspicion is correct. You can show this as follows:

$$X(\omega)=\int_{-\infty}^{\infty}x(t)e^{-j\omega t}dt\tag{1}$$

$$X^*(\omega)=\int_{-\infty}^{\infty}x^*(t)e^{j\omega t}dt\tag{2}$$

Now if $x(t)$ is purely imaginary we have $x^*(t)=-x(t)$, and, consequently,

$$X^*(\omega)=-\int_{-\infty}^{\infty}x(t)e^{j\omega t}dt=-X(-\omega)\tag{3}$$


Filter and crop an image I am a big fan of fft convolution. So, I mostly try to avoid non-fft convolution.

Suppose, I am trying to convolute an image with another image. Suppose the first image is $512\times 512$, and, the mask image is $520\times 255$. How would we convolve them in the frequency domain?

Two major issues are to be addressed here,

to apply fft, images must have dimensions which are the multiple of power of 2. Our image and mask may or may not have dimensions of power of 2
to apply the fft, the dimensions of the image and the mask must be equal.

So, in order to address both issues, we can do any one of the following,

we can take the max of the width and heights of both image and the mask, round them to next power of 2 and, pad them.
we can take the sum of the widths and heights of the image and mask, round them to next power of 2 and, pad them.

Questions:

So, in both of the cases we are actually obtaining oversized output images with padding. Am I right?
Suppose, I pad an image and get it FFT-filtered (Sharpened/Blurred/... and so on). So, I obtain a filtered image with padding. Now, I need to remove the padding to obtain the filtered image without padding. Am I correct?
If YES, should I just crop the image?
Or, is there any better way possible altogether from the ground up?

Suppose, this is the original image ($512 \times 512$),

And, suppose, the following is the filtered image with padding ($1024 \times 1024$),

Now, I want to obtain this,

How should I remove padding from this image?
Should I just crop, or, is there any better technique to get around this problem from start?

First, let's start with the wrong part of your post:

  to apply fft, images must have dimensions which are the multiple of power of 2.

This is incorrect. FFTs do not have to be a power of 2.  The only requirements for applying an FFT are that you can decompose the signal length into its prime factors and that you have an FFT that can deal with these prime factors.

  to apply the fft, the dimensions of the image and the mask must be equal.

This is correct. To do the point-wise multiplication in the frequency domain, it's easier to do when both image and kernel FFTs are the same size.

  we can take the max of the width and heights of both image and the mask, round them to next power of 2 and, pad them.

The power of 2, as pointed out before, is incorrect.

the width and heights of both image and the mask, round them is also incorrect.

Your second option:

  we can take the sum of the widths and heights of the image and mask, round them to next power of 2 and, pad them.

is closer to the truth, but still incorrect.

Suppose your image is $N\times M$ and your kernel is $K \times L$. To use the FFT to do linear convolution you need pad both of them to be $(N + K - 1) \times (M + L - 1)$. If you don't you'll get time aliasing from the circular convolution.

Now, to answer your questions:

  So, in both of the cases we are actually obtaining oversized output images with padding. Am I right?

Yes, to avoid circular convolution results (and the associated "time" / "spatial" aliasing), you need to pad.

  Suppose, I pad an image and get it FFT-filtered (Sharpened/Blurred/... and so on). So, I obtain a filtered image with padding. Now, I need to remove the padding to obtain the filtered image without padding. Am I correct?

Yes. If you're interested in an in-place replacement (so that the resulting image aligns with the original unfiltered image), then you need to remove the padding.

  If YES, should I just crop the image? Or, is there any better way possible altogether from the ground up?

Yes, just crop the image at the right place.


Stopping Viola Jones Face Detection early I want to try to improve an image processing technique for finding faces in images. I want to stop the Viola Jones face detector before it finds all faces, and do histogram equalization and other techniques on "candidate" faces. For example, say we have an image with one face. I would like to cut off facial detection before it has determined there is only one face, i.e. where there might be 3 prospective faces and Viola Jones has not yet taken out those regions which are not actually faces. Is this possible, and if so how do I do it?

Viola-Jones works in stages.  Only areas of the image that have a high enough "maybe" score get passed to the next stages.  If you implement your own Viola-Jones detector, then you should be able to stop between each stage and apply color/lighting corrections before going on to the next stage.

As an alternative to implementing your own, you might take a look at the openCV implementation.  You might be able to modify it for your needs.

The source is here for the cascade detector, which is the openCV Viola-Jones implementation.  You will need to examine the area after the comment //---------------------------------------- Classifier Cascade

That is where the actual detector is implemented.  I think the interesting places are operator and ocl_detectMultiScaleNoGrouping, and I really think operator is where you would need to start.