What is the cross-correlation of a WSS signal and the LTI filtered signal? I remember a result about the cross-correlation function of a stationary signal and that same signal, filtered through a linear time-invariant filter. There might have been an assumption that the signal was Gaussian. What is the cross-correlation function?

$$E\left\{ x(t) y(t) \right\} = \,\,\,?$$

where $y(t)$ is the output of an LTI filter with the input $x(t)$.

If the linear system (does not have to be time-invariant) is bounded-input bounded-output (BIBO) stable (meaning that if $|x(t)| \leq M$ for all $t$ for some finite positive
number $M$, then there is another finite positive number $M^\prime$ such that
$|y(t)| \leq M^\prime$ for all $t$), then when the input is a finite-variance
random process (a.k.a. second-order proecess), one can interchange the order of
the expectation operation and the integration/sum for the convolution.  That is, if
$$Y(t) = \int_{-\infty}^{\infty} h(s)X(t-s)\,\mathrm ds$$
where $h(t)$ is the impulse response of the LTI system, then
$$\begin{align*}
E[Y(t)] &= E\left[\int_{-\infty}^{\infty} h(s)X(t-s)\,\mathrm ds\right ]\\
&= \int_{-\infty}^{\infty} E\left[h(s)X(t-s)\,\mathrm ds\right ]\\
&= \int_{-\infty}^{\infty} h(s)E\left[X(t-s)\right ]\,\mathrm ds\\
\end{align*}
$$
and so if the input is a WSS process with $E[X(t)] = \mu_X$ for all $t$,
$$E[Y(t)] = \mu_X\int_{-\infty}^{\infty} h(s)\,\mathrm ds = \mu_Y$$
is also a constant.  In fact, $\{Y(t)\}$ is also a WSS process.  If we define the crosscorrelation function $R_{X,Y}(\tau)$ as
$R_{X,Y}(\tau) = E[X(t-\tau)Y(t)]$, then
$$\begin{align*}
R_{X,Y}(\tau) &= E\left[X(t-\tau)\int_{-\infty}^{\infty} h(s)X(t-s)\,\mathrm ds\right]\\
&= \int_{-\infty}^{\infty} h(s)E[X(t-\tau)X(t-s)]\,\mathrm ds\\
&= \int_{-\infty}^{\infty} h(s)R_X(\tau-s)\,\mathrm ds.
\end{align*}$$
In short,
$$R_{X,Y} = h*R_X.$$
For completeness,
$$R_Y = \tilde{h}*h*R_X = (\tilde{h}*h)*R_X = R_h*R_X$$
where $\tilde{h}(t) = h(-t)$ for all $t$ is the time-reversed impulse response
and $R_h = \tilde{h}*h$ is the autocorrelation function of the
deterministic signal $h(t)$.  Translated to the frequency domain, this gives the
power spectral density relationship
$$S_Y(f) = |H(f)|^2 S_X(f).$$


Problem understanding steps towards defining the convolution sum I have a problem understanding the reasoning behind a step that was taken to characterize a LTI system.

So, we were told the following:

For each integer $k$, we have the following function: $$   \delta_k[n] =
\begin{cases}
1,  & \text{for $n=k$} \\
0, & \text{otherwise}
\end{cases}$$

Thus, if we write $x_k = x[k]$, we have that: $$x[n] = \sum_{k\ = \ -\infty}^{+\infty} x[k]\delta_k[n]$$ where x[n] is a discrete signal.

Finally, we can exploit the linearity of the system as follows: $$H\{x[n]\}\  = \ H\left\{ \sum_{k\ = \ -\infty}^{+\infty}x_k \ \delta_k[n] \right\}$$
$$=\sum_{k\ = \ -\infty}^{+\infty}x_k \ H\left\{\delta_k[n]   \right\}$$ where $H$ is the system.

It is this last step that confuses me a lot. How does linearity explain that we have $\sum_{k\ = \ -\infty}^{+\infty}x_k \ H\left\{\delta_k[n]   \right\}$ ? It only makes sense if we consider $x_k$ as a scalar (or a constant if you prefer), but we defined $x_k$ as a signal, so it can't be a constant. So, why is this correct?

PS: Just for clarification, we have that $\delta_k[n] = \delta[n-k]$ where $\delta[n]$ is the Kronecker-delta function.

The point of the whole exercise it to prove that you can calculate the output of an LTI system by convolving the input with the impulse response of the systems. This is done in the following steps.

Rewrite the input as a sum of delay delta impulses
Put that sum through the systems
Exploit LTI, swap the order of applying the system and the sum. Response of sum = sum of responses !

Let's go through the last step a little slower. We have
$$H\{x[n]\} = H\left \{ \sum_{k=-\infty}^{\infty} x[k]\cdot \delta [n-k]\right \}$$

First let's sawp summation and system, assuming that the response to the sum is the same as the sim of the responses (per LTI definition)

$$H\{x[n]\} =  \sum_{k=-\infty}^{\infty} H\left \{ x[k]\cdot\delta [n-k]\right \}$$

Now comes the tricky part: you have two time-like variables in the thing you need to transform: n and k. k is just the counting index for the different basis functions. n is actually the transform variable. You can think of this as having an infinite number of functions and you need to calculate the response for each one of them (and then sum the results). However, the response of each individual function is independent from the other ones. Hence, you can pull anything that's only dependent on k out from the transform. Different way of of looking at it: if H{} where an integral your integration variable would be $\partial n$ not $\partial k$. So we have

$$H\{x[n]\} =  \sum_{k=-\infty}^{\infty}  x[k]\cdot H\left \{\delta [n-k]\right \} = \sum_{k=-\infty}^{\infty}  x[k]\cdot h[n-k]
$$

where $h[n]$ is the impulse response of the system


Solving simple transfer function I was hoping that you can help me understand how to find poles and zeros in this transfer function.

$$
H(s)=\frac{1}{(s-s_{\alpha1})(s-s_{\alpha2})}
$$

I know that this function does not have any zeros, because of the 1.

From the solutions I can see that the poles are $s_{\alpha1}=-1+j$ and $s_{\alpha2}=-1-j$ and to tell you the truth, I have no idea how to get solution for the poles.

I tried to solve it using quadratic equation, but the only thing I get is a big mess.

Thanks!!

Because $s-s_{\alpha1} = 0 $ and $s-s_{\alpha2} = 0 $ so poles: $s_1,s_2 = s_{\alpha1},s_{\alpha2} $
What numerical values you have for $s_{\alpha1},s_{\alpha2} $ are the What numerical values you are looking for.


Where can I find information about how to implement IIR filters on micro-controllers? I am trying to implement an IIR filter on a LPC1769 (I don't know if it's relevant). I already have the H(z) so all the design part is prety much solved. What I am looking for is information about:

How to select an appropiate word lenght for the coefficients, acummulator, etc.
How to select the most appropiate structure (parallel,cascade, direct forms, indirect forms....) (or some guidance at least)
Some understanding on how to avoid overflow (scalling the coefficients or signal)
A deeper explanation of the fdatool parameters than the one given in real-time digital signal processing.
etc....

I've been using the fdatool of matlab to play around with the quatization of the filter. But I don't know the implications of these or what they actually mean.

I've also been reading a few chapters of this book: real-time digital signal processing : implementations and applications sen m. kuo, but I am looking for more detail in these implementation aspects. So could you tell me of some book, application note, paper or whatever that explains these issues?

I am triying to implement an A-weighting filter (the transfer function is in the link). The application is mainly taking data from the adc, perform the A-weighting and taking the RMS value of this signal.

Here are 3 points that you must consider :

1 - What is the minimal latency you can tolerate? This is important, as it's important to know whether to need process every sample as fast as possible or if you can perform batch processing which is typically more efficient but creates latency .

2 - No matter if you implement the IIR in fixed-point or floating-point, it is always a good idea to split the filter in sub-filters. This minimizes quantization issues that can occur even if you use double-precision arithmetic. One effective strategy is to cascade the filter in second-order sub-filters. Check out https://www.mathworks.com/help/signal/ref/sosfilt.html

3 - If you implement in fixed-point, it is a good idea to create a fixed-point model in Matlab/Python/Scilab/tool of your choice and compare it to the "ideal" floating-point filter before implementing it in your microcontroller.


Practical diversity combining at receivers (Maximal Ratio Combining & Equal Gain Combining) I am trying to perform practical diversity combining implemented in GNU Radio by using Ettus N210 USRPs, one as a transmitter and two as receivers.

To compare Equal Gain Combining (EGC) and Maximal Ratio Combining (MRC) techniques, I simply add the same received signal from two different receivers constructively for EGC and multiply each branch with a gain derived from square root of its signal-to-noise ratio (signal+noise/noise) before adding for MRC.

Contrary to the literature, for a few cases MRC is not performing better than EGC. I have noticed that when there is no signal, the noise level for both receivers is not same. In theory, it is assumed that all branches have equal mean power and noise level.

Now, I am confused. Should I attenuate/amplify the received signal (signal+noise) from each receiver to make the mean signal power same across both receivers or should I attenuate/amplify the received signal (signal+noise) in a way that noise level become same in the signal from both receivers before performing EGC.

Once I am done with EGC, should I simply multiply the weight at each branch to perform MRC?

Since, in theory, it is assumed that long term average SNR and noise remain same across all branches, first, I configured my USRPs to have same noise levels and then placed them in a way that average SNR across both branches also remain the same.
Finally, simply combining them or first multiplying instantaneous SNR factor with each branch leads to EGC and MRC, respectively.


Rayleigh fading channel as a scalar quantity In most cases rayleigh channel is a vector. can a rayleigh channel be a scalar quantity i.e., can we have only one element in rayleigh fading channel.

Rayleigh channel is defined between transmit and receive antennas, if there is one transmit and one receive antennas then there is one Rayleigh fading coefficient.

If there is one transmit and multiple receive antennas then there is a vector or fading channel.

If there is one receive  and multiple antennas then there is a vector or fading channel.

If there are multiple trasmit and receive antennas then you have a Rayleigh fading matrix represting the wireless channel


Audio signal comparison for automatic singing evaluation I am trying to compare audio waveform from input waveform (A) to preset waveform (B)from in the device.

