Difference between sampling rate and baud rate and their relation to channel capacity The sampling rate says that   one Hz can carry at least 2 samples, and the bit rate is obtained by multiplying the sampling rate by the number of bits per sample. Thus, many samples and many bits can be contained in one Hz.  Shannon, on the other hand, says that the maximum bit rate that can be transmitted over a channel is twice the bandwidth. I am confused. Could you please help me?

First let's estabish the context.
Let's say you want to transmit a sequence of numbers $\lbrace a_k \rbrace_{k=1}^{N} = a_1, a_2, \ldots$ Using PAM, this sequence is transmitted with the signal $$s(t) = \sum_{k=1}^N a_k p(t-kT_p),$$ where the pulses $p(t-kT_p)$ form an orthonormal set, that is, $$\int_{-\infty}^\infty p(t) p(t) \text{d}t = 1$$ (which implies that the pulse $p(t)$ has unit energy) and $$\int_{-\infty}^\infty p(t-jT_p) p(t-kT_p) \text{d}t = 0$$ whenever $k \neq j$, which allows the receiver to recover the transmitted sequence using a correlator or a matched filter.
Second, let's see what Nyquist and Shannon say.
First of all, the pulses in $s(t)$ are transmitted at a rate $R_p = 1/T_p$. This is called the "pulse rate" or the "symbol rate". You want $T_p$ to be short, to transmit information faster. However, this increases the bandwidth $B$ of $s(t)$. The relationship is $$B = \alpha \frac{R_p}{2} = \alpha \frac{1}{2T_p}.$$
The value of $\alpha$ is the "spectral efficiency" of the pulse $p(t)$. The minimum possible value is $\alpha=1$, when $p(t)$ is a sinc pulse. In this case, $B=R_p/2$ or $R_p=2B$. This is what it means to say that "$2B$ pieces of information can be transmitted per Hz of bandwidth available."
If this was all there was to it, then an infinite amount of information could be transmitted in any amount of time. All one has to do is let the $a_k$ be arbitrary real numbers, with infinite precision. Even if one is constrained to choose from a constellation, such as $M$-QAM, one could choose $M$ to be as large as needed.
In practice, however, noise, distortion, and other imperfections prevent transmission of $a_k$ numbers with infinite precision. The simplest model of a communications system assumes that the received signal $r(t)$ is corrupted with white Gaussian noise with power spectral density $N_0/2$: $$r(t) = s(t) + n(t).$$
Transmission under the effects of noise restrict the choice of symbols $a_k$: they must belong to a finite set whose elements cannot be easily mistaken for one another. For example, if we constrain the symbols to take values either $-0.5$ or $0.5$, and the noise magnitude is very unlikely to be larger than $0.5$, then transmission can be done very reliably; however, we pay a price in two different ways:

Power: since the symbol amplitudes cannot be made too small, a minimum amount of power is required;
Rate: since we're limited in the choice of symbol amplitudes, we can't transmit as fast as we wish.

In terms of a $M$-QAM constellation, the consequence is that we can't increase $M$ arbitrarily; when $M$ is made too large, the system cannot transmit information reliably any more.
This is where Shannon comes in. He was interested in the relationship between the entropy of the information source that produces the transmitted sequences $\lbrace a_k \rbrace$, the bandwidth, and the signal to noise ratio in the received signal. He definitely did not say that the maximum bit rate is twice the bandwidth; rather, the maximum bit rate is given by the bandwidth, the SNR, and some specific properties of the sequence $\lbrace a_k \rbrace$. The study of these properties is known as "coding theory".
Finally: is it a coincidence that a signal with bandwidth $B$ must be sampled at rate $2B$, and that a signal with bandwidth $B$ can transmit symbols at rate $2B$? It is indeed not a coincidence. You want to transmit symbols at rate $R_p$; in other words, you want a signal that has specific amplitudes at times $kT_p$:

Think of these amplitudes as signal samples:

The blue signal is the transmitted signal $s(t)$, and since its samples have rate $R_p$, then it has bandwidth $B=R_p/2$.


Deriving the Convolution Kernel of the Inverse of a Signal Let $y$ be the inverse (in the sense of convolution) of $x$, i.e.

$$x \star y = \delta$$

Context: $x[n]$ is a discrete signal defined for $n = 0,\ldots, N$. We can assume $x[n] = 0$ if $n \not\in [0,\ldots, N]$.

How to compute the inverse of $x$ (for the convolution)?

\begin{align}(x\star y)[0]&=x[0] y[0] = 1\implies y[0] = \frac{1}{x[0]}\\
(x\star y)[1]&=x[1] y[0] + x[0] y[1] = 0\implies y[1] = -\frac{x[1]}{x[0]^2}\\
(x\star y)[2]&=x[2] y[0] + x[1] y[1] + x[0] y[2] = 0\implies y[2] = -\frac{x[2]}{x[0]^2} + \frac{x[1]^2}{x[0]^3}
\end{align}

I tried to find $y[3]$ manually, but I don't see a general formula appearing...

Since Convolution is a Linear Operation the easiest way to derive the coefficients is using a Matrix Representation of the problem.

Take the kernel $ y $ and build from it the Convolution Matrix.
Then write the vector representation of the Delta Funciton $ d $ and solve:

$$ Y x = d $$

In case the problem is solvable (Namely there is an inverse operator of $ Y $) you'll get you answer.
Of course you can always use Least Squares Solution.

Full code for building Convolution Matrix was given in Estimating a Signal Given a Noisy Measurement of the Signal and Its Derivative (Denoising) and Deconvolution of Synthetic 1D Signals - How To?


Why doesn't the sampling theorem work in this case? In this example, under "Discussion", we can see that the Nyquist sample rate will result in samples with the value zero only. But why is this? I thought the sampling theorem said that if you sample with the Nyquist sample rate or higher, then the original signal can be exactly recovered? But in this case it can't? So does that mean the sample theorem is not actually valid?

The formula given is

$$ F_s > 2F_{\mathrm{max}}$$

not

$$ F_s \ge 2F_{\mathrm{max}}$$

otherwise it is possible that you could be sampling exactly on the zero crossing points. In practice you would choose sometime a bit higher. For example, CD audio uses 44100 Hz to represent signals up to 20000 Hz.


Filtering Frequency Bands Out of a Signal Let's say I have a 2 second data set taken at 220Hz sample rate and I would like to filter out the frequency bands associated with the EEG Spectrum:
$$\begin{align}
\Delta:& [1,3]\text{ Hz}\\
\theta:& [4,7]\text{ Hz}\\
\alpha_1:& [8,9]\text{ Hz}\\
\alpha_2:& [10,12]\text{ Hz}\\
\beta_1:& [13,17]\text{ Hz}\\
\beta_2:& [18,30]\text{ Hz}\\
\gamma_1:& [31,40]\text{ Hz}\\
\gamma_2:& [41,50]\text{ Hz}
\end{align}
$$
What would be the most simple approach to do this ?

What would be the most simple approach to do this ?

I'm a GNU Radio person. So the simplest approach was this GNU Radio companion-designed flowgraph:

Generating this flowgraph file generates a python program that uses GNU Radio to do the signal processing, on as many CPU cores as there are.
You can then run it as
$ ./bandpass_filters.py -f {input file containing float32 samples, one after the other}

To decrease the complexity of all these bandpasses, I first used a low pass to reduce the sampling rate of the signal by half – that's Ok, since you are only interested in frequencies below half of half of your input sampling rate.
I tried this with a 200MB "dummy" file, creating eight 100MB output files – this, to and from a temp directory, took a whole 13 seconds.
I'm pretty hopeful it's fast enough! I benchmarked this a bit, and it seems the source is able to push through about 20 million samples per second; the slow part is the writing to the 8 files.
Note that this is a quick and easy solution – a proper, speed-optimized solution would probably use

decimation, because you're output sampling rate is still 110 Hz, which contains a maximum of 9 Hz (that's a waste of processing power, and hence, speed)
another pre-decimation by 2 through low-pass filtering to 27.5 Hz for the four filters below that rate

However, really, at your modest amount of samples (16GB = 4 Billion 32bit floats?) this isn't really necessary, if you wanted to grab a coffee, anyways.
That flowgraph was really nice and easy to design. Here I am, shamelessly plugging GNU Radio's Guided tutorials if you want to learn to design such signal processing flow graphs yourself; in fact, it's even more fun to use GNU Radio in live flowgraphs, i.e. in systems where the in- or output (or both) are physical devices (microphones, ecg sensors, sonars, radio frontends). Exactly the same flow graph can work with a sound card, if you replace the "file source" by an "audio source" (both come with GNU Radio) and specify a sound card-compatible sampling rate via the -r flag.
Now, since you probably know better than me how to work with Mathematica, here's what I'd recommend you implement:

(first, one time) You use the GNU Radio companion to generate a python program that does the signal processing for you.
You use Mathematica's BinaryWrite() to write your raw samples to a file as 32 bit floating point numbers, e.g. BinaryWrite("/tmp/samples.dat.f32", your_sample_vector, "Real32").
You use Mathematica's abilities to call executables on your machine to execute that python script. Now, I'm no Mathematica expert, so this is all Google-based-guessing: RunProcess({"/path/of/python/executable", "/path/of/generated/python/program", "-f", "/tmp/samples.dat.f32"})
You read (BinaryRead) the resulting /tmp/samples.dat.f32_1_3.dat (and …_4_7.dat etc) in Mathematica.

Optionally, you really consider doing less workflow in Mathematica and use the sheer mass of signal processing things that come with GNU Radio or its module ecosystem.


What is the proper way to filter out unwanted noise from accelerometers? I'm a beginner in signal processing and have a hard time to filter out environmental noise from the floor vibration data collected using accelerometers.

I have tried using butterworth filter (IIR) and sinc function (FIR) along with some window functions. Yet, I do not see obvious change in the time domain plot nor in the frequency domain plot by just using above-mentioned low pass filter.

The frequency domain plot does show some changes when I use the filter along with the window function. But, when I use inverse fft of the filtered data, the magnitude of the original data is reduced. I am not sure if the filter/ window function is applied correctly.

fs = 1652 Hz
fc = 208 Hz

I'm not sure the right way to determine the cutoff frequency since I am not sure what is acceptable range for passband/stopband ripples.

Here is my sinc filter:

def fir_lowpass(data, fc, fs, nfft):

# Convert to normalized frequency
fc_nor = fc / (fs / 2)

N = nfft

n = np.arange(N)

# Compute sinc filter.
h = np.sinc(2 * fc * (n - (N - 1) / 2))

# Compute hanning window.
w = signal.hann(N)

# Multiply sinc filter with window.
h = h * w

# Normalize to get unity gain.
h = h / np.sum(h)

# Apply filter to data
data = np.convolve(data, h, 'same')

return data

Here is my butterworth filter:

def IIR_lowpass(data, fc, fs):

fc_nor = fc / (fs / 2)              # Normalize cutoff frequency

b, a = signal.butter(5, fc_nor, btype='low', analog=False, output='ba')

filt_data = signal.filtfilt(b, a, data)

return filt_data

Assuming:

That you limit yourself to LTI filters.
That you can characterize both the noise and the signal of interest.

Then:

(a) If you want to detect a signal of interest (e.g. detect footsteps), use a matched filter.
(b) If you want to estimate the value of such signal, use a Wiener filter.

These are "the best" you can do (under a bunch of assumptions). You will need to measure and characterize your signals (i.e. get an estimate of the power spectral density of noise and signal for (b), for (a) you need the shape of the signal itself.)


How to decompose a mixed signal into its original signals? I have generated a number of sinusoidal signals and added them together using the following Matlab code snippet:

%% Time Properties
Datapoints = 512;               % Sampling frequency
T = 1/Datapoints;               % Sampling period
Length = 512;                   % Length of signal
t = (0:Length-1)*T;             % Time vector

%% Signal Properties
signal=0;
for amp=1:10
    for freq=2:5:100
        for ph=0:pi/36:pi
            signal=signal+amp*sin(2*pi*freq*t+ph);
        end
    end
end

After that I get back the original frequencies using a fast Fourier transform:

NFFT=Datapoints;
fs=Datapoints/Length;
f = fs/2*linspace(0,1,NFFT/2+1);
fftsignal=fft(signal);
fax_bins=[0:NFFT-1];
[pks,pos]=findpeaks(abs(fftsignal.^2),fax_bins);

Now pos will give me the constituent frequencies.

In addition to that I need to know the corresponding phases and amplitudes to reconstruct the original signals.

How can I get the original phase and amplitude information?

Edit 1 But pks is not giving me the constituent amplitudes.

Edit 2 For debugging I used the following code snippet:

Fs = 20000;
t = 0:1/Fs:0.01;
fc1=200;
fc2=300;
x1 = sin(pi*fc1*t);
x2 = sin(pi*fc2*t);
x=x1+x2;
x=x';
xFFT = abs(fft(x))/length(x);
[amp,freq]=findpeaks(xFFT);

I get amp=0.6252 and freq=2

An FFT will only decompose a mix of sinusoids to the original constituents if all the constituent sinusoids are exactly integer periodic in the FFT's length.  One of your signals is not integer periodic in the number of samples in your FFT's length, so (1) its "leakage" (really windowing artifacts) will contaminate the magnitudes of other signals, and (2) you will need to interpolate (with a Sinc interpolation kernel) between FFT result bins to get a more accurate estimate of its magnitude.


Can spectral density be a complex quantity? I have a signal ($S(t)$) which is product of a Gaussian ($G(t)$) and a random phase function ($e^{i\theta(t)}$, here $\theta(t)$ is a random function), as shown below

$S(t)=G(t).e^{i\theta(t)}$

If I calculate the auto correlation of such a signal ($E[S^*(t)S(t−τ)]$) it turns out to be a complex quantity and the same goes with the spectral density (Fourier transform of the auto-correlation function). My questions are the following.

If the above analysis valid? as the process is not wide sense stationary.
If the analysis is not valid, is there some way to handle this kind of situation?
If the analysis is valid, What does the complex spectral density signifies?

and the same goes with the spectral density (Fourier transform of the auto-correlation function).

No, that's not the case.

Since the autocorrelation is a hermitian symmetric function for any $S$, its Fourier transform is always real.

  If the above analysis valid? as the process is not wide sense stationary.

If the process is not WSS, then you can't just proclaim $E[S^*(t)S(t-\tau)]$ to be dependent on only one variable (usually, $\tau$), and hence, a (1D) Fourier transform doesn't make much sense.

  If the analysis is not valid, is there some way to handle this kind of situation?

Depends! You might want to define/find coherency times and do Short-Time Fourier Transforms within that.

Your system, in fact, is just a phase shifted impulse response – as such, a phase-delay representation, which might be derived from a Frequency Shift-Delay plane, might be more helpful in analyzing things. You'll find such a Frequency Shift-Delay plane in what is called scatter function in wireless communications, representing the Doppler and path coefficients of a wireless channels in motion.

But in your case: Is trying to understand PSD or PSD-equivalents really useful? Don't you just want to build a parametric estimator for $\mu_G$, $\sigma_G^2$ and $\theta$ instead?


Why entropy is undefined in low contrast image I am trying to measure contrast of image by entropy of histogram of image.

Code for computing entropy

float measureContrast_inImage(Mat imagel) {

    Mat hist;
    /// Establish the number of bins
    int histSize = 256;
    /// Set the ranges ( for B,G,R) )
    float range[] = { 0, 256 };
    const float* histRange = { range };
    bool uniform = true; bool accumulate = false;
    /// Compute the histograms:
    calcHist(&imagel, 1, 0, Mat(), hist, 1, &histSize, &histRange, uniform, accumulate);
    hist /= imagel.total();

    Mat logP;
    cv::log(hist, logP);

    cv::imshow("logP", logP);

    //Performs an element-wise multiplication or division of the two matrices
    float entropy = -1 * sum(hist.mul(logP)).val[0];

    cout << entropy << endl;

    //float entropy = 0;

    waitKey(0);
    return entropy;
}

Lena image

But entropy is -nan(ind)

similarly, entropy is undefined for another Image.

Entropy is 5.29657 for another image.

Is there any problem in my code.

I wonder that why in low contrast image, entropy is undefined. Kindly help me to figure out.

Edited

I computed histogram and get entropy by the following code
and my entropy is 256

float entropy(Mat seq, Size size, int index)
{

    int cnt = 0;

    float entr = 0;

    float total_size = size.height * size.width; //total size of all symbols in an image

    for (int i = 0; i<index; i++)
    {

        float sym_occur = (float)seq.at<uchar>(0, i); //the number of times a sybmol has occured

        if (sym_occur>0) //log of zero goes to infinity

        {

            cnt++;

            entr += (sym_occur / total_size)*(log2(total_size / sym_occur));

        }
    }

    cout << "cnt: " << cnt << endl;

    return entr;
}

void imhist(Mat image, int histogram[])
{
    // initialize all intensity values to 0
    for (int i = 0; i < 256; i++)
    {
        histogram[i] = 0;
    }

    // calculate the no of pixels for each intensity values
    for (int y = 0; y < image.rows; y++)
        for (int x = 0; x < image.cols; x++)
            histogram[(int)image.at<uchar>(y, x)]++;
}

Mat histDisplay(int histogram[], const char* name)
{
    int hist[256];
    for (int i = 0; i < 256; i++)
    {
        hist[i] = histogram[i];
    }
    // draw the histograms
    int hist_w = 512; int hist_h = 400;
    int bin_w = cvRound((double)hist_w / 256);

    Mat histImage(hist_h, hist_w, CV_8UC1, Scalar(255, 255, 255));

    // find the maximum intensity element from histogram
    int max = hist[0];
    for (int i = 1; i < 256; i++) {
        if (max < hist[i]) {
            max = hist[i];
        }
    }

    // normalize the histogram between 0 and histImage.rows

    for (int i = 0; i < 256; i++) {
        hist[i] = ((double)hist[i] / max)*histImage.rows;
    }

    // draw the intensity line for histogram
    for (int i = 0; i < 256; i++)
    {
        /*line(histImage, cvPoint(bin_w*(i), hist_h),
        cvPoint(bin_w*(i), hist_h - hist[i]),
        Scalar(0, 0, 0), 1, 8, 0);*/
        line(histImage, cvPoint(bin_w*(i), hist_h),
            cvPoint(bin_w*(i), hist_h - hist[i]),
            Scalar(0, 0, 0));
    }

    // display histogram
    //namedWindow(name, CV_WINDOW_AUTOSIZE);
    //imshow(name, histImage);
    return histImage;
}

int main () {
    Mat image = imread("lena.png", CV_LOAD_IMAGE_GRAYSCALE);
    int hist[256];
    imhist(image, hist);
    Mat histImg = histDisplay(hist, "hist");

    float en = entropy(histImg, image.size(), 256);
}

I believe that it is still wrong.
Please correct me .

You are getting a nan returned because you are performing an operation that yields undefined result.

Based on your code, the only option I can see is the log operater that theoretically yields $-\infty$ when you attempt to take the logarithm of a 0-valued pixel. C++ will thus consider it not a number. Other functions I have seen for calculating entropy such as this example in octave, just ignore the zero pixels to solve these cases.


Jumping from DSP concepts to practice...recommend books? I have fairly good understanding of DSP theory and algorithm design (and also MATLAB and C). What are the best books or what other advice could be given for learning how to actually implement DSP algorithms i.e. DSP programming (e.g. fixed-point/floating-point considerations etc.)

I have also considered getting a DSP evaluation board for learning purposes but am lost as to which one to get, any advice?

The USRP with Gnu Radio software is a good and relatively inexpensive (not that cheap, but much cheaper than almost anything else out there) DSP platform.  Gnu Radio is also open source, so if you develop something good you can submit it to the project.


Trying to understand downsampling and then upsampling I am trying to understand how this works, specifically, what the DTFT of each step looks like in each step of the chain (for understanding). I am not looking for an answer like because the input signal has a frequency of $\pi/3$, it will be reconstructed correctly.

I plotted the DTFT of each step in MATLAB and simply don't understand the result. (DTFT of $x_d[n]$ and $x_r[n]$)

I would also like to know what the transfer functions of the downsampler and upsampler are.

The following pictures shows the DTFT of $x[n]$, $x_d[n]$ and $x_e[n]$ respectively

I don't know where the frequency components in the upsampled signal ($x_e[n]$) are coming from. I have noticed though, that it is what you would get if you convolved the original signal ($x[n]$) with an impulse train with period $2\pi/3$.

PS: Why is there no DTFT tag on this site?

Downsampling can be defined as follows:

$s_{out}[n]= s_{in}[n\cdot M]$ (taking every M-th signal sample)

After applying Z-transform we obtain:

$\bar{s}_{out}[z]=\dfrac{1}{M}\sum_{m=0}^{M-1}\bar{s}_{in}\left[ z^{1/M} \cdot w_{M}^{m} \right]$, where $w_{M}=e^{-2\pi i /M}$