I try to lookup everywhere. Most likely someone give such a source code as follow :

    function compireAudio(){
    var audio1 =    "http://soundjax.com/reddo/86502%5Ealarm.mp3";
    var audio2 =    "http://soundjax.com/reddo/44368%5EALARME.mp3";
    var i,j,d;
    var matching = 0;
    var t = 0;
    var i,j,d;
    var matching = 0;
    var t = 0;
    var audio1Arr = Array();
    var audio1Len = audio1.length;
    for (i = 1; i<=audio1Len; i++) {
        //reverse so its like a stack
        d = audio1.charCodeAt(audio1Len-i);
        for (j = 0; j < 8; j++)
        {
            audio1Arr.push(d%2);
            d = Math.floor(d/2);
        }
     }
     var audio2Len = audio2.length;
     for (i = 1; i<=audio2Len; i++) {
       //reverse so its like a stack
       d = audio2.charCodeAt(audio2Len-i);
       for (j = 0; j < 8; j++)
       {
         if(d%2 == audio1Arr[t])
         {
             matching++;
         }
         d = Math.floor(d/2);
         t++;
       }
     }
   var avarage = Number(matching)/((Number(t)+Number(audio1Arr.length))/Number(2))*Number(100);
    alert('The Matching with the two audio is '+avarage+' %.');

And most people said this is not correct because it compare the stream instead of the audio content.

Anyway, I have think of compare every samples after normalization (ensure that amplitude is not a factor here). However, a phase is another problem here.

I have solution is cross-correlation between A and B and find the biggest convolution as output. In this case, they should be in same phase.

In short.
WAS Question is What are the common methods for two waveform comparison ?

UPDATE Maybe I should ask What are the common characteristic of a audio waveform for comparing audio waveform ?

P.S. I read this Question and answer. The bottom method seems work, but the comment said No.

Please update your question to mention which kind of comparison you want to perform. You seem to want your comparison to be robust to amplitude changes and time delays, but what about:

Audio coding. Should a signal and a MP3 version of itself be considered as similar or different?
EQ and filtering. Should a signal and the same signal with the high frequencies shifted by 10° be considered similar or different?
Dynamics change. Does applying a time-varying gain factor to one signal make it different from the original?
Time and pitch shifting. Should a signal and a slightly accelerated version of itself be considered as similar or different?
Noise and recording conditions. Should two versions of the same signal recorded on different equipment, microphone type, etc. be considered as similar or different?
Timbre. Should two versions of the same melody played by two different instruments be considered as similar or different?

Depending on your application (assessing degradations through a signal processing chain, identifying copyrighted work in an audio stream, shazam-like song identification, plagiarism detection on music...), you might want your comparison to be robust to some of these transformations.

Your approach consisting in normalizing the signals is not very robust. If you have two identical signals $x$ and $y(t)$, and if you change $x(0) = 2 \max_t y(t)$, the normalized version of $x$ is twice smaller than the normalized version $y$ and the comparison will fail.

Using cross-correlation to look for a time delay is a sound approach.


Aliasing due to the Convolution of Gaussian Functions For two discrete-time sequences $f[n]$ and $g[n]$, their linear convolution $(f*g)[n]$ is also given by
$$f*g = \mathcal{F}^{-1}(\mathcal{F}(f) \cdot \mathcal{F}(g)),$$
where $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the DTFT (discrete time Fourier transform) and IDTFT (inverse discrete-time Fourier transform) respectively.

Let us consider two Gaussian sequences $f[k]$ and $g[k]$ in the momentum space and we want to convolve them using the formula above.

import numpy as np
import matplotlib.pyplot as plt

L = 17.5
N = 48
dx = L/N
dk = 2*np.pi/L
x = np.arange(N)*dx - L/2
k = 2*np.pi * np.fft.fftfreq(N, dx)

alpha = 1
beta = 1

def f(k):
   return np.sqrt(np.pi/alpha) * np.exp(-k**2/(4*alpha))

def g(k):
   return np.sqrt(np.pi/beta) * np.exp(-k**2/(4*beta))

fs = np.array([f(k_) for k_ in k])
gs = np.array([g(k_) for k_ in k])

fftconvfg = np.fft.ifft(fs * gs)

plt.plot(ks, np.fft.fftshift(fftconvfg))

I get the following plot.

Question

How can I demonstrate that my plotting is (or, is not) subject to the aliasing artifacts? And if there are aliasing artifacts, how can I get rid of it?

Short answer: For the convolution of two Gaussian sequences there will be aliasing.

If the sequence produced as the result of linear convolution has infinite domain of support, then there will always be aliasing in the Fourier domain implementation of the linear convolution using the circular convolution implied by the DFT (discrete Fouerier transform) which is used to represent the theoretical DTFT that provides the theorem in your question.

Therefore, if you want an alias free implementation of the theorem in your question, then the resulting signal (result of the true linear convolution) must have a finite length.

That being said, for many practical cases, the amount of distortion due to aliasing will be small, so small to be lost behind any ADC step size or even smaller than represantable by the floating point number system that's used to implement the algorithm. In such cases it can practically be considered to be a sufficient realization of the true lienar convoluion.


Regarding phase and frequency components in the signal As I am not a signal processing student, I have limited understanding of the concepts.

The general formula for sine wave is

$$s(t)= A \sin(\omega t-p)$$

where $A$ is amplitude $\omega$ is frequency and $p$ is phase.

Is there any possibility that the phase of the signal remains constant even if the frequency of the signal is changed?
If we observe a real-time signal the signal can be a periodic sine wave signal but not smooth curve as it is accompanied by noise. In such cases how the phase of the signal is varied? Can the noise effect the phase of signal?

I find it easier to understand the phase relationships by observing complex tones rather than sines and cosines, which are related using Euler's identity:

$e^{j\theta}= \cos(\theta)+j \sin(\theta)$

$e^{j\theta}$ is simply a phasor with magnitude 1 and angle $\theta$. If that angle is changing with time at a constant rate ($e^{j\omega t})$ the real portion will be a cosine and the imaginary portion will be a sine. If that phasor that is spinning at a constant rate does not start at magnitude = 1 and angle = 0 at time $t = 0$, then we have a phase offset ($e^{j\omega t + \phi})$ where $\phi$ here represents the starting phase at $t = 0$.

Thus a single impulse in the frequency domain is a single spinning phasor in the time domain as shown in the figure below. From that figure the I axis represents the real part of the spinning phasor, which would be a cosine wave, or a sine wave if the starting phase at $t=0$ was $\pi/2$ since $\cos(\theta+\pi/2)= \sin(\theta)$

So to answer your questions with reference to the explanation above:

  Is there any possibility that the phase of the signal remains constant even if the frequency of the signal is changed?

By definition of frequency being phase changing with time, the only signal that has a phase that remains constant is a DC signal. So if you stop the phasor from spinning (at any angle) this would be frequency = 0 (DC).  You can have a constant phase offset for all frequencies, which is the phase you are referring to- this simply means at t=0 every frequency in your scenario has the same phase; but then as time proceeds each frequency would have a phase that would change at a different rate (the slope of the phase vs time IS the frequency!).

  If we observe a real-time signal the signal can be a periodic sine wave signal but not smooth curve as it is accompanied by noise. In such cases how the phase of the signal is varied? Can the noise effect the phase of signal?

Yes of course it does effect the phase, and this is what is called "phase noise". Refer to the plot above of the spinning phasor in the time domain, and now imagine noise being added to it. The noise would appear as a random vector on the end of the phasor (proportional to how strong the noise is), and that vector would change randomly with time. Often this is shown as a round "noise cloud" that is on the end of the vector signifying its statistical distribution with time (but of course at any given instant the noise vector would be at one particular amplitude and phase). This noise cloud (If a complex Gaussian noise) would equally effect the amplitude and phase of the signal, causing both AM (amplitude modulation) and PM (phase modulation), what we would refer to as AM noise and Phase Noise.


What is the amplitude range of a spectrogram? Is there an upper and lower bound to a spectrogram representation of an audio signal? The range of the amplitude values of the audio signal is from 0 to 1.

Recall how a spectrogram is computed. You place a window of length $N$ around a sample of interest in your signal $x$ and then compute its DFT. This corresponds to one column of the spectrogram. In symbols, a fixed column is given by a DFT operation:

$$
X(k) = \sum_{n=0}^{N-1} x(n) w(n) e^{\frac{-j 2 \pi k n}{N}}.
$$

I believe you want an upper and lower bound on the magnitude of the spectrogram. An obvious lower bound is $0$. To find an upper bound:

\begin{eqnarray}
|X(k)| &=& \left | \sum_{n=0}^{N-1} x(n) w(n) e^{\frac{-j 2 \pi k n}{N}} \right|\\
&\leq& \sum_{n=0}^{N-1} |x(n)||w(n)|\cdot 1\\
&\leq& \sum_{n=0}^{N-1} 1\\
&=& N.
\end{eqnarray}
where I assumed that the window function has a maximum value of $1$.

This means that the magnitude of any element in the spectrogram matrix must be between $0$ and $N$ where $N$ is the window length.

Note that there are other ways of scaling a spectrogram which will change these calculations. Some people plot squared magnitude instead of magnitude, or $20\;\log_{10}(\cdot)$ of the magnitude. Some people like putting a $\frac{1}{N}$ normalization. Matlab provides 'psd' and 'power' options and scipy provides 'density' and 'spectrum' options which scale things differently.


Generate fourier transform signal I'm newbie in DSP (maybe the question title is misleading for this reason, I apologize for this, please feel free to edit it) but not in programming. I want to experiment with Fourier (or FFT) transformation. Basically I need the generator which will give me signal looking like this:

I would also like to add/remove to/from this signal some noise, This is example for sine wave but I hope you get the point.

Is this possible for someone whom only knowledge about Fourier transformation is that is composition of sines and cosines with the aim of being able to create signal of any shape? I've been searching for some code snippet where I can change parameters and see how the signal is changing but my google-fu was not enough for this. Preferred language is python (where i can import module and call function) or R or something multi platform (open source) and simple and fast. Thank you.

You can create a sine wave in numpy.  Think of it as if you were sampling data from a signal, only you are sampling a mathematical function:

sample_rate = 1024
dt = 1.0/sample_rate
t = np.arange(sample_rate)*dt  # 1 second of samples
freq = 5
amp = 1.0
sine1 = amp*np.sin(2*np.pi*freq*t)
sine2 = .5*np.sin(2*np.pi*15*t)
sinsum = sine1 + sine2


What kind of effect is this and what is it caused by? In the screenshot below I have two audio signals. The upper one is the audio signal that is created  with python and then played, whilst the lower one is the "same" signal but recorded via microphone on an android device. Tried my best at aligning them in audacity but I am unsure where the "fade-in" comes from. The same thing happens on the other end of the signal, but in reverse of course.

Also, there is no other noise in the recording, as I am recording from the audio output relayed to the input and not directly from the speakers (so no background noise).
Don't know which terms to search for. Any help would be much appreciated.

Looks like pre-ringing. That happens whenever you apply a linear phase filter to a signal.  In your case that could be the sigma-delta converters or any type of sample rate conversion or other unknown processing that the operating system or the drivers doo.


Can we choose a sampling frequency to remove unwanted noise at a specific frequency? I am studying for my exam in signal processing. In one of the old exam papers I am told to find a sampling frequency, which will remove 80 Hz noise. The filter the exam question is based around has an impulse response $h[n]$ of
\begin{bmatrix}
-3&0&-3&0
\end{bmatrix}
My approach is to just choose the sampling frequency at a multiple of the noise frequency. In this case, I would choose it to be 320 Hz, since 4 samples gives a sum of 0. When doing a DFT with the known impulse response it would also result in a net 0 frequency response.
Does this seem like a sustainable approach?

Or another way, visually, with the plot of the magnitude response of your filter shown below you see that it goes to zero at $\pi/2$.

So what you could do is use the relationship between the angular frequency $\omega$ in [radians/sample] and the relative frequency $f$ in [cycles/sample]
$$
\omega = 2\pi f = 2\pi \left(\frac FF_s\right)
$$
Then compute $F_s$ for $F = 80 \ \rm Hz$ and $\omega = \pi/2$. You then get
\begin{align}
\frac \pi 2 &= 2\pi \left(\frac{80}{F_s}\right)\\
\implies \frac 12 &= 2\cdot \left(\frac{80}{F_s}\right)  \\
\iff 80 &= \frac 12\cdot \left(\frac{F_s}{2}\right)
\end{align}
Giving you the $320\ \rm Hz$.
In MATLAB
You can plot the magnitude response from the frequency response using MATLAB's freqz function where you can specify the frequency grid from $0$ to Nyquist (i.e. $\pi$). With your resulting $320\ \rm Hz$, you can run the function with this $F_s$ and see that $80\ \rm Hz$ is indeed suppressed.


How can I do QPSK in MATLAB Could you please explain the step below, what is the benefit of it and how does it work. .
I do qpsk demodulation, do I started doing the meshgrid,

%%%%%%%de-qpsk
 sM = 2;
 [x,y] = meshgrid((-sM+1):2:(sM-1),(-sM+1):2:(sM-1));
 alphabet = x(:) + 1i*y(:);
 datadeqpsk=[];

 for k=1:N
    datadeqpsk(k) = alphabet((dataqpdk_be(k) alphabet)==min(dataqpdk_be(k)-alphabet)); % dataqpdk_be is the data
 end

datadeqpsk_R=real(datadeqpsk);
datadeqpsk_I=imag(datadeqpsk);
dataend=[];

Till now, all is ok, so what the benefit of below step? how does it work? could you please explain?

     for ij=1:N
         dataend((ij-1)*2+1)=(datadeqpsk_R(ij)+1)/2;
         dataend(ij*2)=(datadeqpsk_I(ij)+1)/2;
    end
      R_D_data_out=[R_D_data_out dataend];

thank you very much

Your code is performing the QPSK demodulation. it's assigning the odd value for real value and even values for Imaginary valued Or, the opposite, try to re-check. And Then make the QPSK demodulation.

About you question, how to do the qpsk demodulation,  you can do it as below:

% suppose you have data for 48, 2-bit symbols
       data = randi([0 1],96,1);
 % Create a QPSK modulator System object with bits as inputs and Gray-coded signal constellation
       qpsk_mod = comm.QPSKModulator('BitInput',true);
 % Modulate and plot the data
       modData = step(qpsk_mod , data);
       scatterplot(modData)

Go to MATLAB and check the help of comm.QPSKModulator, then you can explore the demodulation easily.

you can also check this web page for deep understanding https://www.gaussianwaves.com/2010/10/qpsk-modulation-and-demodulation-2/

Good Luck


Brief introduction to wireless communication and communication terminologies In the Question Design of equalizer for wireless communication the diagram of a basic communication system is presented with informative points and discussions. My Questions are:

Considering the channel modeled as a FIR filter having a single input and single output. If the impulse response (channel coefficients) take values in the complex domain then am I looking at a Rayleigh Fading channel?
What constitutes a multi-path channel? Would multi-path channel be typically not a single input and single output?
In the diagram, the input to the channel is Tx. Is Tx digital? I mean can it take bits or multiple level bits? Suppose, we want to send an image or a real-valued time series through a communication system, then is it transmitted in real values or digital in wireless communication?
At which stage is digital modulation performed? Then the output of the modulation should be digital. So, is this digital output the Tx that is sent to a channel?
What is the difference between modulation and encoding? When we are doing encoding, we are basically transforming the information from one representation to another. Isn't this what modulation also does?

Well, these are very fundamental questions, and I think they are written in any good book. But maybe it's hard to extract the exact information from them, which you ask for. So, let me briefly answer:

In wireless, the channel is usually modelled as a FIR filter with impulse response $h(\tau)$. But, this channel might change over time (due to e.g. mobility), to make it a time-variant impulse response $h(t,\tau)$. Moreover, the channel is often modeled as a tapped delay line, $h(t,\tau)=\sum_ia_i(t)\delta(\tau-\tau_i(t))$ where $a_i(t)$ are time-varying gains of the taps and $\tau_i(t)$ are the time-varying delays of the taps. The fact that $a_i(t)$ has complex values does not make it Rayleigh fading. Rayleigh fading describes the temporal distribution of $a_i(t)$. If over a long time range the distribution of each path gain is a Gaussian, it becomes Rayleigh fading. The doppler spread of the Rayleigh fading describes how quickly the value $a_i(t)$ changes over time. In the special case of block-fading $a_i(t)$ is a constant (the channel remains static) for the whole transmission. Then, when a new transmission starts, a new channel is generated. Then, the realizations $a_i$ follow a Gaussian, to make it block Rayleigh fading.

A multipath channel corresponds to the tapped-delay line model. Have a look at one of my articles about multipath propagation, maybe it helps. The number of inputs and outputs depends on the number of TX and RX antennas of your system (talking about MIMO systems here) or the number of users. For single-antenna systems there is just one input and one output.

The input to the physical channel is analog, if you model it as a tapped delay lines. However, you can assume transmission in passband, then all signals are real. Or, you can do an equivalent analysis in baseband, where all signals are complex but the math is easier to treat. For the description of transition from baseband to passband, you can look at another article of mine.
Though, there are other channel models such as BSC or BEC, though these are not used to model the transmission of an analog signal over the wireless physical channel, but abstract away from that.

Modulation is performed after encoding ;-). Modulation means to generate a signal which is adapted to the channel, a a function of the transmit bits. It can e.g. include mapping the bits to complex QAM constellation points.

Encoding takes your bits you want to transmit (the payload) and adds redundancy via a channel code. This way, if at the receiver side, one bit is lost, it can be recovered by appropriate decoding.


Calculating the inverse filter for the (exponential) sine sweep Method I am currently working on my Bachelor-Thesis about Real-Time convolution and Impulse Response Measurements. After reading different papers about the (exponential)sine sweep method I didn't find an answer on how to calculate the inverse filter for deconvolving the impulse response.

In the papers I read, it is described as the Time Reversal Mirror and after that some kind of scaling.

Could someone explain, how to calculate the inverse filter for a given sine sweep?
If you could add an example, or an algorithm it would be appreciated

Assuming that your Exponential Sweep Sine was generated using the formula:

$$x(t)=\sin\left(\frac{2\pi f_1 T}{R}\left(e^{\frac{t R}{T}} -1\right) \right)$$

where:

$f_1, f_2$ - Initial and final frequency of the sweep

$T$ - Duration of the sweep

$R = \ln\left(\frac{f_2}{f_1} \right)$ - Sweep rate

Then the inverse filter is calculated by scaling the amplitude of time reversed $x(t)$ by:

$k = e^{\frac{tR}{T}}$

Which will result in an exponentially decaying sweep:

$f(t) = \frac{x_{inv}(t)}{k}$

Example in Python:

#!/usr/bin/env python

from __future__ import division
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

def dbfft(x, fs, win=None):
    N = len(x)  # Length of input sequence

    if win is None:
        win = np.ones(x.shape)
    if len(x) != len(win):
            raise ValueError('Signal and window must be of the same length')
    x = x * win

    # Calculate real FFT and frequency vector
    sp = np.fft.rfft(x)
    freq = np.arange((N / 2) + 1) / (float(N) / fs)

    # Scale the magnitude of FFT by window and factor of 2,
    # because we are using half of FFT spectrum.
    s_mag = np.abs(sp) * 2 / np.sum(win)

    # Convert to dBFS
    ref = s_mag.max()
    s_dbfs = 20 * np.log10(s_mag/ref)

    return freq, s_dbfs

if __name__ == "__main__":
    # Sweep Parameters
    f1 = 10
    f2 = 100
    T = 3
    fs = 1000
    t = np.arange(0, T*fs)/fs
    R = np.log(f2/f1)

    # ESS generation
    x = np.sin((2*np.pi*f1*T/R)*(np.exp(t*R/T)-1))
    # Inverse filter
    k = np.exp(t*R/T)
    f = x[::-1]/k
    # Impulse response
    ir = sig.fftconvolve(x, f, mode='same')

    # Get spectra of all signals
    freq, Xdb = dbfft(x, fs)
    freq, Fdb = dbfft(f, fs)
    freq, IRdb = dbfft(ir, fs)

    plt.figure()
    plt.subplot(3,1,1)
    plt.grid()
    plt.plot(t, x)
    plt.title('ESS')
    plt.subplot(3,1,2)
    plt.grid()
    plt.plot(t, f)
    plt.title('Inverse filter')
    plt.subplot(3,1,3)
    plt.grid()
    plt.plot(t, ir)
    plt.title('Impulse response')

    plt.figure()
    plt.grid()
    plt.semilogx(freq, Xdb, label='ESS')
    plt.semilogx(freq, Fdb, label='Inverse filter')
    plt.semilogx(freq, IRdb, label='IR')
    plt.title('Spectrum')
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dBFS]')
    plt.legend()

    plt.show()

And output:

For Robert, here is the spectrum plot for case without amplitude modulation of the inverse filter:

Relevant literature:

  Q. Meng - Impulse Response Measurement with Sine Sweeps and Amplitude
  Modulation
  Schemes

  A. Novak - Nonlinear System Identification Using Exponential
  Swept-Sine
  Signal

  K. Vetter - ExpoChirpToolbox - a Pure Data implementation of ESS
  impulse response
  measurement


Is there a way to calculate modulation depth after the modulation process? Suppose I have some amplitude modulated carrier wave $x$. If I have no access to the original message wave nor the unmodulated carrier, is it some how possible to calculate what is the modulation depth of $x$?

If I have understood this correctly, amplitude modulation depth is given by
$M/A$, where $M$ is peak amplitude of the message wave and $A$ is the peak amplitude of the carrier wave. From $x$ I will know that $M/A$ equals to some $\lambda$ (when I check its peak amplitude), but M and A could be anything.

Please correct me if I have understood something incorrectly.

If you look at the envelope of an AM signal then for 100% modulation the amplitude of the envelope will vary between 0 and 2 * Ac, where Ac is the amplitude of the carrier. It's a linear relationship, so for e.g. 50% modulation the amplitude will vary between 0.5 * Ac and 1.5 * Ac. So as long as you know that you have a full scale modulation signal then you can derive Ac from the max and min envelope amplitude (mean of Amax and Amin) and from there you can calculate the modulation depth (Amax / Ac).


Expression for discrete fourier transform of linear ramp I am trying to compute a single coefficient of the DFT of a linearly ramping sequence, $x[n]=an$ where $a$ is a constant that changes from sequence to seqeunce. I have looked at loads of DFT transform property and pairs tables, but this is a simple case I have not been able to find an analytical expression for. I don't want to have to compute the entire DFT of the sequence just to extract a single coefficient.

My first approach was to notice that multiplication by $n$ in the time domain is like differentiation in the $k$-domain. However, since the DFT of $a$ is an impulse function at $k=0$, I haven't found a way to come up with an analytical expression to take the derivative.

In my application, the number of samples in the sequence, $N$, and the coefficient of interest, $k$, will frequently change. It's just such a simple case that, given $a$, there must be a way for me to predict the value of the DFT of the sequence at a particular value of $k$ without actually doing the transform. Is the only way to do this to numerically create a look-up table of all possible $N$ and $k$ of interest, and then scale by $a$?

There is indeed a better way; you can derive the analytical expression for the DFT of a ramp. Let's start with the discrete-time Fourier transform (DTFT) of a finite length sequence (i.e., the sequence is zero outside the interval $[0,N-1]$):

$$X(e^{j\omega})=\sum_{n=0}^{N-1}x[n]e^{-jn\omega}\tag{1}$$

We also need the DTFT correspondence

$$nx[n]\Longleftrightarrow j\frac{dX(e^{j\omega})}{d\omega}\tag{2}$$

The DTFT of the constant signal $x[n]=a$, $n\in [0,N-1]$ (and zero otherwise) is

$$X(e^{j\omega})=a\sum_{n=0}^{N-1}e^{-jn\omega}=\begin{cases}aN,&\quad\omega=0\\\displaystyle a\frac{1-e^{-jN\omega}}{1-e^{-j\omega}},&\quad\text{otherwise}\end{cases}\tag{3}$$