By substituting $z=e^{2\pi i f_n}$ ($f_n $ is the normalised frequency: $f_n = \dfrac{f}{f_s} \in [0, 0.5]$

This yields following relationship between input and output spectrum:

$\hat{s}_{out}[f_n]= \dfrac{1}{M}\sum_{m=0}^{M-1}\hat{s}_{in}\left[ (f_n-m)/M \right]$

You can understand that as a fact that your spectrum is being copied M times (even outside of your original sampling frequency) and stretched by factor of M.

Upsampling on the other hand is governed by following equations:

$s_{out}[n]=\begin{cases} s_{in}(n/M) & \text{n} \; mod \;\text{M} =0 \\ 0 & \text{otherwise}\end{cases}$

After applying Z-transform:

$ \sum_{n} s_{out}[n]\cdot z^{-n} = \sum_{n} s_{in}[n]\cdot z^{-nM}$, which gives: $\bar{s}_{out}[z]=\bar{s}_{in}[z^M] $

By doing same substitution of a unit circle instead of z variable:

$\hat{s}_{out}[f_n] = \hat{s}_{in}[M\cdot f_n] $

You can imagine that this equation is squeezing your spectrum by factor of M, bringing at the same time whatever is outside of your sampling frequency range.

Plots you've obtained are great example how it works!

P.S. There is a dft tag.


How to include phase in a sinusoidal Kalman Filter I start with the equation for sinusoidal motion with an offset and differentiate to get the 2nd order ODE describing the motion of the object.

\begin{align}
x &= A\sin(\omega t + \phi) + O\\
\dot{x} &= A\omega \cos(\omega t +\phi)\\
\ddot{x} &= -\omega^2 A \sin(\omega t +\phi)\\
\implies \ddot{x} &= -\omega^2(x-O) = -\omega^2x + \omega^2O
\end{align}

In matrix form this becomes:

$$
\begin{bmatrix}
\dot{x} \\
\ddot{x} \\
\dot{O}
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 & 0\\
-\omega^2 & 0 & \omega^2 \\
0 & 0 & 0 \\
\end{bmatrix}
\cdot
\begin{bmatrix}
x \\
\dot{x} \\
O
\end{bmatrix}
$$

$\implies$ the system dynamics matrix $\mathbf F$ is:

$$
\begin{bmatrix}
0 & 1 & 0\\
-\omega^2 & 0 & \omega^2 \\
0 & 0 & 0 \\
\end{bmatrix}
$$

then by performing $$
A(t) = \mathscr{\mathbf L}^{-1}\left(s\vec{\mathbf I} - \vec{\mathbf F}\right)^{-1}
$$

I get the following discretised $\mathbf A$ matrix:

$$
\begin{bmatrix}
\cos(T_s \omega) & \dfrac{\sin(T_s \omega)}{\omega} & -\cos(T_s \omega) + 1\\
-\omega \sin(T_s \omega) & \cos(T_s \omega) & \omega \sin(T_s \omega) \\
0 & 0 & 1 \\
\end{bmatrix}
$$

This works to filter a sine wave of a know frequency with an offset, it predicts the correct phase, amplitude and even if a change the amplitude of the signal being filtered half way in it corrects for this very quickly and fits to the true signal. However if I shift the phase of the true signal part way through it does not correct for it. How might I include phase $\phi$ in my Kalman filter model such that it can still track the signal with sudden changes in phase?

If I try with my current filter I get the following result:

The answer is that the filter I detailed can correct for the phase, the issue was that the values in my Q matrix were too small and therefore it was overly trusting of the model and was correcting for the change in phase very slowly. By increasing the values in the Q matrix this model can quickly correct for the change in the phase of the oscillator!

The output now looks like the following:


Calculating the network function of a filter Given a system, that behaves as a 1st order filter with network function $H(s)$. We input:
$$v_1 (t)=1+3\cos(10^4 t)$$

And we obtain as output:
$$v_2(t)=1+1.5\cos \left(10^4 t -\dfrac{\pi}{3}\right)$$

Say what kind of filter it is and find its network function $H(s)$.

I'm trying to solve this problem and I have solved the first part by saying that the filter is a low-pass filter because the continuous term is conserved in the input as well as in the output. And I have supposed that the network function will be of the form:
$$H(s) = \dfrac{\kappa}{s+a}$$

My question is what do I have to do in order to compute the constants $\kappa, a$ according to the filter characteristics?

I have tried to compute the amplification for the following frequencies:
$$\text{Amplification}(\omega=0) =\dfrac{1}{1} = 1$$
$$\text{Amplification}(\omega=10^4) = \dfrac{1.5}{3} = 0.5$$

But I don't know if this is correct since the phase of cosine in the output is different.

You need to consider the system's frequency response

$$H(j\omega)=\frac{\kappa}{j\omega+a}\tag{1}$$

Now you know that $$H(0)=1\tag{2}$$ and $$|H(j\omega_0)|=\frac12\tag{3}$$ (with $\omega_0=10^4$). Note that $H(0)$ is real-valued, whereas $H(j\omega_0)$ will generally be complex-valued. From (1) and (2) you immediately get $\kappa=a$. Furthermore, since the filter is stable we know that $a>0$ (i.e. the pole must lie in the left half-plane of the complex $s$-plane). Combining (3) with (1) gives

$$\frac{a}{\sqrt{\omega_0^2+a^2}}=\frac12\tag{4}$$

from which it should be easy for you to compute $a$.

Note that you don't need the phase value to specify the system. What should be checked is if the phase value is the actual phase value of the system that you just computed. In order to do this you need to verify that

$$\arg\{H(j\omega_0)\}=-\frac{\pi}{3}\tag{5}$$

Luckily that's the case, otherwise the assumption (1) would be wrong. I leave the proof of (5) up to you.


Derive DFT of $x((n+1)/2)$ If $X(k)$ is the $N$-point DFT of $x(n)$, and
  $y(n)= x\left(\frac {n+1}{2}\right)$ for odd $n$, and $0$ for even $n$.

  What is the $2N$ point DFT of $y(n)$ in terms of $X(k)$?

So far, I've noticed that $y(2n+1)=x(n+1)$, but I'm not sure what to do with that information in terms of the DFT summation. Can someone help me out here?

I would follow a signal block diagram based solution for this problem.

First as suggested in the comments, it's very helpful to investigate a few values of $y[n]$ and $x[n]$ for some $n$ :

$$
\begin{align}
y[0] &= 0 ~~~,~~~ y[1] = x[1] \\
y[2] &= 0 ~~~,~~~ y[3] = x[2] \\
y[4] &= 0 ~~~,~~~ y[5] = x[3] \\
y[6] &= 0 ~~~,~~~ y[7] = x[4] \\
\end{align}
$$

with some experience, or by a direct investigation of this above sequence it can be seen that the system that produce $y[n]$ from $x[n]$ is the folowing :

$$ x[n] \longrightarrow \boxed{ \uparrow 2} \longrightarrow  w[n] \longrightarrow \boxed{ z^1} \longrightarrow y[n]$$

where the up arrow indicates an expansion by $2$ and the $z^1$ indicates a left shift (advance) by one sample.

Write down the DTFT relations between those signals $x[n]$,$w[n]$ and $y[n]$ :

$$
\begin{align}
W(\omega) &= X(2\omega) \\
Y(\omega) &= e^{j\omega} W(\omega) \\ \\
Y(\omega) &= e^{j\omega} X(2\omega) \\
\end{align}
$$

And relate $2N$ point DFT $Y[k]$ of $y[n]$ to $N$ point DFT $X[k]$ of $x[n]$.

$$
\begin{align}
Y[k] &= Y(\omega_k) = Y(\frac{2\pi}{2N} k) &, ~~ k=0,1,...,2N-1\\
Y[k] &= e^{j \frac{2\pi}{2N} k} X(\frac{2\pi}{N} k) &, ~~ k = 0,1,...,2N-1\\
Y[k] &= e^{j \frac{\pi}{N} k} X[k] &, ~~ k = 0,1,...,2N-1
\end{align}
$$

It can be seen that the length $2N$ DFT sequence $Y[k]$ consists of two fold replica of length $N$ DFT sequence $X[k]$ weighted by $e^{j \frac{\pi}{N} k}$. Note that $X[k]$ is periodic in $k$ by $N$, hence repeats twice during $k = 0,1,...,2N-1$.


Does the output signal from the FIR filter have transient response? In matlab, I have shown the signal passed through the FIR filter in the time domain. Then the beginning of the signal is like the transient response. Is the result appropriate?

Probably. Except for 1-tap filters, if your signal is of finite support, FIR filtering needs to assumes  values outside the support (before and after the values you were given)), and those values are often considered to be $0$. Hence, the extended sequence $(\ldots,0,0,x[0],x[1],x[2]\ldots)$ may produce some artifacts or transients in the output around time $[0]$, depending on the continuity and smoothness of $x$ close to $k=0$.


Super Resolution in Frequency Domain Using Compressed Sensing To be noted that I'm very new to this topic, I would like to understand the fundamentals of how to get Super Resolution in Frequency Domain estimation using the Compressed Sensing Model.
I am also looking for some references and Python/Matlab code that can help me.
Thanks a lot in advantage and happy new year,
Luca

You can employ Compressed Sensing / Sparse Representation for Super Resolution in Frequency Domain.
One way to do so is solving the problem:
$$ \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| F \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| \boldsymbol{x} \right\|}_{1} $$
Where the $ {L}_{1} $ norm is sparsity inducing regularization and $ F $ is the inverse DFT matrix.
Solving this, quite simple, optimization problem will yield Super Resolution of the DFT.
Super Resolution means, in that context, being able to resolve frequencies which are closer than what the observation time suggests:

In the above you can see the DFT of a sum of 2 sines with the given frequencies. The Gaussian model is using $ {L}_{2} $ for regularization (Which is basically damped zero padding).
You may see that the $ {L}_{1} $ could resolve the 2 sines even when they are only 0.5 [Hz] apart with an observation windows of 1 [Sec].
This is pretty nice...


Computing Discrete Convolution in terms of unit step function Consider :
$$
x[n]:=\begin{cases}
1&\text{if $3\leq n\leq 8$}\\
0&\text{if otherwise}
\end{cases}
\quad\text{and}\quad
h[n]:=\begin{cases}
1&\text{if $4\leq n \leq 15$}\\
0&\text{if otherwise}
\end{cases}
$$
I wish to compute the convolution $\displaystyle y[n]=(x*h)[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]$. Because I dislike the method of graphing, I wanted to take another approach using unit step function :
$$
x[n]=u[n-3]-u[n-8]
\quad
\text{and}
\quad
h[n]=u[n-4]-u[n-15]
$$
By using distributivity property of the discrete-time convolution, we get :
\begin{align*}
(x*h)[n]&=(u[n-3]-u[n-8])*(u[n-4]-u[n-15])\\ \\
&=\underbrace{u[n-3]*u[n-4]}_{(y_{1}[n])} -\underbrace{u[n-3]*u[n-15]}_{(y_{2}[n])}-\underbrace{u[n-8]*u[n-4]}_{(y_{3}[n])}+\underbrace{u[n-8]*u[n-15]}_{(y_{4}[n])}
\end{align*}
Therefore we have four convolutions for which we can see that :
$y_{1}[n]$ = $\begin{cases}\displaystyle\sum_{3}^{n-4}1&\text{if $7\leq n \leq 11$}\\0&\text{if otherwise}\end{cases}=\begin{cases}n-6&\text{if $7\leq n \leq 11$}\\0&\text{if otherwise}\end{cases}$
$y_{2}[n] = \begin{cases}\displaystyle\sum_{3}^{n-15}1&\text{if $18\leq n \leq 33$}\\0&\text{if otherwise}\end{cases}=\begin{cases}n-17&\text{if $18\leq n \leq 33$}\\0&\text{if otherwise}\end{cases}$
$y_{3}[n] =\begin{cases}\displaystyle\sum_{8}^{n-4}1&\text{if $12\leq n \leq 16$}\\0&\text{if otherwise}\end{cases}=\begin{cases}n-11&\text{if $12\leq n \leq 16$}\\0&\text{if otherwise}\end{cases}$
$y_{4}[n]=\begin{cases}\displaystyle\sum_{8}^{n-15}1&\text{if $23\leq n \leq 38$}\\0&\text{if otherwise}\end{cases}=\begin{cases}n-22&\text{if $23\leq n \leq 38$}\\0&\text{if otherwise}\end{cases}$
I must note that I am afraid I don't know if the upper bound in the sum of $(1)$, $(2)$, $(3)$, and $(4)$ are right or wrong. I think what's left to do is compute :
$$
y[n]=y_{1}[n]-y_{2}[n]-y_{3}[n]+y_{4}[n]
$$
I hope someone can show me how from here we can prove that the answer is :
$$
y[n]=\begin{cases}
n-6&\text{if $7\leq n\leq 11$}\\
6&\text{if $12\leq n\leq 18$}\\
24-n&\text{if $19\leq n \leq 23$}\\
0&\text{if otherwise}
\end{cases}
$$

As a first remark, I would really suggest to also familiarize yourself with the graphical way of convolving two functions or sequences. This will also make your calculations a lot easier because you get some intuition on the nature of the solution.
Second, your definitions of $x[n]$ and $h[n]$ in terms of the step sequence $u[n]$ are not correct. E.g., the given $x[n]$ equals $1$ for $n=8$, but $u[n-3]-u[n-8]$ is zero for $n=8$.
It is formally correct to represent the solution as the sum of four individual convolutions of shifted step sequences, so what you could do is derive a general solution of the individual convolutions of the form
$$u[n-n_0]\star u[n-n_1]=(n-N_0)u[n-N_1]\tag{1}$$
I leave it up to you to derive the correct values of $N_0$ and $N_1$ in terms of $n_0$ and $n_1$. Once you have obtained that result, you can just compute the solution to the original problem as the sum of four terms of the form $(1)$.