From $(2)$ and $(3)$ we get

$$\begin{align}a\sum_{n=0}^{N-1}ne^{-jn\omega}&=j\frac{dX(e^{j\omega})}{d\omega}\\&=-a\frac{Ne^{-jN\omega} (1-e^{-j\omega})-(1-e^{-jN\omega})e^{-j\omega}}{(1-e^{-j\omega})^2},\quad\omega\neq 0\end{align}\tag{4}$$

which is the DTFT of the given sequence. For $\omega=0$ we simply have

$$a\sum_{n=0}^{N-1}n=a\frac{N(N-1)}{2}\tag{5}$$

Now we only need to remember that for finite length sequences, the DFT is a sampled version of the DTFT with $\omega_k=2\pi k/N$. From $(4)$ and $(5)$ we obtain

$$\text{DFT}_N\{an\}[k]=\begin{cases}\displaystyle a\frac{N(N-1)}{2},&\quad k=0\\-a\displaystyle\frac{N}{1-e^{-j2\pi k/N}},&\quad k\in[1,N-1]\end{cases}\tag{6}$$


Zeros inserted in QAM Constellation I need to ask about insertion of Zero Symbols after QAM modulation in a OFDM system. Such that the 4 QAM has a constellation with 4 corner and 1 center symbol indicating zero. What advantages or disadvantages woild that offer ?

Is this done for guard bits or as pilot or what ?

So I need to add more information to the question. I just had a tutorial on OFDM system on MATLAB and what my Tutor told was that in addition to the Cyclic Prefix, we may also include these zero symbols after we have modulated the bits with QAM. Which appear as Zero Points in Constellation. He was saying that they have the purpose of controlling the ISI as well. But I was confused because ISI was controlled by Cyclic Prefix. ?

Also I am not sure about what effects this would have on OFDM Spectrum.

As others have mentioned, it sounds like you are talking about the null carriers in an OFDM system:

These bins are used to mitigate leakage from adjacent RF bands 2. Perhaps your tutor meant ISI in the sense that if adjacent RF bands were to spill into your OFDM signal, you would have some interference with your symbols on the edge of your bandwidth, but this is not ISI this is just, well, I (interference). The null carrier at zero is inserted because this carrier will experience significant interference in direct digital down conversion architectures (i.e. cannot recover modulated data when brought to DC at baseband).

In your case all the red/blue are QAM modulated data (if you use QAM for your pilot) the rest are left zero. The frequency spectrum would be exactly what you see in the figure above. It would be the same overall bandwidth, but with a fraction of the bandwidth used. It can be confusing because we use the term "guard band" when we mean frequency guard band which is different than what we mean by the cyclic prefix time domain "guard band" which mitigates the effects of ISI.


The desired frequency-response specifications The complex function $ D (e^{-jw})$ is defined on the domain of approximation $\Omega$ .In most cases the domain $\Omega$  is the union of several disjoint frequency bands which are separated by transition bands where no
desired response is specified .We denote the union of all passbands by $\Omega^p$ and stopbands by $\Omega^s$. If the designed filter is to have real-valued coeffcients only the domain $\Omega\cap[0 ,\pi]$ is
considered,what is the type of coefficients which used in the domain $[0 ,2\pi]
$ or$ [-\pi,\pi]$?

Thanks in advance.

If you have to consider the whole frequency spectrum range in $[0,2\pi]$ for the design of the discrete-time filter, without assuming any type of symmetry, then you are considering the most general case of the filter and its coefficients will be complex valued. And nothing more can be said about them, unless you impose further constraints on the impulse and frequency responses of the filter.


why is the spectrogram from matlab different from audacity Something related

Here are the spectrogram plots:

Matlab version:

Audacity version:

the darker parts towards the right of the figure show speech. Silence was added to the left of the figure in audacity. This silent portion does not show any thing in the spectrogram as plotted by audacity which is as expected. but the same silence portion is covered with blue dots in case of the matlab plot. Why didn't Matlab plot it with all white (since it is silence added through software and is not supposed to have any energy)?

---EDIT: response to pichenettes' answer ---

I think you are right about the first explanation

The previous spectrogram was obtained by changing the time domain wave form view into spectrogram view immediately after I added silence. Actually in my case I had added silence at the end of the file to make it's length exactly equal to 2500ms, then cut and pasted to the beginning. So in my case it actually did show the time domain waveform of silence as well after insertion. But it seems that it does not actually include the 0's of silence in the wavfile until I export it. When I see the spectrogram of the wavfile after an export and a re-import, I an see blue dots in the audacity version as well

And teh difference in colour is probably due to the difference in colour palettes used, as pinchette had pointed out in his other answer

I see two possible explanations.

Explanation 1:

Audacity is a non-destructive/non-linear audio editor. It primarily manipulate edit decision lists of audio signals, rarely the signals themselves. So when you nudge a block of audio 1.4s to the right, it does not fill the gap with 0 (silence). Instead, it just stores somewhere "this block plays 1.4s from the beginning of the timeline". For example, I have created here a 30s sine wave tone and nudged it at time offset 27s. Notice that the gap between the origin and the waveform is empty rather than showing a flat line indicating a 0 amplitude.

Now, if you apply the "Mix and render" command, all editing and mixing commands will be actually performed on the signal.

Explanation 2:

Digital silence is a sequence of 0. This causes a problem during the computation of the spectrogram because this involves the computation of log (0) to get to the dB scale. Audacity might use a special color in its palette for digital silence; while matlab might add a very low noise.


how does this equation correspond to smoothing? Please help me understand smoothing of data. This is a follow up to my previous question posted here. Especially the top answer by Junuxx where he says a way of smoothing a function $f(x)$ is:

$$
f'[t] = 0.1 f[t-1] + 0.8 f[t] + 0.1 f[t+1]
$$

here we can see that for every point in $f[x]$, we are taking a weighted average of that point and its two adjacent points, to get a smoothed version of $f[t]$ called $f'[t]$.

A paper on speech enhancement explains that an equation of the form

$$
y[i] = a[i]y[i-1] + (1 - [i]) x[i]
$$

helps us get the value of y as a recursive smoothing of x. Here $a[i]$ acts as a smoothing parameter and it is itself calculated as

$$
a[i] = \alpha + (1 - \alpha)p[i]
$$

where $p[i]$ is calculated elsewhere and alpha is a constant. $y[i]$, $a[i]$, and $x[i]$ are all arrays with $i$ elements.

How can I relate this equation of $y[i]$ with the equation of $f'[t]$? Both of them are for smoothing data, however equation for $f'[t]$ contains weighted average of consecutive points in the array for $f[x]$ itself while the equation for $y[i]$ does not contain consecutive data points for $x[i]$. How can we comprehend this equation as a smoothing of data in $x$?

If this question is not relevant when the equations are taken out of context then I will be happy to provide more details.

The first equation you give is the difference equation for a lowpass FIR filter, or a linear filter with an impulse response that is finite in duration. I'll write it a bit differently (so that it is expressly discrete in time and causal):

$$
f_s[n] = 0.1 f[n-2] + 0.8 f[n-1] + 0.1 f[n]
$$

$f_s[n]$ is the smoothed version of the discrete-time input sequence $f[n]$, generated by passing $f[n]$ through an FIR filter with the coefficients $[0.1, 0.8, 0.1]$. The frequency response of this filter is as follows:

As it turns out, it's not a very good lowpass filter. As the name implies, a lowpass filter should pass low-frequency content while removing higher frequencies. This provides the "smoothing" action that you're looking for, as "jagged", non-smooth features are associated with high frequencies since they change rapidly with time.

Your second equation is an example of a lowpass IIR filter, a linear filter whose impulse response is infinite in duration. The filter's difference equation is:

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

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

As an example, if $\alpha = 0.8$, the frequency response of the filter is as follows:

which is a better filter than your FIR example; it yields much better attenuation of frequencies toward the upper end of the band. Even though it might not be obvious by looking at the difference equation (because of the feedback from the filter output back to its input), it effectively performs smoothing on the input due to its lowpass nature. I'm not sure if this description will be particularly meaningful to you for your application, but these are pretty fundamental signal processing concepts; some study of introductory DSP texts could help fill in the gaps.

Edit: By request, here's a plot that shows both responses on the same axes, illustrating the relatively poor attenuation provided by the FIR example filter:


Is it possible to calculate a z-transform for a filter calculated with Parks-McClellan? I would like to know whether this can be done and if then how would this feat be acomplished in Matlab given a filter calculated with Remez?

h = remez(...);

If this is not possible then I would appreciate an explanation whether this can be done approximately (approximations of the zeros and poles) or whether this is simply mathematically pointless because a broken rational complex polynomial representation does not exist.

If your impulse response is $h(n)$ then the $\mathcal{Z}$-transform of this impulse response (i.e. the transfer function of the filter) is

$$H(z)=\sum_{n=0}^{N-1}h(n)z^{-n}$$

where $N$ is the filter length, i.e. the number of coefficients.

Note that since it is an FIR filter, all its poles are at the origin of the $z$-plane. The zeros are of course the zeros of the polynomial with coefficients $h(n)$.


Coherent Sampling vs Windowing for 50/60 applications Let's say I need to sample 50/60 Hz signals to find RMS. I can use for example 1kHz sampling frequency. That means 100ms sampling is 5 cycles of 50Hz and 6 cycles for 60 Hz. That's coherent sampling. But does this account for frequency variations? I don't have zero cross detection in hardware, maybe if I do it in software (microcontroller) it may help. I've heard that windowing the signal help reduce discontinuities in case I don't exactly get a integer multiple of the frequencies during sampling. Should I use windowing? What other approaches can be implemented?

That's coherent sampling.

Assuming this is a AC power application: it's only coherent sampling if you derive the sampling clock directly from the AC frequency using a PLL (Phase Locked Loop) or equivalent.

But does this account for frequency variations?

Depends on how you derive the sampling clock. If you use a PLL it will track the AC frequency, if you use an independent clock, the frequency will indeed drift.

Should I use windowing? What other approaches can be implemented?

If you really just want to measure the RMS of the signal, you don't need coherent sampling. Actually incoherent is better. You just need to make sure you have enough samples.


BPSK vs $\pi/2$ BPSK, what the benefits? I have come across modulation techniques that have $\pi/2$ phase shifts from the regular constellations.

For example $\pi/2 $ BPSK vs constellation the regular BPSK.

Other examples include $\pi/2 $ QPSK vs. QPSK

I am wondering what are the benefits of phase shifts like this?

Thanks looking forward for your answers

The signal always changes phase at bit transition times. This makes it easier to achieve bit synchronism at the receiver. In plain BPSK, long runs of 0 or 1, with no phase transitions occurring at the boundary
between two identical bits, the phase-lock loop that is generating the
signal telling the receiver where the bit transitions are, can lose lock; or the initial acquisition is more difficult. Also, in satellite channels with nonlinear amplifiers such as traveling-wave-tube (TWT)
amplifiers, smaller phase transitions are better. In particular, if the RF carrier phase makes
a transition of $\pi$ radians (meaning the carrier envelope makes a
transition through the origin), the TWT amplifier output contains
all the sidebands that were carefully suppressed in the uplink. So,
filtering in the satellite is required to suppress these sidebands.
Hence, $\pi/2$-BPSK and $\pi/4$-QPSK are preferred over their plain cousins. Note 4, not 2, with QPSK: $\pi/2$-QPSK does not avoid the
transitions through $0$, and adds complexity without necessarily
adding benefits.