effect of high pass filtering and low pass filtering on background color Consider we have an image like this:

After applying a Low Pass Filter (say Butterworth), the result is :

And after applying a High Pass Filter (Butterworth), the result is:

The text book (Gonzalez) explains everything clearly. But I still cannot catch the point why the background in LPF is white(gray) and in HPF is black. according to the book, it is because of convolution that averaging a neighbor pixel to calculate the center pixel. I know how convolution works, sliding a small kernel, but why the result from LP filtering has different background in comparison with HP filtering?
I think for HPF, it is because of High Pass Filtering in frequency domain is analogous to differentiation in the spatial domain? so all constants become zero (black)? what about Low pass filtering?
All pictures are from Gonzalez's Digital Image Processing,3rd Edtion,2008.

It's very simple.  White = 255, Black = 0.  After high pass filtering, the pixels that don't change relative to their neighbors (low spatial frequency) get filtered out--meaning the pixel value is reduced.  Thus, any solid color areas will be replaced by blackness, and only points of pixel value transition will be white.


Confused about the Fourier Transform of $e^{at}u(t)$ This is the problem at hand:

I'm unaware of why we didn't have to say anything about $\omega$ like that it should be also greater than $0$, I know it's variable...but it's multiplied by $t$ whatsoever.
Also, if $a$ is complex, I don't understand why do we require that $Re \{a\} >0$ only, without mentioning anything about the imaginary part (it's mentioned further down in the book).
Finally, why isn't it $Re\{a\}>=0$? is it because that may reduce to the Heaviside step function if the imaginary part of $a$ is also 0?
Thank you.

$\omega$'s value does nothing for the convergence of the integral. $a$'s does.

$a$ is not complex, it's real only.

If $a$ is 0, the integral doesn't converge.


DFT without knowing the exact value of sampling rate? I have a sequence of (real) numbers that represent the magnitude of a certain natural event. I know that the samples are equispaced in time, but not the exact value of the spacings. So does that mean I do not know the sampling rate? And, how should I proceed on to doing a DFT on this dataset?

P.s. Let's say the dataset represents a historical record of the heights of high-tides on an unknown planet with one or more natural satellites. I know that the data points are equispaced, but I don't know the exact value of the spacing.

If you don’t know the sample spacing, you don’t know the sample rate.

If you know that samples are uniformly spaced, you can do a DFT and get relative frequencies. You won’t be able to map one dft bin to a specific frequency (in Hertz) unless there is some additional side info (perhaps a known reference hidden within the data).


DFT Vs DCT - spectrum differences When you are transforming a signal in to a frequency domain using both DFT and DCT, say for a function sin(x), how the spectrum will be? will both frequency values (of DCT and DFT) same? Can anyone explain the values of the frequencies for the function sin(x) (obviously the frequency should be 1/2π)

The DCT for a signal $x(t)$ on $[0,T]$ is roughly the same as the DFT for the function $x(|t|)$ on $[-T,T]$, that is, mirrored at the vertical axis. Since that function is even, all sine coefficients will be zero and can be removed from the result.

Now $\sin(|t|)$ is different from $\sin(t)=-\sin(-t)$ for negative $t$, so the resulting coefficients will be quite unrelated. If the sampling interval ends at $t=T=2\pi$, the DFT will only exhibit one peak around frequency $f=\frac1{2\pi}$.

But because of the kink at $t=0$ in $\sin(|t|)$, and a second at the end of the sampling interval, the DCT will exhibit a general background of size $f^{-2}$ in the coefficient for frequencies $f>0$.


Laplace transform of product of signal and impulse train I'm reading 'Discrete Time Control Systems' book by Ogata and came across a few statements (specifically, (3-1) and (3-2)) which I have not been able to understand.

It is said that any continuous signal can be sampled and the output represented as
$$y(t) = \sum_{n=- \infty}^{+\infty}x(nT)\delta(t-nT) $$

Now taking laplace transform
$$\begin{align}
Y(s) &= \sum_{n=- \infty}^{+\infty}x(nT)\mathscr{L}\{\delta(t-nT)\} \\
  &= \sum_{n=- \infty}^{+\infty}x(nT)e^{-nTs} \\
\end{align}$$

Now I have a confusion:

Is the $\delta(t)$ function

the dirac delta function, so that $\mathscr{L}\{\delta(t-nT)\} = e^{-nTs} $ but then the signal representation makes no sense as there is infinite amplitude in the output signal at multiples of $nT$
or is it the unit impulse function (value $1$ at $t=0$ and value $0$ everywhere else) in which case how exactly has $Y(s)$ been evaluated?

$\delta(t)$ is indeed the Dirac delta impulse, which is not defined by its values (because it is no ordinary function), but by its properties under an integral.

The first expression in your question is a standard model for sampling a continuous-time signal, and since it's a mathematical model it does not represent actual physical sampling. You can replace the Dirac delta impulse by other impulse-like functions. If you use a rectangular impulse you get a zero-order hold. However, if you're not interested in any effects of non-ideal sampling, multiplication with a Dirac comb is a convenient and useful model for sampling.


Balance brightness and contrast between two images I have two different images for the same object. Each one was taken from different device (scanner, mobile camera, professional camera ..etc).
the process of  registeration one of them to the other was done correctly. Now I need to do a subtraction between to them for comparison purpose. The problem that is the contrast and brightness are pretty different. I tried these two methods (I tried them on the whole image(globally) and on pairs of the two images (locally)):

Gray-Scale method:

ImageB= ImageA * (Mean_Of(ImageA)/Mean_Of(ImageB));
Color method:
Convert it to HSV then:

ImageB[S]= ImageA[S] * (Mean_Of(ImageA[S])/Mean_Of(ImageB[S]));
ImageB[V]= ImageA[V] * (Mean_Of(ImageA[V])/Mean_Of(ImageB[V]));

return to RGB again

However, the two methods did not give me the results I want. They are good and work well on many images but not in all cases.

Any suggestions?
Many thanks

The images you are working with are probably with gamma compensation, you should remove it first by applying inverse gamma.
More accurately, each camera has its own tone curve mapping. You can build a lookup table that converts one from another.
After you removed the non-linearity source, (whether you used method 1 or method 2), as a first order approximation, you can map colors from different cameras by using a 3x3 matrix applied on RGB.

Now you can compare your images.


How to manipulate attenuation in Parks MClellan algorithm? I am using GNU Radio C++ source code in my app (pm_remez()) for generating FIR filters coefficients.
Input parameters here are: bands[], des[] and weight[].
How I can manipulate stop band attenuation in pm_remez() function?

The stopband attenuation is controlled indirectly by the relation between passband and stopband weights. A large weight means a small error in the corresponding frequency band. E.g., if the passband weight is $1$ and the stopband weight is $10$, then the approximation error in the stopband will be $10$ times smaller than the approximation error in the passband. The actual value of the error (i.e. passband ripple and stopband attenuation) cannot be prescribed in the original version of the algorithm. The approximation error is an outcome of the algorithm, just like the filter coefficients. The only thing that can be prescribed is the ratio between the errors in the different frequency bands.


Power spectrum of 2D image: result interpretation Let's say I have an image, which is NxN pixels in size. If I calculate the fourier power spectrum of this image, I get NxN values, with the highest distinguishable frequency at +- N/2 in each direction.

Let's also say I wanted to get a histogram of power over frequency, because I want to know which frequency carries which amount of power.
Then I'd create concentric circles, with the origin of the fft. Bigger radius (and bigger distance from the origin), higher frequency, etc.

But now the questions start:

If the input image is not a square, but a rectangle, the highest frequency still has to be N/2 and M/2, respectively, and again both + and -.

In my opinion, it does not make any sense to still use circles for calculating the histogram. I'd have to use ellipses, wouldn't I?

Imagine a circle/ellipse now, going in each direction to the max frequency.
Then there are four areas (area of the square minus the area of the circle) in each corner, from which I don't get any contribution to this histogram, because the absolute distance is larger than my max frequency.

Where's the error in that?

And last, I read several books and articles in the past few weeks. I read a passage about the same problem, but they didn't use circles/ellipses, they used a crooked square, i.e. a rhombus, with each corner on one of the axis.
The problem is, I can neither find the passage again nor remember the actual text (for full text search).
Does this approach make any sense? What about the four triangle-shaped areas in this case?

Regards and thanks for your help.

it does not make any sense to still use circles for calculating the histogram. I'd have to use ellipses, wouldn't I?

That's right. For radial averaging on a rectangular grid, the locus of equidistant points from the origin would look like an ellipse instead of a circle.

  there are four areas (area of the square minus the area of the circle) in each corner, from which I don't get any contribution

That's right, and there's nothing we can do about that. If you care about those high frequencies, you should re-acquire your image at a higher sampling rate.

  they used a crooked square, i.e. a rhombus, with each corner on one of the axis. ... Does this approach make any sense? What about the four triangle-shaped areas in this case?

In this case the locus represents points where the sum of the two spatial frequencies is a constant, i.e. instead of $\sqrt{f_1^2 + f_2^2} = \mbox{constant}$ we follow the points where $|f_1| + |f_2|=\mbox{constant}$. There's no reason why one is better or worse than the other. They are both valid methods of collapsing a 2D function of two spatial frequencies into a 1D function of one "combined" spatial frequency.

As with the circle/ellipse, there's nothing we can do about the bins outside the triangle area that never get analyzed.