Note: in plain BPSK, the signals are $\pm\cos(\omega_0 t)
= \{\cos(\omega_0 t), \cos(\omega_0t + \pi)\}$, and so,
at bit transitions, the phase either stays the same, or changes by
$\pi$. in $\pi/2$-BPSK, the signal set is $\pm\cos(\omega_0 t)$
for one bit and $\pm\sin(\omega_0 t)$ for the next. Thus,
regardless of whether we have two successive identical bits
or two different successive bit, the signal phase changes by $\pm \pi/2$
since the signal changes from a cosine to a sine.  Figure out for
yourself why $\pi/2$-QPSK has transitions through $0$ while
$\pi/4$-QPSK does not; then think of offset QPSK a little.


Can we really use equal gain combining in the receiver in the real world? Can we really use equal gain combining(EGC) in the receiver in the real world?
Because we have to know the phase of elements in the channel matrix first before we do the EGC,however,the channel should be randomly generated every second,so I wonder that if we can do EGC in the real world or not?

Or we can  do the EGC,but the result is just a little bit better of not doing EGC?

Because we have to know the phase of elements in the channel matrix first before we do the EGC,however,the channel should be randomly generated every second,so I wonder that if we can do EGC in the real world or not?

Yeah, so you continuously have to estimate the channel. You need to do that anyway, no matter whether you use any combining method or not, because phase is typically important to the receiver!

So, think of EGC as a kind of channel knowledge-based beamforming.

Therefore: Yes, EGC is used. The difference to MRC is that there's no need to estimate the receive path's SNR – and depending on the signal and receiver model, that SNR estimate might be non-trivial (e.g. when you don't know the different noise figures of your receivers, or when you can technically easily know the phase of your receivers, but not the strength of the signal, e.g. because it's very weak, before you used combining).


If a K-NN ($ k $ Nearest Neighbors) Algorithm Performs Very Well for Low $ k $, Can Something Be Inferred About the Data Set? I am running kNN  on a very small data set for binary classification. Each class has 100 samples.
I am getting the best performances for $k=1$ and $k=3$. Can I deduce information about my data set from kNN performing very well for small values of $k$?

If the result is consistent with a large test set than it means your training data is dense and well define the degrees of freedom of the problem.

If the training set was small it means there is a good separation between your classes in space.


Why Is Non Linear Least Squares Method from MATLAB and Alglib Gives Different Results on the Same Data? i'm trying to rewrite my Matalab prototype for some DSP to C++ and encountering a displeasing problem. I'm trying to fit data to a function $y = a * (\pi / 2 + arctg(b * x))$. In Matlab it works well with the following code:

fo = fitoptions('Method', 'NonlinearLeastSquares',...
                'Lower', [0, 0],...
                'Upper', [Inf, max(x)],...
                'StartPoint', [1, 1]);
[curve2, gof2] = fit(x,y,ft);

In curve2 I get a and b coefficients, which are pretty well.
I try to reproduce this in C++ with Alglib's lsfit using Levenberg–Marquardt algorithm:

alglib::real_1d_array y, c;
alglib::real_2d_array x;
... // filling x and y arrays with data
c = "[1.0, 1.0]";
alglib::ae_int_t maxits = 0;
alglib::ae_int_t info;
alglib::lsfitstate state;
alglib::lsfitreport rep;
alglib::lsfitcreatef(x, y, c, 1e-4, state);
alglib::lsfitsetcond(state, 1e-5, maxits);
lsfitsetbc(state, "[0.0, 0.0]", "[+inf, 5.0]");
alglib::lsfitfit(state, func);
lsfitresults(state, info, c, rep);

At the output (c) I get completely different coefficeints which are doesn't have any correlation with the same from Matlab.

x and y are the same in Matlab and C++ program. The only difference I see is the optimization algorithm: in Matlab I use trust-region method when in C++ I use Levenberg–Marquardt algorithm.

Could you explain me that strange behaviour?

When you solve Non Linear Least Squares problem of a non convex cost function the end solution (Which is guaranteed to be a Local Minimum) will depend on:

Method of Minimization.
Method Parameters.
Starting Point.

In the case above you set the starting point to be the same for both.
Yet in AlgLib you use the method of Levenberg Marquardt (Classic for Non Linear Least Squares).
MATLAB used to use Levenberg Marquardt as its default in the past.
Yet in recent versions it uses more modern method called Trust Region.
The trust region based methods limit their step size to be more conservative.
Basically like the LM they approximate the problem using Quadratic Function yet in each step the find the minimum with constrain about the domain (The Trusted Region) hence usually ends in the closest local minima while LM might skip it.
If you share the values of a and b in a single simulation we'll be able to look into the Cost Function (As it depends on 2 parameters).
If you are after more modern methods to solve the problem in C (Which I recommend), you should look for Trust Region solvers:

Gnu Scientific Library (GSL) - Non Linear Least Squares Fitting.
Ceres Solver - Solving Non Linear Least Squares.
OptimLib - Conjugate Gradient (It's not optimized for Non Linear Least Squares, yet can be used).

If you look in GitHub for Non Linear Least Squares and Trust Region you'll find more.
Few tips:

The problem isn't convex hence a good strategy is perturbate the starting point. This is a poor man Global Optimization strategy.
If you're sensitive to the result (Must have the global minimizer) you should start with a rough grid search to get better starting points.
At top level if needed use full Global Optimization framework like you have in MATLAB.


Zero-padding the middle of a signal I sample a signal at a certain frequency for a finite amount of time to get a sequence $$(x_n)_{n=1}^N = (x_1, x_2, ... , x_N)$$

with the intention of analyzing its power spectral density by calculating its DFT.

I decide that the tail end of the sequence shouldn't be analyzed (that is, I want to discard all $x_n$ for $n\gt a$ for some $a\leq N$), but I want the DFT to have the same number of samples (in this case, $N$ samples).

One approach would be to truncate the sequence at $x_a$, remove the mean, apply a Hamming window to this truncated sequence, pad the truncated sequence with $N-a$ zeroes, and then calculate its DFT.

Question: Let's say that the section I don't want to analyze is an interval in the middle of the sequence. How do I calculate the DFT of this and make sure the DFT has $N$ samples? Should I remove the mean, apply two Hamming windows (one each on each section I don't want to discard), set the elements in the middle that I want to discard to zero, and calculate the DFT?

This is related to my other question here.

If there's a region in the middle of the sequence that you want to discard, you can't just set those samples to zero and compute the DFT, because the DFT would reflect the fact that your sequence suddenly jumps to the value zero, and then suddenly jumps back again to non-zero values. So the zero values would significantly contribute to the result of the DFT, which is almost the opposite of what you want to achieve.

What you should do is consider the two relevant parts at the beginning and at the end as two realizations of the underlying process, and compute two DFTs, which can then be averaged, just like in Bartlett's method. Of course, each sub-sequence can (and should) be windowed and zero-padded before computing the DFT.


Remove a known wav file from recorded file I looked at many similar questions but still cannot find an answer to my problem.

I'm playing an audio file N, while recording a man talking (b).
The problem is that the microphone gets the played sound N together with the man talking, b. so actually I have y=N*+b in my recording, where N* is the played sound with a short time delay and a bit noise (due to the recording and the phone speaker).

My goal is to totally remove the N* from the recorded file. It seems like an easy problem cause I have the original played file N, but I find it more complicated.

What I tried:

finding the delay between the N and y, invert N and add to y. Seems like it's not working and the result is even worse probably because N* has different amplitudes and a phase shift.
Tried the same thing but subtracting in the frequency domain between y and N. Still no good.

Any suggestions?

Assuming that there are no non-linear effects between the sound N source and the microphone (such as AGC), you might try to estimate the impulse response of the channel between the source for sound N and the microphone data.  The better estimate you can make of that impulse response, the more of sound N that you can remove.  It’s unlikely that the negative of the impulse response is just the data inverted. That's due to the fact that the amplifier, speaker, microphone, and room acoustics (etc.) probably have frequency dependent phase delays and/or multiple reverb/echoes.

Channel estimation in noise is a non-trivial problem.


How to use convolution operator in matlab? I = imread('13.jpg');
  Ir = I(:,:,1);
  Ig = I(:,:,2);
  Ib = I(:,:,3);

   %# Create the gaussian filter with hsize = [5 5] and sigma = 2
   G = fspecial('gaussian',[10 10],1);
  %# Filter it
  Ig = imfilter(I,G,'same');
  %# Display
  imshow(Ig);

  PSF = fspecial('gaussian',7,10);

  Irc = conv2(Ir,PSF,'same');
  Igc = conv2(Ig,PSF,'same');
  Ibc = conv2(Ib,PSF,'same');

Why am I getting this error?

       Warning: CONV2 on values of class UINT8 is obsolete.
       Use CONV2(DOUBLE(A),DOUBLE(B)) or
       CONV2(SINGLE(A),SINGLE(B)) instead.
       In uint8.conv2 at 11
       In Gaussian_Filter_image at 19
      Warning: CONV2 on values of class UINT8 is obsolete.
       Use CONV2(DOUBLE(A),DOUBLE(B)) or
       CONV2(SINGLE(A),SINGLE(B)) instead.

    In uint8.conv2 at 11
       In Gaussian_Filter_image at 20
       Undefined function 'conv2' for input arguments of type 'double'
       and attributes 'full 3d real'.

      Error in uint8/conv2 (line 18)
       y = conv2(varargin{:});

       Error in Gaussian_Filter_image (line 20)
       Igc = conv2(Ig,PSF,'same');

Here I want to perform convolution of image and gaussian filter but I am unable to do that ?

I have even tried the modification after using

   Ir = double(:,:,1);
   Ig = double(:,:,2);
   Ib = double(:,:,3);

then also I am getting the same error?

You are getting the error because image I is of class (data type) 'uint8', and the arguments to conv2 must be of class 'single' or 'double', i. e. floating point. You can get rid of the error by converting I to double.  The best way to do that is to call I = im2double(I), which will re-scale the pixel values to be between 0 and 1.

But looking at the big picture, you should not use conv2 at all, and use imfilter instead, which does take 'uint8' images.  Look at the help for it.  By default it does correlation, but it can also do convolution.


Modulation of real-valued signal I have one doubt according to modulation of a real signal.

there are so many examples where the real signal was modulated with frequency bigger then 0.

Why the signal could not be putted in f=0?

is there connection with a sending power?
( Antenna send bandpass signal $\rightarrow$ signal is matched $\rightarrow$antenna send low P $\rightarrow$power will be reflected, Tx get P again)

It all depends on the transmission medium (the so-called communication channel). Some media do not tolerate components at DC (i.e., zero frequency) or near DC. A typical channel where a carrier frequency greater than zero is required is the radio channel.


Butterworth-like filter with no zero at Nyquist? I have an old digital synthesizer that has a Butterworth-like 2nd order filter with resonance. However, with high frequency tones it appears there is no steep cutoff near Nyquist (I am sure it is the Nyquist as the tones fold beyond that point), as is usual with the bilinear-transformed Butterworth filter. This poses few questions.

Is the attenuation caused by aliasing or something different, and what is it called?
Are there digital alternatives to the Butterworth or bilinear transform that do not suffer from the attenuation near Nyquist, or is this a usual feature?
What filter could an old resource-strained synthesizer have used to achieve a response without the attenuation?

Clarification: By Butterworth-like response it is meant that the response appears similar at lower frequencies, I have not confirmed that there is no ripple or anything like that. Nonetheless, a 2nd order filter of some sort.

Is the attenuation caused by aliasing or something different, and what is it called?

The two standard methods of mapping an analog filter to a digital one are either the bilinear transform or the impulse invariant transform. Neither does a particularly good job at preserving the transfer function at higher frequencies. The bilinear transform will map whatever happens at "infinite" analog frequency to the digital Nyquist frequency. For any analog lowpass filter (Butterworth or not), this will indeed create the steep roll off at Nyquist.

Are there digital alternatives to the Butterworth or bilinear transform that do not suffer from the attenuation near Nyquist, or is this a usual feature?

There are plenty of alternatives. Impulse invariant is one of them, although it to tends overshoot by 6dB at Nyquist. One heuristic method is to start with bilinear, and then move the zeros from z = -1 towards z = 0 until you have something you "like" by whatever criteria is important for your application.

What filter could an old resource-strained synthesizer have used to achieve a response without the attenuation?

A biquad is a biquad is a biqaud. Computational cost is not a function of the specific filter coefficients. There is really no difference between a Butterworth or any other second order filter.


Removing lines of certain length from an image in MATLAB I want to remove the portion of the line (highlighted by the yellow circles)  where pixels form a straight vertical line segment.

If I use the strel

SE = strel('line', 5, 90);
edgesR = imerode(In_R,SE);

then the remaining lines are also being erased.

Is there any other approach here?

This is a program which removes the line segments in a single column of length which is more than the threshold

function [ret_im] = remove_vertical(in_image)

[NZ_r, NZ_c] = find(in_image);

ori = in_image;

%scan through each row of location of non-zero pixels
% If three (threshold) or more pixels are placed in a single column, turn them to zero
thresh = 3;
for i=1:1:size(NZ_r,1)-1

    probe = NZ_c(i); % Set the initial search element
    len = length(find(NZ_c(:)==probe));
    if( len > thresh)
        for j=1:1:len
            if((NZ_r(i+1) - NZ_r(i)) == 1)
                in_image(NZ_r(i),NZ_c(i)) = 0;
            end
        end
    end
end

ret_im = in_image;

end


Solving Wiener Hopf integral equation for causal filter of predictor Given a stochastic signal $x(t)$ with autocorrelation function $R_{xx}(\tau)=\mathrm{exp}(-
\alpha|\tau|)$, $\alpha>0$.

I want to predict $x(t+\lambda)$,$\lambda>0$  by $x(t-\tau)$, $\tau\ge0$ through a linear etimation,

i.e. I want to find estimator $\hat{x}(t+\lambda)=\int_0^{\infty}h(x)x(t-x)\mathrm{d}x$.

By the principle of linear minimal mean squared error, I can get Wiener Hopf equation of the integral form:
$$R_{xx}(\tau+\lambda)=\int^{\infty}_{0}h(x)R_{xx}(\tau-x)\mathrm{d}x $$
$$\implies \mathrm{exp}(-\alpha|\tau+\lambda|)=\int^{\infty}_{0}h(x)\mathrm{exp}(-\alpha|\tau-x|)\mathrm{d}x$$
I tried to use Fourier transform to solve this integral equation for $h(x)$

$$\frac{e^{j\omega\lambda}}{\\\alpha^2+\omega^2}=H(\omega)\times \frac{1}{\\\alpha+j\omega} $$

But I am not sure the right-hand side is correct or not. Is it legitimate to view right hand side as convoution of $h(\tau)$ and $\mathrm{exp}(-\alpha\tau)\mathrm{u}(\tau)$, where $\mathrm{u(\tau)}$ is a unit step function?
The textbook by Papoulis and Garcia use the factorization teqnique (it's on p594 of Papoulis 4th ed and p616 of Garcia 3rd ed) but I cannot see how to factor the transfer function to causal and non-causal part. I guess the answer of $h(x)$ is $c\times\delta(x)$ for some constant. Moreover, I want to find the mean sqaure error of this predictor.

Your Wiener-Hopf equation is correct, but you can't solve it by simply taking the Fourier transform of both sides. The reason for this is that the equation is only valid for $\tau\ge 0$ due to the causality requirement. Taking the Fourier transform of both sides implies that the equation is valid for all values of $\tau$. You have to take another route and proceed as follows.

The error signal is given by

$$e(t)=x(t+\lambda)-(h\star x)(t)\tag{1}$$

We want to minimize the MSE given by

$$\epsilon=E\{e^2(t)\}=R_{ee}(0)=\frac{1}{2\pi}\int_{-\infty}^{\infty}S_{ee}(\omega)d\omega\tag{2}$$

where $R_{ee}(\tau)$ is the autocorrelation function of $e(t)$, and $S_{ee}(\omega)$ is its power spectral density, which is the Fourier transform of the autocorrelation function.

Using $(1)$ it is a rather straightforward exercise to show that $(2)$ can be written as

$$\epsilon=\frac{1}{2\pi}\int_{-\infty}^{\infty}\left[|H(\omega)|^2S_x(\omega)-2\text{Re}\left\{H(\omega)S_x(\omega)e^{-j\omega\lambda}\right\}+S_x(\omega)\right]d\omega\tag{3}$$

The trick is now to rewrite the integrand of $(3)$ by completing the square, which will explicitly show how $H(\omega)$ must be chosen. But first we need to factor $S_x(\omega)$. Since $S_x(\omega)$ is non-negative and real-valued we can factor it as follows:

$$S_x(\omega)=F(\omega)F^*(\omega)\tag{4}$$

where $F(\omega)$ corresponds to a real-valued causal and stable system with all its poles and zeros in the left half-plane (i.e., $F(\omega)$ is a minimum-phase system). Note that $|F(\omega)|=|F^*(\omega)|=\sqrt{S_x(\omega)}$.

Using $(4)$, we can write the integrand of $(3)$ as follows:

$$|H(\omega)|^2S_x(\omega)-2\text{Re}\left\{H(\omega)S_x(\omega)e^{-j\omega\lambda}\right\}+S_x(\omega)=\left|H(\omega)F(\omega)-F(\omega)e^{j\omega\lambda}\right|^2\tag{5}$$

Clearly, we can minimize $(5)$ and $(3)$ by choosing $H(\omega)=e^{j\omega\lambda}$, which just means that we shift the input by $\lambda$ towards the future to get $x(t+\lambda)$, but that's unfortunately not how prediction works. We need to restrict our system to be causal. Note that if $H(\omega)$ is causal then also $F(\omega)H(\omega)$ is causal (because $F(\omega)$ is also causal). This means that the best we can do to minimze $(5)$ and $(3)$ is to use the term $H(\omega)F(\omega)$ to cancel the causal part of $F(\omega)e^{j\omega\lambda}$. So we need to compute the causal part of $F(\omega)e^{j\omega\lambda}$.

From the given autocorrelation function $R_{xx}(\tau)$ we get

$$S_{xx}(\omega)=\frac{2\alpha}{\alpha^2+\omega^2}=\frac{\sqrt{2\alpha}}{\alpha+j\omega}\cdot\frac{\sqrt{2\alpha}}{\alpha-j\omega}\tag{6}$$

We choose

$$F(\omega)=\frac{\sqrt{2\alpha}}{\alpha+j\omega}\tag{7}$$

corresponding to a causal and stable system. The inverse Fourier transform of $F(\omega)e^{j\omega\lambda}$ is given by

$$\mathcal{F}^{-1}\left\{F(\omega)e^{j\omega\lambda}\right\}=\sqrt{2\alpha}e^{-\alpha(t+\lambda)}u(t+\lambda)\tag{8}$$

where $u(t)$ is the unit step function. The causal part of $(8)$ is simply obtained by multiplication with a step function:

$$\sqrt{2\alpha}e^{-\alpha\lambda}e^{-\alpha t}u(t)\tag{9}$$

The Fourier transform of $(9)$ is

$$\mathcal{F}\left\{\sqrt{2\alpha}e^{-\alpha\lambda}e^{-\alpha t}u(t)\right\}=e^{-\alpha\lambda}\frac{\sqrt{2\alpha}}{\alpha+j\omega}=e^{-\alpha\lambda}F(\omega)\tag{10}$$

From $(5)$ we have to choose $H(\omega)F(\omega)$ to equal $(10)$, the causal part of $F(\omega)e^{j\omega\lambda}$. So we finally obtain

$$H(\omega)=e^{-\alpha\lambda}\tag{11}$$

as the optimal prediction filter frequency response, which corresponds to the impulse response

$$h(t)=e^{-\alpha\lambda}\delta(t)\tag{12}$$

So for the given statistics of the input signal, the optimal predictor only uses the current input value to predict the future value $x(t+\lambda)$. It does not use any past values of $x(t)$.


PSK Bandwidth for SNR I'm conducting experiments on several n-ary Psk signals of differing symbol rates. I'm comparing detection levels of different signal detection approaches on sets of PSK signals. I'm using a cosine-roll-off filter with a roll-off-factor of between 0.2 to 0.5.

To determine the appropriate signal to noise ratios, I had to decide which bandwidth I was gonna use. For reasons of practicality, I'd like to use the approximation symbolrate[bd] equals bandwidth[Hz]. It was proposed in chapter 6.2.4 of the book:

Communications and Information Systems; Michael John Ryan,Michael Frater

Do I act reasonably and doesn't it skew my evaluation with inappropriate error levels?

I'm using a cosine-roll-off filter with a roll-off-factor of between 0.2 to 0.5.

Shamelessly copying from wikipedia's Raised-Cosine page, here's a plot of frequency responses of Raised-Cosines filters with different roll-offs:

As you can see, the cutoff frequency for all of these is always half the symbol rate, $\frac12 f_S = \frac12 \frac1{T}$, but the amount of energy outside of that depends on the roll-off. The closer you get to the perfect sinc shape in time ($\beta\rightarrow0$), the sharper that gets.

Notice that this leads to a (double-sided) bandwidth between the cutoff frequencies of $B=2\frac12f_S=f_S$!

However, even for the extreme case of a roll-off factor (which, by the way, is exactly defined to be the difference in occupied bandwidth compared to the minimal bandwidth given by above cutoff frequency) of $\beta=1$, you'd only get $B=2f_S$, so you actually consider the worst case for raised cosine filters – which is a justifiable simplification.

Notice that I'm presuming the book you're citing has a bandwidth definition that covers only a single side!


Bitmap alpha bevel algorithm? I am looking to create an algorithm which adds a bevel effect to a bitmap using it's alpha as a bump map.

How would I go doing something like this? I have tried specular lighting but I only get the highlight and not the shade.

Here is the effect I am talking about(made using Photoshop):

All of these were done using size: 30px(The depth of the bevel from the contour of the bitmap), angle 130, altitude 50.

From left to right, top to bottom:

Chisel Hard bevel
Chisel Soft Bevel
Smooth Bevel
Chisel Hard with soften: 16px - a blurred bevel?

I am trying to create each of these effects, how would I go about creating the basic bevel? and what do I need to get to each of these from that bevel

This can be accomplished with a convolution of the distance transform.

Use a distance transform on the edge of the mask. Then threshold this distance transform to remove values beyond some distance. I think the secret to getting the shading is to convolve the distance transform result with a kernel that looks something like this:

[ -1.0  -1.0  -1.0
  -1.0   0.0   0.0
  -1.0   0.0   1.0 ]

This should get you started in the right direction:

#include "opencv/cv.h"
#include "opencv/highgui.h"

using namespace cv;
using namespace std;

int main() {
    Mat mask, dist, bevel;
    mask = Mat::zeros(200, 400, CV_8U);
    rectangle(mask, Point(30,30), Point(180,180), Scalar(255), -1);
    circle(mask, Point(30,30), 50, Scalar(0), -1);
    circle(mask, Point(180,180), 50, Scalar(0), -1);
    circle(mask, Point(300,100), 75, Scalar(255), -1);
    imshow("1",mask);

    //find edges and invert image for distance transform
    Canny(mask, dist, 50, 150);
    dist = 255-dist;
    distanceTransform(dist, dist, CV_DIST_L2, CV_DIST_MASK_5);
    threshold(dist, dist, 20, 20, CV_THRESH_TRUNC);
    blur(dist, dist, Size(3,3));
    dist.convertTo(bevel, CV_8U);
    equalizeHist(bevel, bevel);
    imshow("2",bevel);

    //convolve with secret sauce
    float d[] = {-1,-2,-3,
                 -2, 0, 0,
                 -3, 0, 1 };
    Mat kernel(3, 3, CV_32F, d);
    kernel = kernel - mean(kernel)[0];
    filter2D(dist, dist, CV_32F, kernel);

    //normalize filtering result to [-1, 1]
    double maxVal;
    minMaxLoc(dist, NULL, &maxVal);
    dist = 128 * dist / maxVal;

    //convert and display result
    dist.convertTo(bevel, CV_8U, 1, 128);
    bevel = bevel.mul(mask)/255;
    imshow("3", bevel);

    waitKey(0);
}


The use of FDM for a single message signal Suppose I have a single message signal, m(t), that is subjected to AWGN upon transmission.

Is there any way that FDM (or any multiplexing strategy) can be used to ultimately improve the SNR of the demodulated signal at the other end of the transmission line. I'm open to using any sort of modulation technique.

Thanks.

You can simply send the message multiple ($N$) times and if all N messages were received at the same SNR you would coherently average the messages for a processing gain in SNR equal to $10\log_{10}(N)$ in dB.  This is effectively trading bandwidth for SNR as you are using more resources to send the same message. To coherently add you remove the complex carrier phase for each message prior to adding in the average. If the messages were not received at the same SNR (such as if it was a fading channel) you would optimally weight each message by the SNR of each message prior to averaging.
This latter point is similar to what occurs in a matched filter in that each sample within a symbol duration is optimally weighted by the SNR for that sample prior to averaging over the symbol duration.

See @MBaz's good comments under the OP's question clarifying that there really is no actual SNR gain if you consider the total signal power of all messages sent, since the total signal power would need to increased to realize the gain listed above.


About the variance of the coefficient of narrowband noise when using signal-space representation Assume that n(t) is a white Gaussian noise process with zero-mean and power spectrum density $N_0/2$. By using the signal-space representation, it can be expressed as:

$$n(t) =\sum_{j=1}^N n_j \phi_j(t)+n_o(t)$$
where $n_j = \int_{-\infty}^{\infty}n(t)\phi_j(t)\,\mathrm dt$ , $\phi_j(t)$ is the basis function and $n_o(t)$ is a part which cannot be expressed in terms of this basis function.

If the noise process $n(t)$ is a narrowband noise of bandwidth $B$ centered on frequency $f$, can it still be expressed in this way?

My understanding is that only after the matched filter, the part $n_o(t)$ is removed.
So the narrowband noise also can be expressed in this way.

Is my understanding correct?

In addition, the coefficient $n_j$ must be a random variable, how to derive the variance of the $n_j$ in the case of narrowband noise ?

Since the basis functions are narrowband, the formula still works and the variance of each $n_j$ remains $N_0/2$..

In other words, assume the bandwidth of the basis functions is $B$. You can then think of them as filters that ignore any frequency outside of $B$. So, they "don't care" if there is noise outside of $B$ or not.

The expression for $n_0(t)$ will change, though.


How to do Finger Tracking which is occluded by an object? Basically I want the hand gesture recognition even when I have an object in hand by using a camera (MS Kinect) in front.

My fingers are not visible directly to the camera due to the object. For example if I could project some electromagnetic rays which can pass through the object but not through the hand so that I get some trace of hand in the image. Even using markers on the fingers can be an option. The material of the object can be plastic or something of this sort.

Any help or reference material would be highly appreciated. Thanks

Use a hand model (using say the kinematics) in order to hypothesize the pose of the articulated movement. For example you can look into:

http://cvrlcode.ics.forth.gr/handtracking/

http://users.ics.forth.gr/~argyros/mypapers/2012_06_cvpr_twohands.pdf

http://robohow.eu/_media/special/bib/isuvr13.pdf

After robustly tracking in presence of the occlusions, you could use a machine learning method such as HMM to recognize the gesture.


just getting started in Signal Processing - easy question I am reading Cycle Analytics for Traders by John Ehlers and need help.

in the first section "Transfer Response" he refers to:

so if I had $4$ prices say $8+8+8+8$ then average would be $8$ and output/input $= 8/32 = 1/4$

but not equal to $1/4 + 8/4 + 8/4 + 8/4$ so I think I've misunderstood?

I don't want to infringe copyright so all I can do is guide you to amazon where you can view this section.

The symbol $z^{-1}$ refers to a delay of one sample interval, no matter how long this interval actually is. It comes from the $\mathcal{Z}$-transform: if $X(z)$ is the $\mathcal{Z}$-transform of a data sequence $x[n]$, then $X(z)z^{-1}$ is the transform of the delayed sequence $x[n-1]$.

The transfer function (i.e. the quotient of the $\mathcal{Z}$-transforms of output and input sequences) is given by

$$H(z)=\frac14\left(1+z^{-1}+z^{-2}+z^{-3}\right)\tag{1}$$

which in the time (sample) domain means that the output sequence $y[n]$ and the input sequence $x[n]$ are related by

$$y[n]=\frac14\left(x[n]+x[n-1]+x[n-2]+x[n-3]\right)\tag{2}$$

This is of course the average of the current and the past three data samples. If, as in your example, all data samples have the value $8$, then the output is

$$y[n]=\frac14\left(8+8+8+8\right)=8$$


Would it be feasible to implement the Sliding DFT on a digital signal processor? The Sliding DFT generates the spectrum for every input sample. Using the FFT, the spectrum is generated only after a certain number of samples (N) are obtained. It appears to me as though the cycles to perform the Sliding DFT would be much larger than the cycles required to do an FFT, to generate an N point DFT. Would it be recommended to use the Sliding DFT to generate an N point spectrum? The motivation for me to consider using the Sliding DFT to generate an N point spectrum is time resolution.

Could anyone give me an idea of the challenges in implementing the Sliding DFT on a DSP?

The paper Sliding DFT control algorithm for three-phase active power filter implements a sliding a sliding DFT on a TMS320F2812 from Texas Instruments, and discusses some advantages, and An Enhanced Interpolated-Modulated Sliding DFT for High Reporting Rate PMUs does similar work on FPGA.

From what I remember, you may have a look at the stability, see  The Sliding DFT.


Decoding a video signal, how to get frame output to matrix for calculations? I'd like to use python to take a video (eventually a stream) and map each frame to an RGB array of pixels.

My first thought was to use ffmpeg bindings, and in memory (if possible) convert it to png to read each pixel. However, I think that this is redundant to the purpose of the decoder.

It seems like this may be the most basic entry into video processing, if someone could point me to the library or approach that would accomplish this I would be grateful.

You can do it with OpenCV library. Here is OpenCV link for Python library.


Point Pattern Recognition Having two different size of sets of points (2D for simplicity) dispersed within two different size squares the question are that:

1- how to find any occurrence of the the small one through the large one?
2- Any idea on how to rank the occurrences as shown on the following figure?

Here is a simple demonstration of the question and a desired solution:

Update 1:
The following figure shows a bit more realistic view of the problem being investigated.

Regarding the comments the following properties apply:

exact location of points are available
exact size of points are available

size can be zero(~1) = only a point

all points are black on a white background
there is no gray-scale/anti-aliasing effect

Here is my implementation of the method presented by endolith with some small changes (I rotated target instead of source since it is smaller and faster in rotation). I accepted 'endolith's answer because I was thinking about that before. About RANSAC I have no experience so far. Furthermore the implementation of RANSAC requires lots of code.

From a computer vision perspective: the basic problem is estimating a homography between your target point set and a subset of points in the large set. In your case, with rotation only, it will be an affine homography. You should look into the RANSAC method. It is designed to find a match in a set with many outliers.  So, you are armed with two important keywords, homography and RANSAC.

OpenCV offers tools for computing these solutions, but you can also use MATLAB. Here is a RANSAC example using OpenCV. And another complete implementation.

A typical application might be to find a book cover in a picture.  You have a picture of the book cover, and a photo of the book on a table.  The approach is not to do template matching, but to find salient corners in each image, and compare those point sets.  Your problem looks like the second half of this process - finding the point set in a big cloud. RANSAC was designed to do this robustly.

I guess cross-correlation methods can also work for you since the data is so clean. The problem is, you add another degree of freedom with rotation, and the method becomes very slow.


Can a (vanilla) Kalman Filter's Observation Matrix $H_k$ depend on the state vector $x_k$? A vanilla Kalman Filter allows for a time varying observation matrix $H_k$. Is it allowable for $H_k$ to be a function of the system state $x_k$ in a vanilla Kalman filter?

First, am I correct that this is not allowed in the vanilla Kalman filter? If it is allowed, then I don't understand how measurements would then act to correct the parameters of $x$ that predict them via $H_k$. I'd appreciate a more explicit description of this than my own hunch though!

Second, am I correct that the Extended Kalman Filter & the Unscented Kalman Filter do explicitly allow for this situation? If I am, do these filters allow observations to correct the terms of $x$ that are predicting them?

The non-linear system model that the unscented and extended Kalman filters work with (omitting control) is the following:
$$x_{k+1} = F(k, x_k, v_k) \\ y_k = H(k, x_k, n_k)$$
where $F$ is the state transition function and $H$ is the observation function.  The linear Kalman filter is just the special case where $F$ and $H$ have the form:
$$ F(k, x, v) = A_k x + v \\ H(k, x, n) = B_k x + n$$ for sequences of matrices $A_k$ and $B_k$.  In particular, the functions $F$ and $H$ are linear.

The following is a bit of a pedantic distinction (particularly in the non-linear case) but while the expression $F(k, x_k, v_k)$ clearly depends on $x_k$, the function $F$ itself does not.  In the non-linear case this is no limitation.  In the linear case, if the matrix $B_k$ depended $x_k$ then the function $H$ would no longer be linear.


What is the relation between eigenvalues and state-space response in control systems? I understand the mathematics behind it but I want to know what happens physically in a real-life system. How do the eigenvalues come into the picture from a non-mathematical (physical) point of view? In linear algebra, a matrix is basically a linear transformation and eigenvectors are those which, on being applied this transformation, just stretches themselves (that is without any rotation, I'm limiting only to 2D space). So when I consider a physical system (I'm only considering a homogeneous or unforced system to make it simpler) what do these eigenvalues or eigenvectors physically represent and how do they relate to their mathematical interpretations?

Just because eigenvalues apply in linear algebra doesn't mean that the same meaning applies everywhere. Think of the Fourier transform: does it only apply to the study of heat transfer? Same here. The A matrix is the companion matrix, which means its eigenvalues represent the roots of the system. For example:
$$\begin{align}
H(s)&=\dfrac{3}{s^2+2s+3}\tag{1} \\
\mathbf{A}&=\begin{bmatrix}{{0}\quad{-3} \\ {1}\quad{-2}}\end{bmatrix}\tag{2} \\
s_{1,2}&=-1\pm j\sqrt{2}\tag{3} \\
\end{align}$$
I'm afraid I don't quite understand what you mean by "from a non-mathematical (physical) point of view". The eigenvalues are the mathematical roots of the system, they tell about the stability, the time and frequency response. There's nothing non-physical about it. The system itself is physical -- as in it has a meaning -- even if only described on paper, such as the example above. It's just that you don't have to think about them from a single perspective, only.


What is the third coordinate in "Projective Transformation" of an image and how to build the matrix of the inverse transform? I am trying to understand the formalism of the projective transform of 2D image.
it has 9 parameters (a-i) which the 9th is redundant since we use houmogenous coordinates.
this transformation preserves only straight lines.

$\pmatrix{x'\\ y' \\ w'}$  = $\pmatrix{ a& b& c\\\ d& e& f \\\ g& h& i}$  $\pmatrix{x\\ y \\ w}$

I have three questions according this formalism:

Since the redundancy (from 8 to 9) of degrees of freedom is $i=1$ ?
what is the $w$ and $w'$ ? are they parameters? how do i calculate them?
How can I find the parameters so I can transform the image to its original(from the black background image to the square one).

what I want to do is build the matrix to transform this image:

and find matrix $M$

so I get back this image:

All the images  above taken from here

What I want finally do is to find the

M = [a  b c;
     d  e f;
     g  h i];

t_proj = maketform('projective',T);
I_projective = imtransform(I,t_proj,'FillValues',.3);
imshow(I_projective)
title('unprojective - rectangular carpet')

taken from here

You have to understand that projective transformation is not linear.
It's $(x,y) -> (X(x, y)/Z(x,y) , Y(x, y)/Z(x,y))$, where X, Y, Z - linear functions.

$w$ and $w'$ are essential here - they help represent this non-linear transformation as linear - matrix multiplication - and that way this linear representation is a homomorphism (preserve "chaining" of operations).$w$ for input image is usually taken as $1$;
You get $w'$ from matrix multiplication ($w' = gx + hy + iw$ ) and the end result of transform is $(x'/w', y'/w')$

Generally speaking you can't just assign $i=1$, but due to homogeneity you can divide all the matrix by $i$ and get new matrix, representing the same transform with right-left element $1$

I would advise you to read some textbook on projective geometry for computer vision to better understand this subject. Those questions have some deep implications.


Time reversal of input signals in convolution Let $f(t)$ and $g(t)$ be signals.

I can't figure out what happens to the convolution of $f(t)$ and $g(t)$ if

both $f(t)$ and $g(t)$ are reversed
one of $f(t)$ and $g(t)$ is reversed

I tried the change of variable in integral, but couldn't find out the answer.

This question appears to involve confusion about the relationship between convolution and cross-correlation. To illustrate the relationship, consider the two rectangular pulses shown below:

The shorter pulse has width = 1 and the longer pulse has width = 2. For simplicity, assume unity heights for the pulses.
In my answer to a previous question, I used a figure from a textbook 1 to illustrate the convolution of $\mathrm {x_1(\lambda)}$ and $\mathrm {x_2(\lambda)}$, denoted by $\mathrm {x_1*x_2}$.
The convolution integral was given by
$$\mathrm {x_1*x_2 =\int_{-\infty}^{\infty}x_1(\lambda)x_2(t-\lambda)\,\mathrm d\lambda} \tag{1}$$
so $\mathrm {x_2(\lambda)}$ was reflected around the ordinate and is the function that gets shifted. This is clearly shown in the figure in my linked answer.
But it is also possible to reflect $\mathrm {x_1(\lambda)}$ about the ordinate and shift it. Then the convolution would be denoted by $\mathrm {x_2*x_1}$ and the convolution integral would be
$$\mathrm {x_2*x_1 =\int_{-\infty}^{\infty}x_1(t-\lambda)x_2(\lambda)\,\mathrm d\lambda} \tag{2}$$
This is illustrated the following figure that I drew in the same fashion as the figure in my linked answer:

Thus $\mathrm {x_1*x_2 = x_2*x_1}$, i.e., it does not matter which function is reversed for the convolution. Of course, if a function is symmetric about its center, then reversing becomes rather trivial. In this example, both rectangular pulses are symmetric about their centers, so the situation is especially simple.
In corroboration of the above, the following (through their equation 1.5-3) is directly quoted from Blinchikoff and Zverev [2, p. 12]:

Let p(t) and q(t) each be a piecewise continuous function for $-\infty < t < \infty$. The convolution of p and q, denoted by $p*q$, is defined to be the third function
$$\mathrm {v(t) = p*q =\int_{-\infty}^{\infty}p(t-\tau)g(\tau)\,\mathrm d\tau} \tag{1.5-1}$$
and the integral on the right is called the convolution integral. Furthermore, the convolution of q and p, $q*p$, is identical to p*q, that is,
$$q*p = p*q \tag{1.5-2}$$
This is easily shown by letting $\mathrm {x = t -\tau}$, $\mathrm {dx = -d\tau}$ in (1.5-1), yielding
$$\mathrm {p*q =\int_{-\infty}^{\infty}p(x)g(t-x)\,\mathrm dx} = q*p \tag{1.5-3}$$

Therefore, in answer to the OP's second question, it does not matter mathematically if f(t) or g(t) is the reversed function, but the calculation may be simplified, in practice, if either function is symmetric about its center.
For cross-correlation, all that is necessary is to shift one function by $\tau$, as per this photo image from [2, p. 305]:

So cross-correlation calculations are quite similar to convolution calculations, but note the important considerations discussed directly below Blinchikoff and Zverev's Table 7-2. I suggest that the OP think about their first question, i.e., what happens when both f(t) and g(t) are reversed and the overlap integrals are calculated.
1  C.D. McGillem, G.R. Cooper, "Continuous and Discrete Signal and System Analysis", 2nd Ed., Holt, Rinehart and Winston, NY, ©1984, pp. 58-59.
2  H.J. Blinchikoff, A.I. Zverev, "Filtering in the Time and Frequency Domains", Wiley-Interscience, John Wiley & Sons, NY, ©1976.


Remove a frequency and all its multiple (harmonics) In a soundfile, I'd like to remove a frequency and all its harmonics. What is the best solution in your opinion ?

Do a loop with 20 notch filters, one after another.
// f is the fundamental frequency
for (i=0; i++; i<21)
{
file.notchfilter(f*i);
}

Question : does each notch filter create some unwanted side-effects, and the fact that we use 20 filters one after another will create very high unwanted artefacts?
Is it possible to design a filter that will remove 20 frequencies in ONE STEP ?
Then we wouldn't need to do a loop, just one pass wil be okay?
Is it called a Comb filter ?
Other method ?

The choice of algorithm depends on your application scenario - i.e. there is no best solution per se. If your focus is on minimum computational requirements, a comb filter will probably be optimum. If you want optimum sound quality, a notch filter based approach could be a better choice.

To gain an overview about advantages and disadvantages of the different methods, maybe our article is of interest. We compared

comb filters
subband comb filters (only processing the lowpass frequencies) and
notch filters

in terms of disturbance reduction capabilities and degradation of the desired signal.

Answers to Your Questions

Yes, each notch filter creates unwanted side-effects. On the one hand, you will get audible ringing artifacts - especially noticeable for narrow bandwidths. On the other hand, higher bandwidths suffer less from audible ringing but bring about increased degradation of the desired signal.
Comb filters remove evenly-spaced harmonics in one step by placing notches that cover the full frequency range. As a consequence, you have no control about how many teeth they have as the notches reach up to the Nyquist frequency. If you want to have control about the affected frequency range, you can split the frequency range into high and low bands, process only one with a comb filter and add both bands together.
Except for the mentioned subband comb filter, none that I know of.


Understanding and Correcting phase shift in General Linear Phase FIR filters I have been learning about General Linear Phase Filters and the four standard FIR implementations.

I was wondering in the case $\beta = \frac{\pi}{2}$ or $\beta = \frac{3\pi}{2}$ if there was a way to correct the phase distortion caused by the filter? More generally, what does the the $\beta$ term actually do when looking at the signal in the time domain? It has the potential to multiply in a complex constant containing an imaginary term but I don't understand what that means.

A generalized linear phase filter with $\beta=\pi/2$ or $\beta=3\pi/2$ introduces (apart from a constant delay) a frequency independent phase shift of $\pm 90$ degrees. Two standard filters that do exactly that are a differentiator with frequency response

$$H(\omega)=j\omega\tag{1}$$

and a Hilbert transformer with frequency response

$$H(\omega)=-j\,\text{sign}(\omega)\tag{2}$$

Note that both frequency responses are purely imaginary and conjugate symmetric, i.e. $H(\omega)=H^*(-\omega)$, as required for real-valued systems.

In the time domain, the filters do the following:

$$y(t)=\frac{dx(t)}{dt}\qquad\text{(differentiator)}\\
y(t)=\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{x(\tau)}{t-\tau}d\tau\qquad\text{(Hilbert transformer)}$$

So in general these systems introduce severe distortion. It is relatively easy to imagine the effect of differentiation. To get a feeling for what the Hilbert transform does, take a look at this table showing some transform pairs.

These two systems can be inverted. The differentiator can be inverted by an integrator, and the Hilbert transformer can be inverted by another Hilbert transformer and by sign reversal.

Note that these filters are just two common examples of systems with imaginary frequency responses, but they are representative in the sense that they introduce severe distortions in the time domain waveform, just like all other systems with imaginary frequency responses with conjugate symmetry.


What's the Underlying Algorithm for Detail Preserving Upsampling in Photoshop CC? The detail-preserving upsampling algorithm in Photoshop CC is truly amazing. It can upscale any image while preserving details in almost real-time. And I truly wonder how it is implemented.

I've started by searching papers with keywords single image super resolution. And seems the current state-of-art algorithms require an trained external database to guess the missing high frequency information. And I believe it is really hard to train & distribute a complete dictionary, right? Also, I don't think it's possible to train a dictionary with image pyramid in real-time.

After that I've turned to the area of image upsampling, which seems fast enough comparing with super resolution algorithms, but there are so many algorithms, I'm totally lost.

So, does anyone have any guesses or hints on how can that tool is implemented? Thanks!

It's really nothing beyond Bi Cubic Interpolation with "Sharpening".

http://www.lynda.com/Photoshop-tutorials/interpolation-settings/124096/140573-4.html

They haven't updated their interpolation algorithms for ages.


Should I Use a High Pass Filter (HPF) or Low Pass Filter (LPF) to Enhance the Features for Image Splicing Detection? I'm recently learning image splicing detection. I've already known image splicing can bring some abrupt changes to an image,which generate the sharp edges.And sharp edges correspond to high frequency. And I want to extract features from images for splicing detection.
In my opinion,I should use a high-pass filter to enhance the edges,but I've read some papers which mentioned using a low-pass filter for gradient detection techniques applied in discontinuity detection.
Now I am very bothered with this.Which kind of filters should I use?I would really appreciate if someone can help me!

Since you're after a filter which emphasizes abrupt changes you should use High Pass Filter.

The issue is you'd be also sensitive to noise.
Hence one way to do it is to apply High Pass Filter on slightly blurred image (Which actually results in a Band Pass Filter).

One easy choice would be using the Gradient of The Gaussian Filter.
Why is it?
Because it is easy to calculate (Analytically) and efficient to implement because of its separability property of the filter.