When using a DFT, distance can be measured as the number of steps taken on a discrete 2D frequency grid. (See the discussion at the bottom of the page here: http://iubemcenter.indiana.edu/beyondNyquist.html). Constant frequency "ellipses" in this case are really "rectangles" of grid cells centered around the origin. This way all grid points get accounted for and there are no "dead" zones when computing the 1D collapsed PSD from the 2D DFT.


Why do we discard imaginary part of the phase spectrum? Suppose I compute phase spectrum from the fftn function in MATLAB as

image1=imread('cameraman.tif');
figure,imshow(image1);

fourier_transform=fft2(image1);%take fourier transform of gray scale image

phase_spectrum=exp(1j*angle(fourier_transform));
imshow(phase_spectrum, [])

Here, by using imshow(phase_spectrum,[])  we first
discard the imaginary part and then normalize.

But why do we discard the imaginary part of the phase spectrum ? Whether it's spectrum  is same as real part with 90 degree phase shift ?

Is it not useful ? Also,what should be done to get whole phase spectrum not just real part of it?

Phase spectrum by itself is a real function of frequency and, therefore, does not contain any imaginary parts. It is the fourier transform which is in general a complex function of frequency, henceforth, may include real and imaginary parts. But neither the magnitude spectrum nor the phase spectrum contain(s) any imaginary parts.

Phase spectrum is a real function. Whatever methods/functions (correct or not) you use to obtain it does not and shall not change this fact. $H(e^{j\omega}) = |H(e^{j\omega})| e^{j\phi(w)} $ where the $\phi(w)$ contains the phase spectrum and the absolute value part is the magnitude spectrum. Both are real. You shall use a proper function and method in your environment to obtain those signals from $H(e^{j\omega})$


Basic Understanding of Convolution Operation I would really appreciate any help in understanding this.

If $$y[n] = x[n]*r[n]\tag{1}$$ then what is $y[n-n_0]$ for any system?

Is it $$y[n-n_0] = x[n-n_0]*r[n-n_0]\tag{2}$$
or    $$y[n-n_0] = x[n-n_0]*r[n]\tag{3}$$
or    $$y[n-n_0] = x[n]*r[n-n_0]\tag{4}$$

Will that result (answer to the above) change when we now have an LTI system with $h[n]=r[n]$ being its impulse response?

Thank you very much

Since $y[n-n_0] = y[n] \ast \delta [n - n_0]$ and convolution is both associative and distributive, is it not the case that

$$
y[n-n_0] = \left(x[n] \ast r[n] \right) \ast \delta[n-n_0] = x[n] \ast r[n-n_0] = x[n-n_0] \ast r[n]?
$$


Find Magnitude of a DFT signal using a Blackman window Hi I've been given a signal made of a series of cosines. I have taken the DFT of the signal using a rectangular window (blue), hamming window (red) and a Blackman window (black). I have identified two of the cosine components of the signal. I was wondering how can I use the dB magnitude to find the actual magnitude of the frequency. This is so that i can recreate the signal by summing all of the cosines found.
I was also wondering why these frequencies can be seen on the blackman window but not the other 2 windows?
any help would be greatly appreciated!

First of all, this is a (y-axis) logarithmic plot, so the plotted magnitudes are given by :
$$ Y_{dB} = 20 \cdot \log_{10} ( X ) ~~~,~~~ X > 0 .$$
And the linear amplitude is found by inverting it :
$$ X = 10^{Y_{dB}/20}  .$$
There will also be an amplitude scale due to DFT length and the window type, which you can determine by a calibration procedure (if not with a formula). You can also find the formula, or just approximate it.
Coming to the reason why Blackman window was able to show the side-by components of lower amplitudes: among other windows there, Blackman window has the largest main-lobe-width, but the smallest side-lope-peaks; i.e., fast decaying tails.
Large main lobe width, means low frequency resolution : if you have two sine waves very close in frequency together, then they will look like a single component instead of two. Hence spectral resolution (discrimination of two close components of similar amplitude) is reduced.
Small side lobe peaks (or fast decaying tails), means reduced leakage: When there's a mixture of large and small amplitude sinewaves (at arbitrary frequencies) in the signal, then the spectrum of the largest amplitude signal will shadow the spectrum of the remanining smaller amplitude components, especially when they are closer to the main component. All those components falling below the (modulated) window's tail magnitude will be invisible. Hence a window type with small tail magnitude (fast decaying tais) will have less shadow on the remaning spectrum. The effect is more pronounced as you get closer to the main lobe of the largest components (as the tail amplitude will be larger around the main lobe).


What is 'energy compaction' in simple terms? This is my understanding of 'energy compaction' and want to know if it is right.  Take a vector.  The energy of the vector is the sum of the squares of its elements.  If A is the transformation matrix that is unitary, it can be proved that the energy in x and Ax are same.  Energy conservation property.  Energy compaction means that the energy of Ax=y is more concentrated in some elements compared to the distribution of energy in x.  DCT is said to have energy compaction property.  Does that mean, for any x, if A is DCT matrix, energy of y will be more concentrated when compared to the x.  Does this happen to every x or x has to satisfy some properties to get this energy compaction?

Yes, I believe that your understanding of energy compaction is correct.

  Does that mean, for any x, if A is DCT matrix, energy of y will be
  more concentrated when compared to the x.

No, it does not.  All that is needed to prove that such is not the case is to show that there is an $x$ that is not compacted by the DCT.  White noise, for example, would not be compacted by the DCT.  The DCT is useful because in many real life situations the "signal" (e.g. audio, images, videos, etc.) tend to be "pinkish", i.e. tend to have most of their energy in the lower frequencies and so there is a natural assymetry that can be used to our advantage.

Another way to look at this is from an information theory perspective.  If a signal is not "white" (i.e. it doesn't have the same power at all frequencies) that implies that there is some correlation between samples in the time domain, which means that the samples have "information" about the value of other samples.  This mutual information implies that there is redundancy, and thus should be able to reduce the amount of data without losing information.

  Does this happen to every x or x has to satisfy some properties to get this energy compaction?

For frequency-based transforms like the DCT it is clear from the above that what is needed is for the signal to be non-white.  The more non-white it is, the more compaction can be achieved.

There are other kinds of transforms though, which presumably could compact signals based on other features.  That is pretty much the whole point of compressed sensing.


Part 2: Root Locus, Transfer Functions and Unit Step Response? I'm continuing my question referenced here: Part 1 Question / Problem Description
Say I have a new Root Locus shown below
Consider the generic feedback loop, and the transfer function $G(s)$ shown by the following root locus plot.

Where $\mathbf{x}$ denotes the open-loop poles, $\square$ denotes the closed loop poles, and $\circ$ the open loop zeros.
I want to determine if this root-locus produces any of the following output responses to a unit step reference signal:

Question:
I don't fully understand how Matt goes from the closed-loop characteristic equation to solving for the value $K$. Is he comparing coefficients? If so, how can I do that for one real axis pole and two imaginary? I can find the closed-loop transfer function, but it results in a cubic polynomial. Also, his process determines the final settled value, but A) and B) both settle to 1, so how can I differentiate the two?
Thanks!

From the root locus plot, the open loop transfer function has a real-valued pole at $s_0=0$ and a complex conjugate pole pair $s_1=-1.5+ 2j$ and $s_1^*$. There's also a zero at $s=-2$. Consequently, $G(s)$ is given by
$$G(s)=\frac{s+2}{s(s-s_1)(s-s_1^*)}$$
The closed-loop transfer function is
$$H(s)=\frac{kG(s)}{1+kG(s)}=\frac{k(s+2)}{s(s-s_1)(s-s_1^*)+k(s+2)}\tag{2}$$
The final value of the step response equals $H(0)=1$, which leaves options $A$ and $B$.
Looking at the root-locus plot, we can see that the real-valued pole is (a bit) closer to the imaginary axis than the complex-conjugate pole pair. So that real-valued pole is (slightly) dominant. Based on this observation one can guess that step response $A$ is correct because it clearly shows the influence of the real-valued pole (i.e., a smooth increase with no overshoot), whereas step response $B$ just looks like a step response of a second-order system with a complex-conjugate pole pair.
In sum, apart from looking at the final value of the step response given by $H(0)$, you can check if there are one or more dominant poles that determine the behavior of the step response. You should be familiar with the step responses of first and second-order systems; this will allow you to make an educated guess also for higher order systems. This website gives a good overview of step responses of first, second, and higher order systems.


What is an approach in computer vision that can be used for tracking players in hockey games? It differs from human tracking in random environment in that detail that we have fixed environment. And players move pretty fast and it is important to resolve collisions and occlusions, keep track of right person. So please tell me what you think about this problem. Or what paper can be used to implement state-of-the-art approach for human tracking that you are aware of?

Regarding speed and occlusions, I would suggest you to take a look at this work on basketball: Locate and Track People in a Network of Fixed Overlapping Cameras, and notably published in A. Alahi, Y. Boursier, L. Jacques, and P. Vandergheynst. Sparsity Driven People Localization with a Heterogeneous Network of Cameras, in Journal of Mathematical Imaging and Vision, 2011. Then by looking at its citation network, you could find other relevant work.


Detect clusters in an RGB space With OCR one of the first stages in a processing pipeline is to accurately identify what is print and what is paper, generally achieved by some form of binarization.

In the good old days most scans were captured as grayscale images. For clean images with a high level of contrast, i.e. already largely black and white, Otsu's method worked quite well.

However once things like yellowing of paper, brown water stains, notes in blue biro, and other artefacts that appear as dark grey or black in a greyscale image, along with faded print are added into the mix we find ourselves in trouble.

These days almost all scanners and smart phones capture colour images and I would like to take advantage of the essentially 3D nature of RGB colour vs 2D grayscale images to improve the accuracy of the process.

To that end can you please help me with links to methods/algorithms that help to identify the centres of clusters of colour in a 3D RGB space.

At first this would appear to be some form of clustering problem and indeed, the direct answer to this question would be something like K-means clustering to split the colour space in compartments, or rather, sample clouds corresponding to different annotations on a page of text.

However, there is a simpler approach to that. Indeed, strictly speaking, grayscale is a line...or, more or less a cylinder or even an ellipsoid with its longest axis on the "grayscale" line. The line that is formed by setting R=G=B on the RGB colour cube.

Therefore, the problem of "filtering out" anything that is not grayscale ends up being the problem of estimating how far a given colour sample $u(r,g,b)$ is from the grayscale line.

You can get that information by evaluating the distance of a point to a line. Now, for a scanned document that contains both grayscale and some other colour "artifacts" (whatever these may be), the histogram of this distance is going to be bimodal. Because, the pixels that belong to the printed text are going to be close to the line (therefore, small distances) and the pixels that belong to anything that is coloured are going to be far from the line (therefore, large distances). Therefore, it might even be possible to now apply Otsu's method on the histogram of distances and find the threshold. I say "it might even be possible" because the two "classes" might end up being so much separated that Otsu's method is even an overkill.

But, that is not all. Because real world documents, photographed or scanned are almost never strictly speaking grayscale. They are likely to have a tint. The tint is simply a bias towards some colour. This would manifest itself as a different orientation of the grayscale line in the RGB colour cube. Therefore, an additional step might turn out to be that you first have to do some form of normalisation to the image to either ensure that "gray is gray" or that you can discover the sort of "tint" in the image that would imply the orientation of the "grayscale line"

Alternatively, the same task might be easier to handle in a different colour space.

Hope this helps.


How to correctly use sinc interpolation in Matlab? What is the right way to use sinc interpolation for a given discrete signal $x[n]$? Following is the sinc interpolation formula:
$$x(t) = \sum_{n=-\infty}^\infty x[n] \mathrm{sinc}\left(\frac{t-nT}{T}\right)$$
A simple explanation with a real world example of how and when to use a sinc interpolator will be very helpful.

I remember being asked to do something like this in a DSP course and it was confusing because I wasn't sure how to model the continuous variable $t$ in MATLAB where everything is digital. You can still make a nice demonstration of sinc interpolation by choosing a low sample rate and use interpolation to create the interpolated signal which is at a higher rate.
Here is a simple MATLAB script that does what I'm talking about. Hopefully it at least helps you see how to use $x[n]$ to create the interpolated signal.
% Create "continuous time" signal, Fc >> f
Fc = 1e6; % very high sample rate to simulate "continuous time"
Tc = 1/Fc; % sampling period
t = (-0.1:Tc:0.1)'; % time axis
f = 10; % signal frequency
xc = sin(2*pi*f*t); % "continuous time" signal

% Create sampled signal
Fs = 100; % sampling rate
Ts = 1/Fs; % sampling period
ratio = round(Ts/Tc);
tn = t(1:ratio:end); % sampled time axis
xn = xc(1:ratio:end); % sampled signal

% Plot the CT signal and sampled signal
figure
hold on
grid on
plot(t, xc)
stem(tn, xn, 'o')
legend('"Continuous time signal"', 'Sampled signal')

% Create and plot sinc train
sincTrain = zeros(length(t), length(xn));
nind = 1;
figure
cmap = colormap(jet(length(-floor(length(xn)/2):floor(length(xn)/2))));
ax = axes('colororder', cmap);
hold on
grid on

plot(t, xc, 'k', 'LineWidth', 3)
for n = -floor(length(xn)/2):floor(length(xn)/2)
   sincTrain(:, nind) = xn(nind)*sinc((t - n*Ts)/Ts);
   p = plot(t, sincTrain(:, nind), 'LineWidth', 2);
   stem(tn(nind), xn(nind), 'Color', p.Color, 'LineWidth', 2)
   nind = nind + 1;
end
xlabel('t')
ylabel('Amplitude')
set(gca, 'FontSize', 20, 'LineWidth', 3, 'FontWeight', 'bold')

xr = sum(sincTrain, 2); % sum(sincTrain, 2) is the interpolated/reconstructed signal, should be equal to xc

figure
hold on
grid on
plot(t, xc)
plot(t, xr)

error = mean(abs(xc - xr).^2);


Can deep neural networks achieve real-time video analysis? Recently, convolutional neural network based, deep architectures (DNN) such as AlexNet and VGGnet have been very successful in image classification challenges (e.g. ImageNet) and action recognition/video classification tasks. They surpassed hand-crafted feature-based approaches by a large margin in various cases, and thus became the new state-of-the-art approach for many computer vision problems, including but not limited to classification, object detection, and action recognition.

In the state-of-the-art DNN papers, the focus is mostly on achieving a better classification or detection accuracy, whereas the complexity of the model and computational issues stay rather in the background. This made me think whether these DNN-based approaches would be suitable in future real-time action recognition applications (e.g. video surveillance, pedestrian detection etc.) where the test time speed is of utmost importance. It is very important, for example, whether an algorithm can classify the captured video in real-time and respond very quickly to the classification decision to take some pre-determined action.

Are there any DNN-based action recognition algorithms published so far that can satisfy these criteria about real-time operation (possibly even at the expense of classification accuracy)?

I believe for many practical applications, achieving real-time or even faster-than-real-time performance is and will be much more important than providing marginally more accurate classification results compared to the previous state-of-the-art-.

Nvidia seems to have published some white papers comparing DNN inference performance between high-powered CPUs and (of course) Nvidia GPUs.  (one example)

Ballpark seems to be that some systems can meet or exceed typical video frame rate thru-put for some class of DNN classification tasks.  Whether those image sizes and/or DNN architectures and hyperparameters meets your needs is one question, and whether desktop/server-class GPUs fit your desired power envelope is another question.


Fitting new images from a SVD/PCA calculation I'm trying to replicate the ideas from the Eigenface page on wikipedia. From a hundred sample images represented by a data matrix $\bf X$ (where each image flattened to a vector of length $n$, thus $\bf X$ is a $100$ by $n$ matrix), I've computed a SVD decomposition:

\begin{equation}
\bf X = U \Sigma V^{T}
\end{equation}

hence:

\begin{equation}
\bf X X^{T} = U \Sigma^2 U^{T}
\end{equation}

By taking a subset of the largest $q$ eigenmodes, I can approximate the matrix (let $\sigma_1 \ge \sigma_2 \ge \cdots$):

\begin{equation}
{\bf X} \approx \sigma_1 u_1 v_1^{T} + \sigma_2 u_2 v_2^{T} + \cdots + \sigma_q u_q v_q^{T}
\end{equation}

Now given a new vector $y$, which represents an image not in $\bf X$, how do I determine the weighting of the $q$ eigenvectors $\bf U$ to best represent my new image $y$? Except for pathological cases, is this representation unique?

In short, what I'd like to do is this (from the wiki page):

  These eigenfaces can now be used to represent both existing and new faces: we can project a new (mean-subtracted) image on the eigenfaces and thereby record how that new face differs from the mean face.

How do I do that projection?

The "projection" that is referred to is a vector projection. To calculate the projection of vector $\mathbf{a}$ onto vector $\mathbf{b}$, you use the inner product of the two vectors:

$$
\mathbf{a_{proj}} = \langle \mathbf{a}, \mathbf{b} \rangle \mathbf{b}
$$

$\mathbf{a_{proj}}$ in this case is the vector component of $\mathbf{a}$ that lies in the same direction of $\mathbf{b}$. In Euclidean space, the inner product operator is defined as their dot product:

$$
\langle \mathbf{a}, \mathbf{b} \rangle = \mathbf{a} \cdot \mathbf{b} = \sum_{i=1}^{n}a_ib_i
$$

where $n$ is the number of components in the vectors $\mathbf{a}$ and $\mathbf{b}$ and $a_i$ and $b_i$ are the $i$-th component of vectors $\mathbf{a}$ and $\mathbf{b}$, respectively. Intuitively, by calculating the inner product of the two vectors, you find "how much of" vector $\mathbf{a}$ goes in the direction of vector $\mathbf{b}$. Note that this is a signed quantity, so a negative value would mean that the angle between the two vectors is greater than 90 degrees, as illustrated by an alternative definition for the projection operator:

$$
\mathbf{a_{proj}} = |\mathbf{a}| \cos(\theta) \mathbf{b}
$$

where $\theta$ is the angle between the two vectors.

So, given a vector $\mathbf{a}$ and a bunch of basis vectors $\mathbf{b_i}$, one can find "how much of $\mathbf{a}$" goes in each of the directions of each of the basis vectors. Typically, those basis vectors will all be mutually orthogonal. In your case, the SVD is an orthogonal decomposition, so this condition should be satisfied. So, to accomplish what you describe, you would take the matrix of eigenvectors $\mathbf{U}$ and calculate the inner product of the candidate vector $\mathbf{y}$ with each of the matrix's columns:

$$
p_i = \mathbf{y} \cdot \mathbf{u_i}
$$

The scalar value $p_i$ that you get from each inner product represents how well the vector $\mathbf{y}$ "lined up" with the $i$-th eigenvector. Since the eigenvectors are orthonormal, you could then reconstruct the original vector $\mathbf{y}$ as follows:

$$
\mathbf{y} = \sum_{i=1}^n p_i \mathbf{u_i}
$$

You asked whether this representation is unique; I'm not sure exactly what you mean, but it is not unique in the sense that a given vector $\mathbf{y}$ could be decomposed by projection onto any number of orthonormal bases. The eigenvectors contained in the matrix $\mathbf{U}$ are one such example, but you could use any number of others. For instance, calculating the discrete Fourier transform of $\mathbf{y}$ can be viewed as projecting it onto an orthonormal basis of complex exponential vectors of varying frequency.


How to calculate the time delay of a signal by a first-order filter lets assume following signal
$$x=a_1 \sin(2\pi f_1 t) + a_2 \sin(2\pi f_2 t).$$
This is passed through a simple single-pole filter with cut-off frequency $f_c$. How is $x$ delayed by the filter (peak of the cross-correlation)?

In the case $a_2=0$, the time delay between the original and the filtered frequency is
$$T_1=\arctan(f_1/f_c)/(2\pi f_1).$$

In the case $a_1=0$, the time delay between the original and the filtered frequency is
$$T_2=\arctan(f_2/f_c)/(2\pi f_2).$$

But how can we calculate the time delay $T_{12}$ in the case $a_1\neq0$ and $a_2\neq0$?

It is not the group delay: It must depend in some way on the values of $a_1$ and $a_2$, because if $a_1>>a_2$, then $T_{12}\approx T_1$.

Any light on this would be greatly appreciated!

The answer is that the output signal for the case $a_1\neq 0$ and $a_2\neq 0$ is not simply a delayed version of the input signal anymore, as it would be if the input were just a single sinusoid. With $T_1$ and $T_2$ as given in your question, the response to the combined input signal is

$$y(t)=\frac{a_1}{\sqrt{1+(f_1/f_c)^2}}\sin(2\pi f_1(t-T_1))+\frac{a_2}{\sqrt{1+(f_2/f_c)^2}}\sin(2\pi f_2(t-T_2))$$

which cannot be written as $x(t-T_x)$ with some unknown delay $T_x$.

But as you suggested, you could define some 'global delay' by the first maximum of the cross-correlation between the input signal and the output signal. This cross-correlation is

$$R_{xy}(\tau)=\frac{a_1^2}{2\sqrt{1+(f_1/f_c)^2}}\cos(2\pi f_1(\tau-T_1))+\\+
\frac{a_2^2}{2\sqrt{1+(f_2/f_c)^2}}\cos(2\pi f_2(\tau-T_2))\tag{1}$$

For $a_2=0$, there is obviously a maximum of $R_{xy}(\tau)$ at $\tau=T_1$. The same is true for $a_1=0$, which gives $\tau=T_2$ as the location of the maximum. However, for $a_1\neq 0$ and $a_2\neq 0$ I don't see how to compute the first maximum of (1) analytically. A numerical solution is of course possible.


Why is CDMA not enough for wireless communication? I have a question about CDMA, FDMA and TDMA.

I don't understand why we suffer from limited capacity when we have CDMA.

I can understand the limitation in FDMA and TDMA, as time and frequency are limited. However, we can generate as many codes as we want. When the true receiver receives the signal and multiplies by the code, one should receive the original signal even if there is a lot of interference.

So where am I wrong? Why do we suffer from limited capacity when using CDMA?

Thank, I think I made mistake by saying capacity. I specifically meant high date rate.
However, I know multiple access mechanisms.

What I asked is a different question.

Let's say I am using TDMA for 10 users. Each user has to wait 9 users to transmit. If we had 100 users, then each user has to wait 99 users. Very high level of starvation. If you want to decrease starvation, you have to decrease time slot which means decreasing date rate.

Let's say I am using FDMA for 10 users with 100 MHz bandwidth, so each user will have 10 MHz bandwidth so equivalent capacity from Shannon law. If we had 100 users which means 1 MHz bandwidth and much lower date rate equivalently lower capacity. If you want higher date rate, you need to buy more bandwidth.

Let's say I am using CDMA for 10 users, I need 10 orthogonal codes. For 100 users, 100 orthogonal codes.

So, the case is frequency and time are exist in the nature and they are limited.

However, we can generate millions of orthogonal code for millions of user as it is human-invented thing.

Why do we live problem multiple access problems while we have CDMA?

In other words, is there a mathematical or practical limit of usable number of codes ?

If the number of the users in the network exceeds a threshold, Do the codes become to similar which undermines orthogonality?

You're misunderstanding what CDMA, TDMA and FDMA do:

CDMA doesn't increase the channel capacity in any way. It's a MA, MA = Multiple Access mechanism.

In other words, it's just a way of dividing the spectrum among multiple users. No matter what you do, you can't get more data through a channel than the physics and math allows – and Shannon's channel capacity gives a pretty clear upper boundary for that.

All the mechanisms you quote hence only divide the available channel across some measure.

For TDMA and FDMA, as you notice, this is very intuitively understandable.

But to be honest, CDMA is not that different. Instead of sending a single symbol in a single symbol time slot, you have to multiply it onto a longer code sequence.

The maximum number of simultaneous users of CDMA is simply the number $N$ of orthogonal CDMA sequences.

$N$ orthogonal vectors (i.e. sequences) form a base to an $N$-dimensional vector space (which is the signal space). Vectors from an $N$-dimensional space (which is what you're going to transmit to transport one symbol) have inherently $N$ elements.

So, to make (up to) $N$ people talk at once, you need to "blow up" everyone's transmission's symbols (which are $1$-dimensional, each) to $N$-dimensional vectors, increasing the time it takes to send one symbol by a factor of $N$.

That's the same efficiency as using TDMA to split time into $N$ slots (everyone sends $1$ symbol every $N$ slots), or dividing the bandwidth into $N$ channels and using FDMA (everyone can only use $\frac1N$ of the available bandwidth, reducing the individual symbol rate by $N$).


Image Gradient - Flipping Effect of the Convolution Matrix (Kernel) for Edge Detection Prewitt operator is a common operator to compute image gradient in edge detection, and I knew that to apply this Prewitt operator to an image is just to convolve them together, say, the image is represent as $A$, and as the wikipedia page lays out, the Prewitt operator for *x-gradient is $$G_x=\pmatrix{ -1 & 0 & 1\\\ -1 & 0 & 1\\\ -1 & 0 & 1}$$, for y-gradient is $$G_y=\pmatrix{1 & 1 & 1\\\ 0 & 0 & 0\\\ -1 & -1 & -1}$$
So if I want to compute the x-gradient of image $A$, convolve them like this,$$A_x=A*G_x$$, but here comes my problem. $G_x$ is the kernel, when convolve it with $A$, $G_x$ should be flipped which turn it into this, $$G_x=\pmatrix{1 &0 &-1\\\ 1 &0 &-1\\\ 1 &0 &-1}$$
right? If so, the x-gradient will be something like this, $$\nabla A_x=A(x-1,y)-A(x+1,y)$$
, BUT I think it should be $$\nabla A_x=A(x+1,y) - A(x-1,y)$$.
So I presume that either $G_x$ and $G_y$ should not be flipped when convolve them with $A$, or $G_x$ doesn't look like the above, and it should be this,$$G_x=\pmatrix{1& 0& -1\\\ 1& 0& -1\\\ 1& 0& -1}$$
Anyone can help me with this?

This is a common question and there are a few other answers on this site which may be helpful. What you have noticed is correct. To be clear, convolution is defined by this flipping, it gives the convolution operation some nice maths properties. Non-flipping convolution is called cross-correlation.

For these kernels, the only thing flipping changes is the sign of the result. You can safely use the inverse (flip) of the kernels to get the form you desire.


Uncovering data from noise I've recently started experimenting with DSP, specifically, understanding digital signals that are generated by some small applications I am creating. For this quest, I've decided to use The Scientist and Engineer's Guide to Digital Signal Processing. While the provided material is beyond exceptional, the more I read the more confused I become as to how I could apply different forms of DSP in my "computer based" experiments.

Here is what I am currently working on:

I have a capacitative sensor. When the sensor is stretched, the capacitance increases. Contrary is true when it contracts. Image below displays this characteristic:

An interesting observation I noticed is that when I stretch the sensor, while bending it, my signal is buried in noise. This is evident between 800 and 1000 samples in the above photo. Because I am experiencing this issue, I want to uncover the stretch/contraction of the sensor from this noise; hence my goal of studying DSP.

An interesting comment made by the author, is that when ever we deal with Digital Data we will always want to analyze our Time Domains using DFT. The most efficient algorithm to do this today is FFT. I've also concluded ( high level ) that using the mentioned algorithm is ideal as it will give me a sense of what "sin waves" generate the "digital wave" I am presented with. Hence, I should be able to parse out "sin waves" of interest to re-construct missing points of interest hidden in noise.

Now, while I understand that the change to the frequency domain can tell me that "a spike has occurred", what I am really interested in is re-constructing my wave to be less noisy. Hence, I'd like to derive out a consistent wave form between 800 - 1000, alike I have between 400 - 600.

Questions

Since my understanding is still fairly low, is DFT and FFT really the solution to my problem?
Would it be logical to take the FFT , filter out outliers that are generated and cause noise, and re-construct the time domain signal back ( with the results that were acceptable ) using iFFT? That should, in theory, give me back a time domain with less noise, no?
Is there a better section of DSP that I should read that can help me analyze the decouple noise from interesting data in a time domain?

Since my understanding is still fairly low, is DFT and FFT really the solution to my problem?

All the FFT (or DFT) really gives you is a way of visualizing the data in the frequency domain.

Your data may be good for this, except that there will be a large spike at "DC" (the zero frequency bin) because your data has a non-zero mean.

  Would it be logical to take the FFT, filter out outliers that are generated and cause noise, and re-construct the time domain signal back ( with the results that were acceptable ) using iFFT? That should, in theory, give me back a time domain with less noise, no?

If you know the approximate frequency of the variation you're after, yes, you can use the FFT (and IFFT) to filter out frequencies different from the band where the frequency will be.

One simple thing you could try without using an FFT is to use a DC blocker.  That will make your signal approximately zero mean, and will also smooth out some of the roughness you see.

  Is there a better section of DSP that I should read that can help me analyze the decouple noise from interesting data in a time domain?

One place to start might be chapter 15 of that book.

I suspect you'll want to design a bandpass filter around the frequencies of interest and filter out everything except this band. You will still see some effect of the noise, but it should be greatly reduced.

If that isn't good enough, then you may need to try more advanced techniques which don't appear to be covered by the book.


how to handle different durations of audio data? I am new to signal preprocessing, I read about mel_spectrograms, MFCC's. Now I want to apply it and use the CNN model, But the data which I have for practice is having audio of different durations, now because of this, the mel_spectrograms will be of different shapes. for using them as inputs, the model requires inputs to have a fixed shape. So, what should I do to make them have a specific shape??

In this case, usually, Normalization is done.
For example in your training and testing data, you need the same shape, so that you should try something like,
mean = np.mean(X_train_features, axis=0)
std = np.std(X_train_features, axis=0)
X_train_features = (X_train_features - mean)/std

Here X_train_features can be a data frame of spectrogram or mfcc features.
The same thing you can apply for testing features. One thing to remember is the number of columns in shape of training and testing should be the same.
You can also check some Kaggle kernels which are related to Audio processing.


FFT and frequency response There is a system (consisting of several LTI analogue filters) whose input and output signals are being recorded at $20 \ \mathrm{kHz}$ (the sampling period is in the range of 10 times the system's transient response). Suppose I calculate

$$ G = \frac{\mathrm{FFT(output)}}{\mathrm{FFT(input)}}$$

What does that correspond to? Can you call it the frequency response? If it does not have anything to do with the frequency response, can someone provide an (ideally mathematical) explanation why it doesn't? I think, that it is legit if there were no transients present in the measured output. But since they impact the signal significantly I wonder if one can still get useful information from $G$.

EDIT: I guess I should rephrase the question. I know that the frequency response is defined as $$X(\omega) = \frac{Y(\omega)}{U(\omega)}$$. The issue here is that $y(t)$ and hence $Y(\omega)$ are greatly impacted by the system's transient response. May $X(\omega)$ still be called frequency response ?

Similar questions have been asked before, you may search. (I'm sorry but I cannot refer to one here)

You want to identify an analog LTI system from analysing its response $y(t)$ to a known input $x(t)$.

One method of analysis involves the use of continous Fourier transforms of the input and output signals. It can be shown that (for a stable system)
$$H(j\Omega) = \frac{Y(j\Omega)}{X(j\Omega)} $$

Considering a continuous time analog LTI system, this method requires some convenient techniques to compute the Fourier transforms of the input and output signals. This method has practical limitations due to understandable reasons.

With the advent of digital computers and availability of dedicated hardwares, a more implementable approach has been developed during the rise of digital signal processing era, based on the sampling theorem to represent the continuous time signals by their samples, i.e., the continuous time analog system is represented in the form of an equivalent (and in a one to one relation to its continuous time counterpart) digital LTI system, as long as some necessary conditions are met.

An analysis based on the taken samples of input and output signals, however, should be handled with care, as the most critical element of success in this method relies on the fact that the continuous time system's frequency response function $H(j\Omega)$ be bandlimited to the half of the sampling frequency $\Omega_s = 2\pi F_s$, so that it can be digitally represented in an equivalent digital LTI filter with DTFT $H(e^{j\omega})$

If you know that this condition is satisfied, then there are other details such as the fact that Fourier transforms provide information on steady-state behaviour and not on the dynamic transients that will happen in any finite length measurements. (That can be avoided with longer recording durations, as you intent to do)

Also note that the LTI analog system must be stable so that its CTFT $H(e^{j\Omega})$ exists. Most often it will be stable I assume, but it's still important to consider though.

Finaly any noise and/or distortions on both the input and output signals not only prevent high precision computations but can even lead to misleading results.

So I suggest you to begin with a known setup first and apply your method to see if it can also identify that known LTI system. If yes, then you are ok.


What factors determine synchronous signal transmission ability of a channel? So I am Designing a system which would be able to send and receive stereo voice signals.( like a telephone but with full-band stereo voice signals)

here is an example to make my question clearer.
In a music recording studio two microphones are foreseen, one for singer's voice and the other for the guitar's sound. The combination of these two voice signals will result into a Stereo signal.At the same time there is a live conversation with another end which is supposed to receive music studio's signals as stereo and be also able to send back stereo signals.
If Bluetooth were to be used to transmit voice signals to CE Device( smartphone) and VoIP to transmit between two sides of conversation . I have done some studies but apparently Bluetooth has some limitations in case of synchronous transmissions.
what factors should be taken into account to choose a transmission method?

Not sure what exactly you mean by "synchronous". Most digital transport mechanism (such as Ethernet or 802.11) are "asynchronous", i.e. a continuous data streams  are chopped up in packets, the packets are transmitted whenever the connection allows it, the packets are assembled in the correct order at the receiving end and then turned into a continuous data stream again. Since it's not known when exactly the packet can be transmitted, there has to some buffering at the receiving end so there is some data available even if the next packet comes in late. This results in a certain amount of "latency". In the case of Bluetooth + VoIP you are probably looking at several 100 milliseconds to a few seconds of overall latency.

You may also be concerned about inter-channel sync, i.e. do both channels of the stereo stream stay time aligned. This can easily be achieved by interlacing the channels to create a single stream.

This all depends heavily on the specific implementation: Bluetooth supports many profiles. For your application you probably want the A2DP profile. This uses the SBC codec and there are dozens of ways to configure the codec ranging anywhere from high quality to totally unusable. Both parts of the interface need to agree on the codec implementation.

There are also many different VoiP protocols and codecs out there. Most operate in the 10kb/s-40kb/s range which seems way too low for you application, so you may need to create a custom protocol.


Recommendations for learning and getting hands-on Real Time Signal processing I know this is not a technical question, but wanted to get some opinion on how to get some hands-on experience on embedded DSP / real-time signal processing
I am comfortable with DSP and wireless communications simulations in MATLAB /python and wanted to get into embedded DSP development and I'm not able to get any structured way of learning real-time DSP (I regret not taking RTSP class in my masters) . I'm proficient in C and have a very basic level knowledge in  assembly. I'm ok to buy the hardware for my learning curve.
I tried searching for any online training /course from universities which offers such class where we start writing some vector DSP code on the embedded platform (TI or ARM etc).
Any recommendations/insight into this would be highly grateful.
Thanks
Saira

I had taken the course at my university. So most of the materials were provided as part of it. I worked with 2 different DSPs both of which were provided at the lab.
If you are looking to buy hardware, Real-Time Digital Signal Processing by Sen M Kuo and Bob H Lee is a good book with implementation and architecture descriptions of TMS320C55x. It has around 9 interesting application-level implementations as well. But getting all the hardware necessary might be too much though.


Average image over $N$ frames with one frame buffer Can you give some filter that is averaging an image over $N$ samples in a recursive way that only needs one frame buffer?

$\frac{1}{N} \sum_{k=1}^{N} y_{k}(m,n)$

I can only imagine filters that need $N-1$ frame buffers.

Don't fully understand the question; your formula already describes your solution:

Take the first frame into your framebuffer
For every pixel of the second frame, add it to pixel at the same position in the frame buffer.
Repeat 2. for the rest of your frames.
Multiply every pixel of the resulting frame by $\frac1N$.

Done!

(Note: depending on the data type, might be advantageous to multiply every pixel with $\frac1N$ before addition/move to the framebuffer pixel)


Impulse response of a causal system from transfer function in z-domain The transfer function is $$H(z)=(z+1)/(z^2+z+0.5)$$
I need to find the impulse response h[n] of a causal system with x[n] as unit impulse.

I have tried to find the impulse response by the following methods:

partial fraction decomposition and then inverse z transform of transfer function.
homogeneous solution of differential equation.
Matlab iztransf() function
I get different solutions from the above methods and i don't know how to remove iota from the impulse response.What is the correct procedure to find the impulse response?

Let's follow option 1. I'm using Python but same is directly available in matlab. The partial fraction expansion can be found via residue function

>>> import scipy.signal as sig
>>> b = [1, 1]
>>> a = [1, 1, 0.5]
>>> r, p, k = sig.residue(b,a)
>>> r
array([0.5-0.5j, 0.5+0.5j])

>>> p
array([-0.5+0.5j, -0.5-0.5j])

>>> k
array([0.])

This gives us
$$
H(z) = \frac{0.5-0.5i}{z+0.5-0.5i} + \frac{0.5+0.5i}{z+0.5+0.5i}
$$

Now either by computing manually, or using a Z-transform table, or using another computational tool you can show that these are two complex terms can be inversed transformed as

$$
\mathcal{Z}^{-1}\left[\frac{a}{z+a}\right] = \mathcal{Z}^{-1}\left[\frac{az^{-1}}{1+az^{-1}}\right] = -(-a)^k \mathbf{1}(k - 1)
$$
Here $\mathbf{1}(k)$ is the discrete time step function. But since we have complex numbers this looks like our response should be a complex function. However this is not possible since our transfer function has real valued coefficients. Hence even before I compute anything, I know already that imaginary parts will cancel out.

And since $|a| < 1$ this is going to be a decaying sinusoidal. Indeed if we check this with say sampling period $T=1$,

>>> import harold as har
>>> G = Transfer([1, 1], [1, 1, 0.5], dt=1)
>>> impulse_response_plot(G);

As expected, there is a one sample delay coming from $\mathbf{1}(k-1)$


Make a signal that fits another the best possible with a limitation in the 2nd derivative Consider this step function:

The signal that "fits" this should look like the following (in green):

The corners are now smooth because the maximum second derivative allowed is not infinite anymore.

The idea is to find the most close to the original signal given the maximum second derivative allowed, minimizing the difference between the two signals (in orange):

If your system can predict things before it happens it could do something like this:

Minimizing the total difference even more. But also could be achieved if it works in non real time or has some kind of look ahead. But that will be not always possible or desirable.
This is not in the very core of my question, however is also a interesting option.

So, my question is: How can this be done? Is it possible (eg: exist more than one curve that minimizes error to the same value)?

Also to me looks like a filter, is there any filter which definition is something like this (hopefully a digital one)?

Hmmmmmmmmm, interesting question.

Since you want to use the second derivative as your criteria, it would seem that you would want to have the maximum second derivative absolutie value for as short of a duration as possible.  This would suggest piecing together parabolas, matching the first derivatives at the joints.

How to do this algorithmically will take a little more thought.  Somehow you have to identify regions of "violations" and then figure out how big of an interval you have to use to patch it piecewise with parabolas.

Just some thoughts.  If I think of a sure fire way, I'll add it on.

Ced

Followup:

It turns out to be really simple to make a FIR that reacts to a step function with a upward parabola for the first half and a downward parabola for the second half.  It just becomes a matter of selecting the width and the rise rate to set the value of the second derivative.  I haven't done the math yet for the actual value, but here is a program which gives proof of concept.

 h[k] = [ 1  2  3  4  3  2  1 ]  / 16

It appears the second derivative proxy $ y[n-1] - 2 * y[n] + y[n+1] $ has an absolute value of the normalizing factor of the integer step vector, e.g. $1/16$.

import numpy as np

#======================================================================
def main():

        P = 2

#---- Build the Weights

        K = 2 * P - 1

        h = np.zeros( K )

        for k in range( P ):
            h[k] = k + 1

        for k in range( P - 1 ):
            h[P+k] = P - k - 1

        theRecip = 1.0 / sum( h )

        print h, sum( h )

        h *= theRecip

        print h

#---- Set up the Step Function

        x = np.zeros( 20 )

        x[10:20] = 1.0

        print x

#---- Apply Filter

        y = np.zeros( 20 )

        for n in range( 8, 20 ):
            Sum = 0.0
            for k in range( K ):
                Sum += x[n-k] * h[k]
            y[n] = Sum

#---- Show Result

        for n in range( 6, 18 ):
            second = y[n-1] - 2.0 * y[n] + y[n+1]
            print n, y[n], second

#======================================================================
main()

Yet a better answer:

[Moved to a separate answer]


Optimized Ansi C libraries for DSP I am new to DSP? Where is a good resource for open source DSP algorithms?

MATLAB is great at making protos but once we move to C coding, it takes time and we end up making too many mistakes. I would guess there must be good repositories of C codes for various implementation of filters etc. Any pointers would be good.

The rub is that ANSI C is not a terribly good language for real time DSP. Many features that a real DSP chip has can't be addressed in ANSI C. These include address generators, bit-reverse addressing, circular addressing, SIMD, memory segments, multiple buses, fractional data types, fractional multiply, certain fixed point math, zero overhead looping etc. Hence most compilers for DSPs have an extensive set of language extensions (#pragma, inline assembler, new keywords, etc). However these are decidedly non-ANSI and not portable.

Writing a simple fixed point FIR filter in ANSI C can easily be 10 times less efficient than using the proper extensions or assembly code.


Pole-zeros of a real-valued causal FIR system Could someone please help me with the following question?

Below is the magnitude response of a real-valued causal linear phase FIR system of order N = 6. Determine the location of poles and zeros.

I know that for FIR systems all the poles are located at the origin, so we have a pole of order six at the origin. Also from the given diagram, I can say that we have a zero at 0.3pi and one at 0.8pi (both on the unit circle). Now since the system is real-valued, location of poles and zeros should be symmetric w.r.t. the real axis. But I don't know about the two other zeros?

Also, what about the pick in the diagram? Does it mean we have another pole?

Note the difference between the zeros at $0.3 \pi$ and at $0.8 \pi$.

The first one is clearly a zero crossing, much like $abs(x)$ at $x=0$.

At $\theta = 0.8 \pi$, however, the curve is tangent to the horizontal axis, much like $x^2$ at $x=0$. So you have a doulbe zero here.

So your zeros are:

2 zeros at $z = e^{\pm j 0.3 \pi}$
2 double zeros at $z = e^{\pm j 0.8 \pi}$


Projection of 2D signal Paper describes that a 2D signal(say video) can be projected to 1D for time series analysis by using a transformation function fx=Ax+constant and fy=Ay+constant where Ax , Ay represents the (x,y) coordinates of a single trajectory A. Then,these fx vs t can be plotted and various analysis can be performed. Going by their justification,is it logical to think that a coordinate position of an object can only be described by a single coordinate when doing time series analysis and signal compression?

Also,they mention about taking cartesian product.Why is that?

Yes you can re-represent the 2d signal as a 1-d signal and use signal processing algorithms on it. In doing this you can lose some of benefits of representing the data this way.

Eg. If you were looking at compression of a mostly continuous 2-d dataset you may get benefits by treating the data as 2-dimensional and compressing in each dimension. Representing the data as concatenated rows or columns introduces discontinuities that the compression won't handle as well.


MIDI out from freeware sequencer to PD or Max/MSP I'm looking for a way to send MIDI from a cheap or free MIDI sequencer to Pure Data or Max MSP so I can do realtime MIDI processing in PD/MAX. So I guess I have a couple questions:

What is a good cheap/free midi sequencer (that runs in Windows 7) that supports MIDI out to another program
How do I configure MIDI out?

Is this something normal midi sequencers can do? It seems like it should be pretty basic, but I was having trouble finding information on it.

Max/MSP has the midiparse object which is incredibly useful for handling midi data. Midi Ox is a program that will help you with routing occasionally. As far as freeware, you get what you pay for. Hydrogen is an open source a drum machine that will send midi out. You'll just to configure Max or Pd to receive the channel you're sending out of your sequencer.


Is a recording of a chord merely a superposition of recordings of its constituent notes? Is a recording of a chord merely a superposition of recordings of its constituent notes? Assuming the recording of its constituent notes are sufficiently clean.

It's close enough for most practical synthesis purposes.  However in real polyphonic instruments, there tends to be some coupling between the various oscillators that can, for example, synchronize notes that are close enough to whole number ratios.

For other ensemble instruments, the coupling can work via the hearing of the instrumentalists and minute corrections (of embouchure, fingering pitch, and so on).

And for nonlinear effects (like accordion bass mechanic sympathetic vibration or a distortion pedal for a guitar or a distorted amp), the results are definitely non-linear.

But more often than not, just adding the samples is a good enough approximation.


Thanks to  $u(t)$,
 $$
\begin{align}
   \int\limits_{-\infty }^{+\infty }{{{e}^{-t}}\sin(2\pi ft)u(t){{e}^{-j2\pi ft}}\,dt}=\int\limits_{0}^{+\infty }{{{e}^{-t}}\sin(2\pi {{f}_{c}}t){{e}^{-j2\pi ft}}\,dt}
\end{align}
$$
Now
$$0 \le | e^{-t}\sin(2\pi {{f}_{c}}t){{e}^{-j2\pi ft}}| \le | e^{-t}|$$
so your function is absolutely integrable.

The trouble comes for the exponential product rule:

$${{{e}^{-t}}{{e}^{j2\pi {\pm{f}_{c}}t}}{{e}^{-j2\pi ft}}} = e^{-t\left(1\mp j2\pi {{f}_{c}}+j2\pi f\right)}$$
which integrates quite well.


number of samples and Amplitude accuracy,what's the solution to this problem I'm colleceting data samples from a device (temporal signal) in order to analyze it with FFT function in MATLAB R2015b,first i get N=4063 data points with (delta_t=0.000015s) so my frequency resolution is (delta_f=1/N*delta_t)=16.4Hz,this give me pretty good results with amplitude,but when i did another measurment with N=8192 data points (delta_t=0.000061)  to obtain better resolution,the amplitude is not correct,i don't know really how to deal with this issue.so please if any body can help me with this i will be grateful,thank you

ps:i'm new in signal processing

below my code in MATLAB

sig_2=asp_ach_de1;
Fs2 = 1/0.000061;      % Sampling frequency
T2 = 1/Fs2;            % Sampling period
L2 = 8192;             % Length of signal
t2 = (0:L2-1)*T2;      % Time vector
% FFT FUNCTION
Y2 = fft(sig_2(:,2))
Q2 = abs(Y2/L2);
Q1 = Q2(1:L2/2+1);
Q1(2:end-1) = 2*Q1(2:end-1);
f2 = Fs2*(0:(L2/2))/L2; %FREQUENCY
figure(3)
plot(f2,Q1)

I think that what is happening to u is that the frequency which you are trying to measure is in between 2 of your FFT bins which will cause the power to spread between two different indexes of your result FFT.

Another possible solution , only if with higher resolution you are getting a lower value  , is that since as I understand this is a real signal with noise you may be summing more noise , on the first attempt, for each frequency of your FFT and as a result raising all the values of your FFT.

Without knowing the signal you are trying to process it is hard to know what is exactly the problem. also your FFT plots will help give a better answer.


What are the 'phases' in this 2D scaler filter I am pretty new to DSP and have been looking at how image scaling is done in hardware. I and came across a 2D scaler core for an FPGA here. I have read up a bit about polyphase filters, but am confused about what the "16 phases" referred to on page 3 with the following diagram are.

I have not come across this elsewhere and was hoping for some insight as to what is going on here.
My best guess is that they limit the available interpolation points to 16 discreet points in $x$ and the same for $y$, and so any fractional scaling where a sample point falls between $2$ of these $1/16^{\rm th}$ sampling points, will just use the closest point to it?
Any insight would be much appreciated.

My best guess is that they limit the available interpolation points to 16 discreet points in x and the same for y, and so any fractional scaling where a sample point falls between 2 of these 1/16th sampling points, will just use the closest point to it?

Yes, that's it.  Probably on the assumption that within 1/16 of a pixel shift, and with the further interpolation on the output, that filter will be good enough.