Oscillations in a first-order discrete-time linear systems According to my understanding, continuous first order systems do not exhibit oscillatory behavior because a first order system has only one energy storage element and thus oscillation is not possible. I was simulating the below difference equation for a first order discrete system:

$$y[n] =  \tfrac{1}{10} u[n] - \tfrac{7}{10} y[n-1]$$

where $u[n]$ is a unit step function. As can be seen, this equation describes a first order difference equation. However, when plotting $y[n]$, I can see that it oscillates at the start and settles next:

However, when I simulate the following difference equation then this oscillation disappears:

$$y[n] =  \tfrac{1}{10} u[n] + \tfrac{7}{10} y[n-1]$$

I am not sure how first order discrete time systems are different from their continuous counterparts, and why would the sign change the stability of the system. Any explanation?

For continuous-time systems, a pole at location $s_0=\sigma_0+j\omega_0$ will create a time-domain contribution of the form

$$e^{s_0t}=e^{\sigma_0t}e^{j\omega_0t}\tag{1}$$

which is a damped oscillation if the pole is in the left half-plane (i.e., $\sigma_0<0$), and if the pole is not on the real axis (i.e., $\omega_0\neq 0)$. For $\omega_0=0$ there is no oscillation. This is the case for a real-valued first-order system, which can only have a real-valued pole.

For discrete-time systems, a pole at $z_0=re^{j\omega_0}$ ($r\ge 0$) causes a sequence of the form

$$z_0^n=r^ne^{j\omega_0n}\tag{2}$$

If the pole is inside the unit circle (i.e., $r<1$) that sequence will decay in amplitude. Note that the sequence will not oscillate only if $\omega_0=0$, i.e., if the pole is on the positive real axis. If the pole is real, but if it is negative, we have $\omega_0=\pi+2\pi k$, i.e., we have an oscillation with maximum frequency, corresponding to an alternating sequence.

In sum, for continuous-time systems if the pole is on the real axis, it will not cause oscillations. For discrete-time systems, there are no oscillations only if the pole is on the positive real axis. A pole on the negative real axis causes oscillations with maximum frequency. Consequently, a discrete-time real-valued first-order system can exhibit oscillations if the pole happens to be on the negative real axis.

For additional insight consider the exact mapping between the $s$-plane and the $z$-plane:

$$z_0=e^{s_0T}\tag{3}$$

where $T$ is the sampling interval. Inversion of $(3)$ gives

$$s_0=\frac{\ln(z_0)}{T}\tag{4}$$

A $z$-plane pole $z_0=-r$ ($r>0$) on the negative real axis maps to a complex valued pole in the $s$-plane:

$$s_0=\frac{\ln(-r)}{T}=\frac{\ln(r)}{T}\pm j\frac{\pi}{T}\tag{5}$$

The frequency $\omega_0=\pi/T$ corresponds to half the sampling frequency, i.e., the maximum frequency of the corresponding discrete-time system.

As a final note, oscillations have nothing to do with whether the system is stable or not (stability in the bounded-input bounded-output sense). For stability, it only matters if the contributions of the poles to the output signal decay or not.


can MFCC coefficients alone be used to do voice detection? Aim of my project is to differentiate between cases of silence, single active speaker and multiple ( simultaneous ) active speakers using as basic ( but effective ) logic as possible. Is it possible to achieve this without use of extensive machine learning algorithms ??

Actually, i don't need to do speaker "identification". An audio stream will be taken as an input and I have to just say whether that audio consists of a single speaker, multiple speakers or silence... I have to implement this in C (real time), hence machine learning (unsupervised) would increase the computational complexity.

In essence, yes.

You would have to use something like Dynamic time warping to measure the similarities between two given signals, as long as they vary in time or speed.

From what I can remember, you can use DTW for speaking identification. BUT it is important to note that for single word recognition and single speaker identification this method would work but would not be suitable for commercial or large scale applications.

The DTW algorithm calculates the distances from the two signals, based on the Euclidean distances.. Whilst studying for my undergraduate degree, I did a simple word recognition system (In C++) and the machine learning algorithms (for instance, Hidden Markov Models) provided a much higher success rate than DTW.

Feel free to take a look at some of the code that I wrote for this project. Here

double Euclidean_distance(const vector<double> &actual, const vector<double> &training) {

    double distance = 0.0;

    for(unsigned i=0; (i < actual.size()); i++)
    {
        distance = pow((actual[i] - training[i]), 2);
    }

    return sqrt(distance);
}

double Distance(const std::vector< std::vector<double> > &actual, const std::vector< std::vector<double> > &training)
{
    int m = actual.size();
    int n = training.size();

    double cost[m][n];

    cost[0][0] = Euclidean_distance(actual[0], training[0]);

    for(int i = 1; i < m; i++)
            cost[i][0] = cost[i-1][0] + Euclidean_distance(actual[i], training[0]);
    for(int j = 1; j < n; j++)
            cost[0][j] = cost[0][j-1] + Euclidean::Euclidean_distance(actual[0], training[j]);

    for(int i = 1; i < m; i++)
            for(int j = 1; j < n; j++)
                cost[i][j] = std::min(cost[i-1][j], std::min(cost[i][j-1], cost[i-1][j-1]))
                    + Euclidean_distance(actual[i],training[j]);

    return cost[m-1][n-1];

}

(Note this is Academic code, not tested and just illustrate the algorithm)


Which one is the correct graph for $u[2-n]$? $u[n]$ is the unit step function in discrete form. I want to know which one of the following two is the correct waveshape for $u[2-n]$, where $k$ is a constant.

What's the answer? The top one is $u[n]$, the unit step function. So, what's the graph for $u[2-n]$, middle or bottom?

HINT: For which value of $n$ does the argument of $u[2-n]$ become zero? That's where the step occurs. For which values of $n$ is the argument non-negative? That's where your unit step equals $1$. If you think about it for a minute, it should become really easy.


Efficient way to get power levels from an SDR (RTL-SDR specifically) At the moment I can step through the frequency band of interest and take samples to estimate relative power at various frequencies.  But this is slow.  Is it effective to center myself in the band of interest, make sure my sample rate is high, and do an FFT instead?  I'm just wanting to efficiently get relative power levels for frequencies throughout a band.

As far as I understand your question, you want to estimate the PSD in a large bandwidth.

What you're currently doing is tuning to one part of that bandwidth, receive samples at a low sampling rate, calculate the power in that small bandwidth, retune, and repeat.

So, yes, using a larger sampling rate will obviously reduce the time you need to sense a large bandwidth: A larger sampling rate will give you more spectral information.

The DFT->square method is one classical spectrum estimator. So, yes, instead of tuning $N$ times to receive an $N$th of the overall bandwidth $B$ by sampling at a rate of $\frac BN$, you could just as well sample with a rate of $B$, and calculate $N$-point DFTs in the same time. Since sampling the same time at rate $B$ gives you $N$ times the samples of sampling at $\frac BN$, the information would be the same (Nyquist/Shannon).

You might want to look into other spectral estimators, too – but for now, I don't really see a big disadvantage in simply using the DFT with an appropriate Window – and given the fact that neither the analog nor the digital filters on the RTL dongle probably are overly fantastic, I think you can just stick with a Hann window and be fine.


Slice selection in MRI Looking at some examples in the "handbook of MRI Pulse Sequences" by Bernstein, I see the following:

So if I can generate an RF pulse with a bandwidth of 2 kHz and I want a 3 mm slice thickness, the gradient amplitude can be chosen with:

$$
G_z = \frac{2\pi \Delta f}{\gamma \Delta_z}
$$

Plugging the values, this gives a gradient strength of 15.6 mT/m.

Now, my question is that will this gradient select a 3 mm thick slice at the magnet isocenter?

Now say in my MR experiment I am acquiring 128 such slices. Now, if I want to excite the "bottom most slice" than I need to adjust the RF bandwidth and the change in RF BW is given by:

$$
\delta f = \frac{\gamma G_z \delta_z}{2 \pi}
$$

Here now, my $\delta_z = -64 * 0.003 m$. Hence, the new RF BW will need to be

$$
f = \Delta f - 64 \Delta f = -126 khz
$$

I am not sure if this calculation is correct. The RF value seems quite high to be although I am not sure what is the range of modern MR scanner RF generators.

The bandwidth of the RF pulse does not change at all - it is the center-frequency of the RF-pulse that needs to be adjusted. The slice thickness is related to the frequency spread caused by your slice-selection gradient (as you see in the first formula). Hence, to achieve the same slice thickness for each slice, the RF bandwith may not change, but the RF center frequency.


modified Bessel function for Speech enhancement Using MMSE in matlab I am trying to implement an algorithm proposed in an old research paper by Ephraim and Malah for my DSP course project. I hope someone can help me with this difficult question. For the clean signal Amplitude Estimator, they used the following equation:

I understand most of the variables in the equation except the Bessel function stuff. Why did they use a Bessel function in here?
and most importantly, how can I implement it in Matlab?

The Bessel function just turns up as a result of the MMSE estimation problem given the assumed statistical model (cf. Appendix A of the paper). This modified Bessel function is implemented in Matlab as besseli.


given an image and its fourier transform, find another image fft using fourier properties I managed to solve few of those questions , but I'm not manage to answer 3 of them correctly,

the question is: given an image and its fourier transform, find another image fourier transform ,using fourier properties to explain why it would look like that.

Imagine the same problem in 1-dimension and you will be able to find the answer.

Hint: If you know the Fourier transform of a delta function (or time-shifted delta function) and the Fourier transform of a rectangular window, then you should be able to answer the above problem.

P.S. Cannot comment because I need 50+ reputation to do so.


Convolution and Cross Correlation Blackbox Suppose I was handed a blackbox that I know is either a Convolution filter or Cross-Correlation filter, and my job is to find out which it is.
My idea is to pretend that it's Cross-Correlation and attempt to do Template Matching using an image and a crop of the same image. If the output has a peak where a successful template match would be, it's Cross-Correlation. Otherwise, it's Convolution.
Is this a reliable test? Is it possible to also arrive at the same expected output if the blackbox was actually Convolution?


Connection between SNR and the dynamic range of the human ear For signal-to-noise ratio (SNR), SNR=6B dB  is called the dynamic range of
the quantizer.

The dynamic range of the 16-bit quantizer is 6B = 6·16 = 96 dB. Note that the dynamic range of
the human ear is about 100 dB. Therefore, the quantization noise from 16-bit quantizers
is about at the threshold of hearing. This is the reason why “CD quality” digital audio
requires at least 16-bit quantization.

How can we obtain(calculate) the fact that  the dynamic range of
the human ear is about 100 dB?

Where does this quote come from? It does NOT look like a reliable source of information to me.

The dynamic range of the 16-bit quantizer is 6B = 6·16 = 96 dB.

Wrong. The quantization noise for a 16 bit bipolar rounding quantizer is about -101dB dFS. That makes the signal to noise ratio for a sine wave 98dB and for the "average" music file about 8dB-90dB.

Note that the dynamic range of the human ear is about 100 dB.

Grossly oversimplified and misleading. The dynamic range of human auditory systems depends highly on frequency (and unfortunately also and age and how much rock& roll you have played, sigh). Max dynamic range is about 120 dB at the midrange.

This is the reason why “CD quality” digital audio requires at least 16-bit quantization.

16-bit quantization was chosen as a complicated trade off between audibility of quantization noise and cost & technical complexity at the time the standard was created. 16-bit per se is in many cases not sufficient unless proper "dithering" and "noise shaping" is applied as well.
The explanation really misses the main point. SNR and audibility of quantization noise depends A LOT on frequency and the specific signal. 16-bit PCM quantization is a "brute force" approach with most bits being wasted and a few important ones missing. Perceptual codecs such as AAC, MP3, Ogg, etc. do this much better and generate way more "sound quality per bit".
Finally, the question:

How can we obtain(calculate) the fact that the dynamic range of the human ear is about 100 dB?

Through listening test performed over a wide range of diverse subjects and test signals. The hearing threshold is determined as the sound pressure where the subject says "I can't hear this anymore" and the pain threshold corresponds to "ow, that hurts. Stop it!". The pain threshold is particularly difficult to test since it's fairly subjective and the testing can cause physical damage to the test subject. The results of the hearing test are then compiled using statistical analysis. You end up with a set of curves like this: https://en.wikipedia.org/wiki/Equal-loudness_contour .


Remove noise on nonlinear signal This signal represents a drawn line. For some reason, the device showed interference as shown in the figure below:

The question is: how can I remove the interference, making the signal as continuous as possible?

Any suggestions?

It looks like you are plotting a histogram, which is basically an estimate of a probability mass or probability density function. Maybe this is not correct. If your x-axis label is wrong, and this is truly a signal as a function of time, then what you are looking for is an estimation of the actual signal given additive noise. A simpler solution is a frequency-based low-pass filter. People can provide more answers if they know the nature of the signal in your graph.


Which is the total energy of the product of two discrete energy signals Assume that I have two signals with finite energy.

The first, $x_s$
\begin{equation}
E_s = \sum_{i=0}^{N} |x_s(i)|^2
\end{equation}

and the second, $w$
\begin{equation}
E_w = \sum_{i=0}^{N} |w(i)|^2.
\end{equation}

Then, which is the enrgy of the follwoing signal that combine the two ones as follows

\begin{equation}
E_s'=\sum_{i=0}^N |x_s(i) \cdot (1+K w(i))|^2
\end{equation}

In particular assuming $x_s$ and $w$ uncorrelated, is there a way to express $E_s'$ as fucntion of $E_s$ and $E_w$?

Thanks in adavance.

PS: I have already posted a smiliar question here
https://math.stackexchange.com/questions/1560821/which-is-the-total-energy-of-the-product-of-two-discrete-energy-signals

You don't give enough information for an expression of the combined energy in your expression.
$$\begin{eqnarray}
& &\sum_{i=0}^N \left| x_s(i)\cdot\left( 1+K\,w(i) \right) \right|^2 \\
&=&\sum_{i=0}^N \left( x_s(i)\cdot\left( 1+K\,w(i) \right) \right)^\star \times\left( x_s(i) \cdot\left( 1+K\,w(i) \right) \right)\\
&=&\sum_{i=0}^N \left( x_s(i)^\star\,x_s(i) + x_s(i)^\star\,K^\star\,w(i)^\star x_s(i)\,K\,w(i) + x_s(i)^\star x_s(i)\,K\,w(i)+x_s(i)\,x_s(i)^\star\,K^\star w(i)^\star\right)\\
&=&\sum_{i=0}^N \left( |x_s(i)|^2 \times \left(1+2\,\Re(K\,w(i))+|K|^2\,|w(i)|^2 \right) \right)\\
&=&E_s+\sum_{i=0}^N|x_s(i)|^2\left(2\,\Re(K\,w(i))+|K|^2\,|w(i)|^2 \right)
\end{eqnarray}$$

This means you need to know the correlation of the square magnitude of $x_s$ with the real part of $w$ and the square magnitude of $w$. The correlation of just $x_s$ and $w$ does not play any role.


Correcting phase response of a signal I have a sensor that applies a frequency dependant phase alteration to my signals.  I'm trying to find a way to correct these phase responses using an FIR filter.  Essentially, I'm looking for methods to apply an all-pass filter with arbitrary frequency response to correct this.
For magnitude calibrations of my sensor, I can achieve this by drawing out my calibration spectrum, applying an ifft to it, then applying a hamming window.  Then I simply convolve the resulting impulse response with my signal.  This is essentially MATLAB's fir2().  Is there an analogous procedure I can do to correct arbitrary phase responses in my signals?

You can equalize magnitude and phase simultaneously by defining a desired complex frequency response
$$D(\omega)=M(\omega)e^{j\phi(\omega)}\tag{1}$$
with magnitude $M(\omega)$ and phase $\phi(\omega)$ chosen such that they compensate for the given magnitude and phase distortions.
An FIR filter approximating $(1)$ can be designed by using the following error function:
$$E=\sum_kW_k\big|D(\omega_k)-H(\omega_k)\big|^2\tag{2}$$
where $\omega_k$ are frequencies on a dense grid, $W_k$ are some non-negative weights allowing to emphasize certain frequency regions compared to others, and $H(\omega)$ is the frequency response of the FIR filter to be designed:
$$H(\omega)=\sum_{n=0}^{N-1}h[n]e^{-jn\omega}\tag{3}$$
Minimizing $(2)$ with respect to the filter coefficients $h[n]$ can be achieved by solving a system of linear equations. I've written a Matlab/Octave function that does just that: lslevin.m
Note that it may be necessary to add a linear phase (i.e., a delay) to the required $\phi(\omega)$ in order for a causal filter to better approximate the desired complex frequency response. As a guideline, the average delay implied by the desired phase should be about half the chosen filter length $N$. It is worthwhile experimenting with different additional delays because they may have a significant influence on the resulting approximation error.


How to tell whether an image has a low gamma (Finding out imaging system tone map) I am using Emgu C#.

I acquire images from my cameras.

I know when there is no saturation (ie a grey-tone image) by looking at the average hue channel of the image. If zero then pretty sure picture taken at night time.

I can also detect how bright an image is by checking the average Luminescence channel.

What I would now like to do is a value for the gamma an contrast of an image.

Now, I know it is the interpretation of this hypothetical gamma and contrast that determines its value within the context of how I use it. But what factors can i use to attempt to put a scale on these values?

Should I be looking at a different Color Space?

Thanks

You should take a photograph of an object with known color values, such as Macbeth chart:

By using the lowest row of grayscale colors, you will be able to find out the tone-mapping (sometimes incorrectly called gamma) of your imaging system.


Recognition of markes in long audio recording I am about to do a lot of a recordings when performing tests. The thing is that there are some parts I would like to extract easily from my recording. Is there a way to play some special sequence of sounds, for example few beeps, and then search the whole recording for those markers?

In the end I want to obtain file with labeled segments similarly to what Audacity is doing. Until now I was using hand claps and searching for energy spikes in the long recorded signal. Since I am performing measurements of room with starter pistol this is not a good idea. Later this will allow me for very quick and easy analysis to search for groups of similar sounds. For example:

  start  time, end time, location

  100, 399, loc1

  500, 600, loc2

I tried calculating correlation with some pattern, but when room is reverberant it is not performing very good. I get some false markers from threshold correlation. It is because marker is smeared and not exact.

What would be the best method and type of a signal for that task? Any good and fast techniques such as matched filtering?

Addendum:

saying simple what i want:

play the marker sound/pattern at random times

record this signal and other signals

search for these recorded patterns in very long recording
being able to say at which point in time those patterns are
for example cross correlation (isn't working well) could give spikes i can search for
it should be idiot-proof against noise, tonal/impulsive sounds, reverberation and low playback level
it should not be very slow. for example at 10 minutes file 40 seconds is my upper limit

You can use DTMF tones which are very easily decodable with resonator banks. Brief explanation at http://en.m.wikipedia.org/wiki/Dual-tone_multi-frequency_signaling

If you are in a reverberant room then something like a mobile phone will not put too much energy in it and consequently the recording will not be affected too much by the reverberation. Also, DTMFs are relatively high pitch and would stand out from the background.

There are 16 combinations in a typical keypad but then you can alter symbols (e.g. 1213), possibly in a recording, to build more complex phrases or markers. Altering the symbols quicker than the time constant of the room will reduce the effect of reverberation.

If you are using the starter gun to obtain impulse responses from the room, you can play a sequence in the beginning and just after the end of the recording to roughly identify a particular response but I wouldn't recommend relying on the markers for absolute syncing. You can add a mic close to the gun to obtain the pulse start and end and use the room mics to estimate the dynamics of the room (attack, decay and others)

Hope this helps

Updated part of this answer:

The general idea remains the same, tag a point in time by playing a sound and detect its presence later. The question now is, which sounds and how to detect them. Here is an example using DTMF codes.

1) Get a DTMF generator (for example, a DTMF android app such as this one, https://play.google.com/store/apps/details?id=hu.soska.dtmf&hl=en )

2) In a "clean" (non-reverberant) room, record a sequence of sounds, for example "1928". This is the "tag". Please make sure that each number is keyed in for at least 1 second.

3) In the room under measurement, start your recording. When you want to tag a point, have the "1928" sound reproduced and recorded along with the main recording.

4) Build a filter bank of 8 filters at the DTMF frequencies. The resonators are dead simple filters, here is some example code to get your coefficients: http://lifeorange.com/MATLAB/MATLAB_FD.htm

5) Open the file, run the recording through the filters. Their output will vary depending on the power of the tone in the signal. During the reproduction of DTMF tones, the output of the filters will be increasing. Apply a very simple rule: To maintain a "correct match" you would have to see pairs of filters (representing numbers) going "high" and staying high for at least 1 second. Look for this rule in the outputs of your filters. You now have the matches.

Why single tones and resonators? Because they are very easy and fast to decode.

Alternatives?

Record a Pseudorandom generator signal at step 4 and use cross correlation for step 5. Cross correlation though will be much slower because it would have to run repeatedly in a "sliding window" fashion (please see http://en.wikipedia.org/wiki/Overlap%E2%80%93add_method).

An alternative for step 5 could be multimon (please see http://manpages.ubuntu.com/manpages/gutsy/man1/multimon.1.html ), but that would not give you the offset.

These are all approximate solutions though, they will not give you a dead accurate point in time that your "event" occurs, only a rough approximation (which also depends from the length of your recording).

Hope this helps.


Analyzing movement of fringes I have taken a series of images of a fringe pattern at regular time intervals. This fringes are generated by shining a laser beam onto the CCD camera. The laser goes through a lot of optics, which is why it creates these fringes.

There is some vibration in the laser system setup which causes these fringes to move around although not by much, couple of pixels at most. I want to analyze how much these fringes are moving around. Is there an easy way to do this?

One idea that I had, and I don't know if it's right, is to take the FFT of the image, and discard the DC (constant) components. I could analyze just one frequency component. Now, if the fringe pattern is moving about, then that should change the phase of that frequency component. Does that make sense?

Thanks.

Update
As suggested by A_A, I took one frame as a reference and subtracted it from the rest:

Out of the mess of fringes, I can clearly identify one pattern which is moving about. I particular, the diagonal fringes do not seem to be changing from frame to frame. In the end, I just want to identify if there's any one particular optic which is responsible for this vibration, so I could just set up my live camera feed to display the difference between consecutive images and play around with the optics to see if I can reduce the fringes.

...take the FFT of the image, and discard the DC (constant) components.

That would result in removing the mean brightness from an image. Due to the fact that grayscale images have positive numbers between 0-255, what would be returned from that "high pass" filter (that only "cuts" the DC) would be the image with the same grayscale variance but now centered around 0 (instead of its original mean level).

  I want to analyze how much these fringes are moving around. Is there an easy way to do this?

You could obtain the absolute difference between two successive frames and then the sum of all values of that quantitty as a simple metric of how much movement was there.

If you would also like to be able to estimate how much displacement was there between two successive frames then the simplest thing would be to use cross correlation (on successive frames) and track the position of its maximum.


MATLAB : Proper estimation of weights and how to calculate MSE for QPSK signal for Constant Modulus Algorithm I had asked conceptual Question about Constant Modulus Algorithm. I am implementing the simple steps of the algorithm, as I do not have the cma() built in module.

I am considering a Finite Impulse Response (FIR) system  whose true coefficients are $h = \begin{bmatrix}1 &0.45 &-0.2\end{bmatrix}$. The algorithm is presented briefly :

$u_n$ : Output of the FIR filter that is driven by Input signal, $s$ is a white Gaussian signal that drives an FIR process by the command

s=round(rand(1,N))*2-1;
u = filter(h,1,s);

We receive a noisy signal which is the output of the system corrupted by AWGN. Let the noise corrupted signal be

$$x_n = h^T s_n + \eta_n$$
Constructing an equalizer

$$y_n = w^T x_n$$

The cost function to be minimized by gradient descent is

$$J(w) = E\left[\left(\big\lvert y_n\big\rvert^2 -1\right)^2\right]$$

The weight update equation is given by

$$w_{n+1} = w_n - 2\mu e_n  y_n^T x_n$$

where error,

$$e_n = \left(\big\lvert y_n\big\rvert^2 -1\right)$$

Operator $T$ is the transpose and assuming real signals without any imaginary and complex parts.

Below is the code. The algorithm does not give proper result in terms of estimates of the weights.

Problem 1: The Graph for the Mean Square Error (in DB) for the estimated weights returned by the algorithm vs. Signal to Noise Ratio is giving an opposite trend i.e. instead of MSE decreasing with increasing SNR, I am getting the opposite when I re-run the program!! Below is the image of what I mean. The first figure is correct, but with the same code, I ran it again and got the second figure. Why is this & how can I prevent it?

I am not sure if this is due to the initialization of the variable $L$, filter order = number of delays = 2. Somethings are not straight for me which are for the FIR filter of the form : $u(t) = e(t) + 0.45e(t-1) - 0.2e(t-2)$, what is the filter order, smoothing length and the number of weights unknown to be estimated?
 Can somebody please help so that it works?

   clear all
clc

N = 256;
h = [1    0.45   -0.2];

R2 = 2;

mu =  1.0000e-009;
noisedB  =0;

L=2; % smoothing length L+1
ChL=1;  % length of the channel= ChL+1
EqD=round((L+ChL)/2);  %  channel equalization delay

i=sqrt(-1);

Ch=[1 0.45 -0.2]; %Channel
%Ch=[0.8+i*0.1 .9-i*0.2]; %complex channel
Ch=Ch/norm(Ch);% normalize

skip =1
for l=1:6
i          = 1;

 TxS=round(rand(1,N))*2-1;  % QPSK symbols are transmitted symbols
%TxS=TxS+sqrt(-1)*(round(rand(1,N))*2-1);
x=filter(Ch,1,TxS); %channel distortion

n=randn(1,N);   % additive white gaussian noise

n=n/norm(n)*10^(-noisedB/20)*norm(x);  % scale noise power
x1=x+n;  % received noisy signal

%estimation using CMA
K=N-L;   %% Discard initial samples for avoiding 0's and negative
X=zeros(L+1,K);  %each vector
for j=1:K
    X(:,j)=x1(j+L:-1:j).';  %y_n = w^T x_n
end

e=zeros(1,K);
w=zeros(L+1,1);
w(EqD)=1;    % initial condition
while i<=K

e(i)=abs(w.'*X(:,i))^2-R2;                  % initial error
   w=w-mu*2*e(i)*X(:,i)*X(:,i)'*w;  % update equalizer co-efficients
    cma_mse_h(l,i)    = sum((w'-h).^2)/3;
    est_w(i,:) = w;
   w(EqD)= 1;

            i = i+1;

end
noisedB = noisedB + 5;
end

for ii = 1:6

 Error(ii) = 10*log10(mean(cma_mse_h(ii,:)));
end

 plot([0:5:25], Error(1:6));
 grid on;
 xlabel ('SNR(dB)')
 ylabel('MSE_h')

Second Problem : The true FIR channel coefficients are not imaginary and have only real parts. But the estimated weights will have both real & imaginary if I work with real and complex representation. How can I properly calculate the weights & its MSE for this case?

Your code reveals many misconception about what the CMA is supposed to achieve:

your step size mu is much too small; note, however, that the optimal step size can only be found through experiment.
the variabe noisedB appears to be the desired SNR of the received signal. An SNR of $0\,\text{dB}$ as specified by you is very poor (the noise is as strong as the signal), and in such noisy conditions the CMA cannot perform properly.
you have a channel with $3$ taps, and your equalizer length is also only $3$ taps; this is much too short to achieve any reasonable equalization.
you define the MSE as the error between the channel impulse response and the equalizer response. This is wrong, because the equalizer is not supposed to converge to the channel impulse response! The equalizer must equalize the channel, i.e. the concatenation of the channel and the equalizer should give a pure delay.

There are a few other small errors in your code, but instead of correcting them all I show you a very simple and short code example. The figure below the code shows the error curve vs symbol index, and the total impulse response of the channel and the equalizer. As you can see, it approximates a delta impulse (i.e. no distortion). The SNR is $30\,\text{dB}$ for this example. The error becomes much smaller if there is no noise. For worse SNRs, the algorithm becomes unstable.

N = 20000;          % # symbols
h = [1,.45,-.2];    % channel impulse response
h = h/norm(h);
Le = 20;            % equalizer length
mu = .001;          % step size
snr = 30;           % snr in dB

s0 = round( rand(N,1) )*2 - 1;     % BPSK signal
s = filter(h,1,s0);                % filtered signal

% add Gaussian noise at desired snr
n = randn(N,1);
vs = var(s);
vn = vs*10^(-snr/10);
n = sqrt(vn)*n;
r = s + n;          % received signal

e = zeros(N,1);     % error
w = zeros(Le,1);    % equalizer coefficients
w(Le)=1;            % actual filter taps are flipud(w)!

yd = zeros(N,1);

for i = 1:N-Le,
    x = r(i:Le+i-1);
    y = w'*x;
    yd(i)=y;
    e(i) = y^2 - 1;
    w = w - mu * e(i) * y * x;
end

np = 100;           % # sybmols to plot (last np will be plotted); np < N!

subplot(3,1,1), plot(e.*e), title('error')
subplot(3,1,2), stem(conv(flipud(w),h)), title('equalized channel impulse response')
subplot(3,1,3), plot(1:np,s0(N-np+1:N),1:np,yd(N-np+1-Le+1:N-Le+1))
    title('transmitted and equalized signal'), legend('transmitted','equalized'), axis([0,np,-1.5,2])


Typical parameters and stability of Chebyshev filters I design Chebyshev first and second type filters using scipy.signal.iirfilter.

What are typical values of maximum ripple in the passband and minimum attenuation in the stopband? I struggle to find optimal values of those parameters.
I know that IIR filters are often unstable, but is it possible to make it stable by changing values of maximum ripple in the passband and minimum attenuation in the stopband? If not, are there any ways to make it more stable, say, for the 6th order filter?

I design Chebyshev first and second type filters

You are not going to see a lot of effect with this low order filters. If you want to study ripple and attenuation crank it up.

  What are typical values of maximum ripple in the passband and minimum attenuation in the stopband? I struggle to find optimal value of those parameters.

This really depends on your application. You can pick something easy for starters: 1 dB ripple and 40 dB attenuation. That's probably not useful for a real application but at least you can make some graphs

  I know that IIR filters are often unstable,

IIR filters are only unstable if the poles are on or outside the unit circle. Chebyshev are inherently stable unless you made a mistake. Implementation of a stable filter can indeed be unstable. That is heavily dependent on filter representation, topology, numerical accuracy and pole locations.  A 6th order filters in scipy should not be a problem unless you are trying to use fix point or the cutoff frequencies are extremely low.


Does the output of this system have a closed form expression? \begin{align}
x(n)&=\left(-1\right)^n u(n)\\
h(n)&= u(n-1)+u(n-2)\\
y(n)&= x(n)*h(n)\\
y(n)&= \sum_{k=-\infty}^\infty \left(-1\right)^k u(k)\left[u(n-k-1)+u(n-k-2)\right]
\end{align}
Graphically it looks like the output is $u(n-1)$ but I can't do the math to get it into the closed form.

Edit: Ok I pretty much worked it through I'm just iffy on justifying $n \ge 1$

$u(k)$ sets the bottom limit of the summation to $0$ and $u(n-k-1)$ sets the upper limit to $n-1$. $n$ has to be at least $1$ for the summation to make sense and the partial sums = $1$.

You figured out the solution by yourself, so this is just about how to show it formally. First, split the sum into two sums:

$$y[n]=\sum_{k=-\infty}^{\infty}(-1)^ku[k]u[n-1-k]+\sum_{k=-\infty}^{\infty}(-1)^ku[k]u[n-2-k]\tag{1}$$

The unit steps give you upper and lower summation indices, but note that the resulting sums are only valid if the upper summation index is not smaller than the lower one:

$$y[n]=u[n-1]\sum_{k=0}^{n-1}(-1)^k+u[n-2]\sum_{k=0}^{n-2}(-1)^k\tag{2}$$

Now note that you can increase the upper summation index by any even number without changing the value of the sum. The second term on the right-hand side of $(2)$ can thus be rewritten as

$$u[n-2]\sum_{k=0}^{n-2}(-1)^k=u[n-1]\sum_{k=0}^{n}(-1)^k\tag{3}$$

In $(3)$ the argument of the unit step can be changed to $n-1$ because the value of the sum equals zero for $n=1$. Using $(3)$, Eq. $(2)$ can be rewritten as

$$\begin{align}y[n]&=u[n-1]\left(\sum_{k=0}^{n-1}(-1)^k+\sum_{k=0}^{n}(-1)^k\right)\\&=u[n-1]\left(2\sum_{k=0}^{n-1}(-1)^k+(-1)^n\right)\\&=u[n-1]\left((1+(-1)^{n+1})+(-1)^n\right)\\&=u[n-1]\tag{4}\end{align}$$


Help with Audio Signal noise removing using FFT and IFFT in Matlab I dont understand why my modified fft of a audio signal returns complex values after ifft.
I modified my signal by zeroing out all the unnecessary freq.

FFT: Fast fourier transfrom

IFFT: Inverse fast fourier transform

Heres the matlab code and graph

[a,fs]=audioread("Recording.m4a")

plot(a)

n=numel(a)/2  *//*since audio signal is stereo but i want only mono*

b=abs(fft(a))

y=fft(a)

c=y(:,1)     *//*to convert audio signal from stereo to mono*

plot(f,abs(c))

xlim([0 1000])

c(abs(c)<31)=0

c(round(n*1000/fs):end)=0

f=0:n-1

f=(f*fs/n)'

plot(f,abs(c))

xlim([0 1000])

z=ifft(c)

plot(z)

The short answer is do not null out the "mirror frequencies" that are located above $f_s/2$ that match the frequencies you want to keep. If your FFT was generated from a real signal, then when you do the IFFT you will get the real signal back as long as you did not zero out those upper frequencies (as you did).

The DFT (which the FFT computes) returns complex values but for real signals the “positive” frequency values will be the complex conjugate of the negative frequency values which together represent real signals. Due to how the FFT is computed, instead of returning the positive and negative frequencies from $-f_s/2$ to $+f_s/2$, instead those same frequencies are given from $0$ to $f_s$.

Consider Euler’s Identity:

$$2cos(\omega t) = e^{j\omega t}+ e^{-j\omega t}$$

If it wasn’t already clear, the general form of  $Ae^{j\theta}$ represents the magnitude and phase of a complex number as $A\angle{\theta}$ so above we see how the cosine is represented by two phasors each rotating in time in opposite directions (a positive and negative frequency) such that their sum is always on the real axis.

The FFT returns the frequency bins from 0 to one sample less than the sampling frequency: $n =0$ to $N-1$ bins where each bin is spaced by $f_s/N$ with $f_s$ as the sampling rate. Due to the cyclical nature of the FFT this represents the positive and negative frequencies by mapping everything above $f_s/2$ to the negative frequencies, specifically samples $N/2 + n$ for $n\ge N/2$ map to $-N/2+n$. The command fftshift in MATLAB or Octave does this specifically.

For this reason you should NOT null out all the "mirror" frequencies that are located above $N/2$ as without them your signal will be complex. Consider the example above with the cosine.

$$2cos(2\omega_n t) = e^{j2\omega_n t}+ e^{-j2\omega_n t}$$

If sampled such that there were 100 total bins, and $\omega_n$ was the spacing on one DFT bin, then the DFT would have non zero values at the 3rd bin (bin 0 is DC) and the 2nd to last DFT bin which after shifting would be at the expected +/-2 bins. If you removed the higher bin then the resulting waveform would simply be:

$$e^{j2\omega_n t}$$


Confusion in FFT of Resampled Data I have the following code in Matlab:

fs = 1e6;

y = randn(1,100);
y = upsample(y,8);
windowSize = 8;
b = (1/windowSize)*ones(1,windowSize);
a = 1;
y = filter(b,a,y);

fft1 = fft(y);
fft1 = fft1(1:length(y)/2+1);
psd = (1/(fs*length(y)))*abs(fft1).^2;
psd(2:end-1) = 2*psd(2:end-1);
f = 0:fs/(length(y)):fs/2;
figure(1);
plot(f,10*log10(psd))

y = resample(y,7,4);
fft1 = fft(y);
fft1 = fft1(1:length(y)/2+1);
psd = (1/(fs*length(y)))*abs(fft1).^2;
psd(2:end-1) = 2*psd(2:end-1);
f = 0:7*fs/(4*length(y)):7*fs/8;
figure(2);
plot(f,10*log10(psd))

In this code I have a random noise signal and I am upsampling it and passing it through a zero order hold filter then plotting its Power spectral Density (PSD) in figure 1.
After resampling the signal to a new rate when I check its PSD in figure 2 I get some different result.
Could someone explain why is the spectrum of the resampled data different from the original one and what can I do to get the same PSD as the original signal?

As you know you have to use an ideal low-pass filter, and zero order hold is an approximation which needs more compensation. Ideal LPF has a flat response with sudden cut-off and zero stop band but your moving average window has a sinc frequency response. So the deep valleys belong to the zeros of sinc function and the remaining which has high attenuation passed from the stopband of your filter (which is not zero because it's not an ideal LPF).


One approach would be to cast the problem as least-squares smoothing. The idea is to locally fit a polynomial with a moving window, then evaluate the derivative of the polynomial. This answer about Savitzky-Golay filtering has some theoretical background on how it works for nonuniform sampling.

In this case, code is probably more illuminating as to the benefits/limitations of the technique. The following numpy script will calculate the velocity and acceleration of a given position signal based on two parameters: 1) the size of the smoothing window, and 2) the order of the local polynomial approximation.

# Example Usage:
# python sg.py position.dat 7 2

import math
import sys

import numpy as np
import numpy.linalg
import pylab as py

def sg_filter(x, m, k=0):
    """
    x = Vector of sample times
    m = Order of the smoothing polynomial
    k = Which derivative
    """
    mid = len(x) / 2
    a = x - x[mid]
    expa = lambda x: map(lambda i: i**x, a)
    A = np.r_[map(expa, range(0,m+1))].transpose()
    Ai = np.linalg.pinv(A)

    return Ai[k]

def smooth(x, y, size=5, order=2, deriv=0):

    if deriv > order:
        raise Exception, "deriv must be <= order"

    n = len(x)
    m = size

    result = np.zeros(n)

    for i in xrange(m, n-m):
        start, end = i - m, i + m + 1
        f = sg_filter(x[start:end], order, deriv)
        result[i] = np.dot(f, y[start:end])

    if deriv > 1:
        result *= math.factorial(deriv)

    return result

def plot(t, plots):
    n = len(plots)

    for i in range(0,n):
        label, data = plots[i]

        plt = py.subplot(n, 1, i+1)
        plt.tick_params(labelsize=8)
        py.grid()
        py.xlim([t[0], t[-1]])
        py.ylabel(label)

        py.plot(t, data, 'k-')

    py.xlabel("Time")

def create_figure(size, order):
    fig = py.figure(figsize=(8,6))
    nth = 'th'
    if order < 4:
        nth = ['st','nd','rd','th'][order-1]

    title = "%s point smoothing" % size
    title += ", %d%s degree polynomial" % (order, nth)

    fig.text(.5, .92, title,
             horizontalalignment='center')

def load(name):
    f = open(name)
    dat = [map(float, x.split(' ')) for x in f]
    f.close()

    xs = [x[0] for x in dat]
    ys = [x[1] for x in dat]

    return np.array(xs), np.array(ys)

def plot_results(data, size, order):
    t, pos = load(data)
    params = (t, pos, size, order)

    plots = [
        ["Position",     pos],
        ["Velocity",     smooth(*params, deriv=1)],
        ["Acceleration", smooth(*params, deriv=2)]
    ]

    create_figure(size, order)
    plot(t, plots)

if __name__ == '__main__':
    data = sys.argv[1]
    size = int(sys.argv[2])
    order = int(sys.argv[3])

    plot_results(data, size, order)
    py.show()

Here are some example plots (using the data you provided) for various parameters.

Note how the piecewise constant nature of the acceleration becomes less obvious as the window size increases, but can be recovered to some extent by using higher-order polynomials. Of course, other options involve applying a first derivative filter twice (possibly of different orders). Another thing that should be obvious is how this type of Savitzky-Golay filtering, since it uses the midpoint of the window, truncates the ends of the smoothed data more and more as the window size increases. There are various ways to solve that problem, but one of the better ones is described in the following paper:

  P.A. Gorry, General least-squares smoothing and differentiation by the
  convolution (Savitzky–Golay) method, Anal. Chem. 62 (1990) 570–573. (google)

Another paper by the same author describes a more efficient way to smooth nonuniform data than the straightforward method in the example code:

  P.A. Gorry, General least-squares smoothing and differentiation of
  nonuniformly spaced data by the convolution method, Anal. Chem. 63
  (1991) 534–536. (google)

Finally, one more paper worth reading in this area is by Persson and Strang:

  P.O. Persson, G. Strang, Smoothing by Savitzky–Golay and Legendre
  Filters, Comm. Comp. Finance 13 (2003) 301–316. (pdf link)

It contains much more background theory, and concentrates on error analysis for choosing a window size.


Frequency content Signal processing of very large files I recently had to process a very large file to obtain the power spectrum. The file had about 60m measurements (arguably its not that huge, but its getting to that point). Although, in the end (initially I tried with Labview) I managed to find a few tools (Matlab, R) that allowed me to perform the calculation I needed on my PC, I couldn't help but wondering what would be the mathematically correct way to break up that problem, so that I can perform it piecewise.
I guess, my question boils down to:

Can I perform FFT on pieces of a very large signal, and then assemble the results to obtain the FFT of the total signal?

Bonus question:

What happens with the periodogram?

The way I interpret the periodogram is that its an fft on a smaller portion of data. So, would I get a much different result if I took a subsets of the signal (overlapping) and performed fft and the created the periodogram?

Update of information on signal No:2 (500kHz):
In advance apologies for changing the sampling rate. I just got some more results with higher sampling rate because I suspect the noise is so high frequency that cannot be captured.

Signal Description:

The signal comes from measurements in a wind tunnel.
The recorded values are the electrical signal of pitot tube, when the wind tunnel speed is set at a specific value.
The reason, the measurements were taken was because noise was creeping in (probably electrical), the signal and there is a need to investigate the frequencies that contribute so that an appropriate filter is applied (one that does not affect the measurement of wind turbulence).

Metadata:

Data acquisition rate: 500kHz (I just got another set of data with 500kHz)
Data recording rate: 500kHz
Recording duration: approximately 1 min (this is the typical length).

Preprocessing:

The raw signal is recorded with no preprocessing.
For the FFT the only preprocessing involves removing the mean value from the signal (so as to remove the DC component). (A window is not applied).

Sample of data

Following are two excerpts (1000 data points from a 500kHz data).
The first is the data with the suspected source of the noise turned on (it is actually the inverter that power the wind tunnel).  Although the inverter is turned on there wind speed is set to zero.

The second set, is an excerpt with the inverter turned off (obviously the wind speed is set to zero), at the same scale on the y axis.

The following graph is exactly the same as the second with autoscale on the y axis.

Longer interval excerpt 0.06[s] (30000 points)

This is to give a better overview of the signal. I switched to points (instead of lines) to give a more clear picture. The following two graphs are again with the electric noise a) on and b) off

Can I perform FFT on pieces of a very large signal, and then assemble
the results to obtain the FFT of the total signal?

This is the basis for the FFT algorithm, in that a large DFT can be more efficiently done as 2 DFT's each half the size (and so on and so on until the only thing left is 2 point DFT's).
(The efficiency comes about since a DFT requires on the order of $4N^2$ multiplications and $4N^2$ additions, while an $N$ pt DFT derived from two $N/2$ pt DFT requires only $4(N/2)^2 + 4N$ multiplications and $4(N/2)^2 + 4N$ additions.)
Below shows the combining approach for splitting the DFT into two separate DFT's, and the same approach can be continued further back dividing the blocks in half each time to whatever minimum size is desired.

The graphic above is given by the mathematical equation below showing the divide and conquer reduction as the basis for the FFT algorithm.
$$\sum_{n=0}^{N-1}x[n]W_N^{nk} = \sum_{r=0}^{N/2-1}x[2r]W_{N/2}^{rk} + W_N^k\sum_{r=0}^{N/2-1}x[2r+1]W_{N/2}^{rk}$$
Where $W_N^k = e^{-j2\pi k/N}$
Alternatively, individual FFT blocks could be averaged, but in this case it would be the magnitudes that would be averaged to result in a spectral smoothing but that would not improve the spectral resolution the way increasing the total time duration of the DFT as the process above provides. Averaging the blocks would be equivalent to adjusting the video bandwidth on a spectrum analyzer (post detection filtering), in that the noise floor would be smoothed, but not reduced, while increasing the total FFT time duration would be equivalent to adjusting the resolution bandwidth on a spectrum analyzer.


Why do we need >= 95% throughput for LTE Sensitivity measurement according to 3GPP According to 3GPP 36 series for LTE, for LTE Sensitivity measurements we need >=95% throughput. How is this 95% value determined? Why it is not 98% or 88%.
Thanks in Advance!

It's a definition. There's probably no specific reason
other than "the committee agreed on it".
Point is: this is for comparability of other tests. So, the actual number 95% is less important than that all use the same.


aliasing in image processing I know that aliasing occurs when a signal is subsampled. If the sampling rate is lower than twice the max frequency in a signal, aliasing occurs.
How is it in pictures? as far as I know, a sinc-filter is not applicable in photos. so my question is, is aliasing always in pictures?

i see aliasing happening when the image is captured with a sensor that cannot capture higher frequencies.

Exactly!

If i capture a image, is it always aliased?

Not when the image content is sufficiently bandlimited.
That's often the case because of the physics behind things, or simply because the motive isn't higher in frequency.
(Physical effects include motion blur, abberation, depth of field/ring of confusion/bokeh, wave effects on apperture edges, lens imperfections, air density variations and many more. I'd recommend you look for tutorials on how to take really sharp macro and tele photos to learn what effects can reduce your image bandwidth.)
But, sure. Go photograph a fine grid, and you see a Moiré pattern. By the way, your eyes have limited resolution, too, although I'd presume the brain largely adapts that away.

Other types of aliasing are as far as i understand not really aliasing. e.g if a screen has not many pixels, the image is aliased.

You call it aliasing, it aliases things, you describe it as aliasing: It's really aliasing.


Polyphase Filter decomposition. It is not working So, I want to do a polyphase implementation of a filter bank (I had problems with it since the beginning and I already asked for help to the community in the past, indeed: Norm MPEG-1 Layer III (Mp3) Filter banks distortion). Thanks to the support I was able to make it work (without polyphase implementation).

Nonetheless, I tried to do all the analysis filters and synthesis filters by means of polyphase filters banks as follows:

From:

To:

The idea is, basically, based on the first picture, with the low pass filter prototype modulated by a cosine to make frequency shifts of the filter response to cover all the bandwidth (with 32 filters), decomposing each particular filter response into 32 polyphase branches.

(So 32 branches, obtained by means of shifting the prototype by multiplying a cosine. Each of these 32 branches is implemented by means of decomposing each analysis and synthesis filters into 32 polyphase branches).

So, once the first part works (i.e. The filter bank not implemented by means of polyphase filter banks), I do the polyphase decomposition as follows:

1. Polyphase decomposition:

h are the modulated coefficients of each filter of the bank.
hh is the polyphase filter bank coefficients
the $i^{\rm th}$ row of the matrix represents the coefficients of the $i^{\rm th}$ filter of the polyphase bank
M are the polyphase filter branches, 32 in our examplehh=zeros(M,length(h)/M);
for l=0:1:M-1
n=l+1:M:length(h);
hh(l+1,:)=h(n);
end
end

2. The shift ($Z^{-1}$) is performed by a function that performs as follows ("delay" is the delay to apply. e.g. if 0 samples are to be delayed, then the output shall be equal to the input):

in_data: Input sequence
out_data: Output sequence.
delay_samp=delay+1;
out_data=zeros(1,length(in_data));
out_data(delay_samp:length(in_data))=in_data(1:length(in_data)-(delay_samp-1));

3. Ok, now, I apply the filters as follows:

Analysis:

yy=zeros(branches,length(hh(1,:))+fix(length(in_signal)/branches));

for k=1:branches

    %1. Shift
    shifted_input=shift_data(in_signal,k-1);

    %2. Downsample
    temp=zeros(1,fix(length(in_signal)/branches));
    downsampled_data=shifted_input(1:branches:end);

    temp(1:length(downsampled_data))=downsampled_data;
    shifted_downsampled_input=temp;

    %3. Filter
    temp=conv(hh(k,:),shifted_downsampled_input);
    yy(k,:)=temp;

end

And Synthesis:

yy=zeros(branches,branches*(length(hh(1,:))+length(in_signal)-1));

for k=1:branches

    %1. Filter
    temp_conv=conv(hh(k,:),in_signal);

    %2. Upsample
    temp=zeros(1,length(temp_conv)*branches);
    temp(1:branches:end)=temp_conv;%upsample(temp,branches);

    %3. Shift.
    yy(k,:)=shift_data(temp,branches-(k-1));

end

The result is very distorted:

The output with a 1KHz tone as an input:

I guess that I am making something wrong with the upsampling and downsampling processes as there are spectral replicas all over the second example, but I have no clue on what exactly is wrong.

I checked on the internet and I also checked the reference referred in the other post. I mean this book. But I think there is a problem with the way I am coding this not with the theoretical back.

This is the low pass filter prototype.

The following is a working code that uses 32-component polyphase decomposition of the associated 32-channel anslysis and synthesis filterbanks. As I have already commented, the speed gain is not dramatic in this cae due to short signal and filter lengths. However further architectural improvements as well as coding optimizations can provide better results.

% S0 - Load the prototype lowpass filter impulse response h0[n]:
% --------------------------------------------------------------
load h2.mat;         % h[n] is the prototype lowpass filter of length 512
L = length(h);

% S1 - Create the 32 x 512 analysis filter-bank hha[k,n] by cosine modulation from protoype :
% -----------------------------------------------------------------------------------------
numbands = 32;                 % number of banks (channels)
n=0:L-1;

hha=zeros(numbands,L);         % bank of filters hha[k,n] = 32 x 512 array.
for k=0:1:numbands-1
   hha(k+1,:) = h.*cos( ( (2*k+1)*pi*(n-16) ) / (2*numbands) );
end

% S2 - Create the 32-polyphase components hhap[k,m,n] , for each one of 32 analysis filters hha[k,n]:
% ---------------------------------------------------------------------------------------------------
numpoly = numbands;             % polyphase component number = decimation ratio = number of channels
hhap = zeros(numbands,numpoly, L/numpoly);  % hhap = 32 x 32 x 512/32 , 3D ANALYSIS filter bank array

M = numpoly;                    % polyphase system decimation ratio
for k=1:numbands
    for m = 1:numpoly
        hhap(k,m,:) = hha(k,m:M:end);       % create the m-th polyphase component of k-th channel filter
    end
end

% S3 - Design the 32 x 512  synthesis (complementary) filter bank :
% -----------------------------------------------------------------
numbands = 32;                  % number of banks
n=0:L-1;
hhs = zeros(numbands,L);        % bank of filters
for k=0:1:numbands-1
   hhs(k+1,:) = h.*cos( ( (2*k+1)*pi*(n+16) ) / (2*numbands) );
end

% S4 - Obtain the 32-polyphase components hhsp[k,m,n] , for each one of 32 synthesis filters hhs[k,n]:
% ----------------------------------------------------------------------------------------------------
numpoly = numbands;             % polyphase component number = interpolation ratio = number of channels
hhsp = zeros(numbands,numpoly, L/numpoly);  % hhap = 32 x 32 x 512/32 , 3D ANALYSIS filter bank array
M = numpoly;                    % polyphase system decimation ratio
for k=1:numbands
    for m = 1:numpoly
        hhsp(k,m,:) = hhs(k,m:M:end);       % create the m-th polyphase component of k-th channel filter
    end
end

% S5 - Generate the test input signal
% -----------------------------------
N = 2*1024;
wav_in = cos(0.01791*pi*[0:N-1]);        % pure sine tone

% S6 - Apply test signal to the filterbank: ANALYSIS STAGE :
% -----------------------------------------------------------
yyd = zeros( numbands, floor(N/numbands));   % decimated outputs..
M = numbands;
for k=1:1:numbands
    temp = conv([wav_in(1:M:end),0] , hhap(k,1,:));
    for m=2:M
        temp = temp + conv([0,wav_in(M-m+2:M:end)],hhap(k,m,:));
    end
    yyd(k,:) = temp(L/(2*M)+1 : L/(2*M)+N/numbands);
end

% S7 - Apply SYNTHESIS filterbank on the decimated signal :
% ---------------------------------------------------------
ys = zeros(1, N);

for k=1:numbands
    temp = zeros(1, N+L-1);
    for m = 1:numpoly
        temp(m:numbands:end-31) = conv( yyd(k,:) , hhsp(k,m,:) );
    end

    ys = ys + temp(L/2+1:L/2+N);
end
ys = numbands*ys;

% SX - DISPLAY RESULTS:
% ---------------------
L = length(h);
figure,subplot(2,1,1)
stem([0:L-1],h);title('The Prototype Lowpass Filter');
subplot(2,1,2)
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(h,4*L)))));
grid on;

figure
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(1,:),4*L)))));
hold on
for k=2:numbands
    plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(k,:),4*L)))));
end
title('32 CHANNEL FILTERBANK');

figure,subplot(2,1,1)
plot(wav_in);title('input signal')
subplot(2,1,2)
plot(linspace(-1,1,4*N),20*log10(abs(fftshift(fft(wav_in,4*N)))));

figure,subplot(2,1,1)
plot(ys);title('Synthesized Back');
subplot(2,1,2)
plot(linspace(-1,1,4*N),20*log10(abs(fftshift(fft(ys,4*N)))));


mpsk demodulation using digital sample I am trying to receive and demodulate a 8PSK modulated sinusoidal waveform (real not complex) which is constant envelop and changing only in phase of the continuous sinusoid carrier every Symbol Period.
Here symbol 0(000) to 7(111) are mapped to sinusoidal phase of 0:45:360 degree respectively

For further explanation let me set some parameters

Ts = 1 Hz //8PSK Symbol rate
Fc = 8 Hz //center frequency

Now at the receiver I am multiplying with the same carrier and sampling the resultant at lets say Fs = 64 Hz.

I want to find out the transmitted symbol which is actually phase of the received signal from these discrete samples in two cases:

1st assuming receiver carrier phase is in sync with transmitter carrier phase

2nd assuming receiver carrier phase has some fix delta phase difference.

I tried low pass filtering the samples and plotting which does give different levels as per the actual phase modulated but difference between two consecutive phase (45 degree apart) is not linear hence direct thresholding based on average level in symbol period is not working correctly.

Any pointers are welcome.

Not mentioned but from what description you gave it sounds like you have no concern about the Transmitter and Receiver clocks being asynchronous (meaning drifting in phase relative to each other); I assume then that you have confidence in the phase lock of your recovered carrier.

This would not be my preferred approach to do 8-PSK demodulation, but since you are on that path and had specific questions how you could pursue that route, let me offer the following suggestions.
First, when you multiply the two signals (your carrier and the received signal) the result is the sum and difference of the frequencies (and phase) of the two signals. The sum is a 16 Hz signal that you would filter out, and the difference is the cosine of the phase since:

$$cos(\alpha)cos(\beta) = \frac{1}{2}cos(\alpha+\beta)+\frac{1}{2}cos(\alpha-\beta)$$

so
$cos(2\pi f_c t)cos(2\pi f_c t + \phi) = \frac{1}{2}cos(\phi)$ after low pass filtering.

Most significantly this alone has phase ambiguity since it only resolves 0° to 180° uniquely (for instance the output of 0 corresponds to both 90° and 270°).  This can be resolved by using two multipliers; one where you multiply by sine and the other by cosine to get full non-ambiguous 0° to 360° phase demodulation (One output is I and the other is Q and the result is mapping out a circle given by I + jQ.

If you were only interested in phase resolution from 0° to 180°, then I would recommend a saturated multiplication (or specifically and equivalently an X-OR gate) as this would have a linear response of output voltage versus phase after low pass filtering. However in the proper IQ approach above, this would have the deleterious effect of turning the preferred circular response (with all phases then equidistant) into a square.

Probably already obvious but mentioning just in case that it is important to hard limit your signal prior to phase demodulation as the result is also directly sensitive to amplitude (a change in amplitude of the input will lead to a change in amplitude of the output). When you have sufficiently positive SNR, the hard limit operation on the phase modulated signal will actually give you further noise suppression as well (as it eliminates the AM components of the noise).

For other approaches in an all digital receiver, consider using a complex baseband signal and some of the phase detector approaches described in this link High modulation index PSK - carrier recovery.  You are already there with the approach of using two multipliers with a sine and cosine carrier with the addition of a Hilbert Transform to the input signal that is on the 8 Hz carrier.


suppose if we send AM modulated message signal through space how the receiver will be knowing that we applied AM modulation? suppose if we send AM modulated message signal through space how the receiver will be knowing that we applied AM modulation basically how it will be knowing that which modulation is applied for our modulated signal

Communications systems are always designed under the assumption that both emitter and receiver know what "language" they will be speaking to each other.

AM modulation comms are standardized in frequency, channel width for example. So each receiver is materially designed to demodulate AM signals with the appropriate hardware.

When people start trying to find which modulation was applied to a signal its generally involving electronic warfare or reverse engineering.


Cross-correlation of two signals I want to find the direction of a sound source in 2D space. I use two microphones and one sound source for that. I have recorded audio data into two arrays. When a sound is detected, it starts to record data into arrays and stop the recording when no sound is detected. I am confusing about the next step.

How I take the cross-correlation of these two signals using python3?
How I find TDOA using that data?

Your 2-mic array can provide an angle-of-arrival estimate with front-back ambiguity, using cross-correlation to estimate TDOA and then calculating AOA from TDOA and microphone spacing. The front-back ambiguity can be removed with sufficiently directional microphones, or just by constraining the geometry, putting the microphones against walls, for example.
You will also have a multi-path effect. If the space is cluttered, if you lose the direct path, or if the sound source is directional, this multi-path effect will be strong and the errors it produces may be overwhelming.
But even in ideal conditions, the two mics produce a single geo-observable, the TDOA, and thus can calculate only one geometric constraint. In other words, you have only one "known" and thus can solve only one "unknown", not a full 2D position.
If the sound source and the sampling at the microphones can be synchronized, you can also measure propagation delay directly to get distance estimates from each mic - two "knowns". Such absolute distance estimates are also very powerful for resolving the multi-path problem.
If the source and receivers cannot be synchronized, you will need a third microphone to produce a second independent TDOA measurement, a second "known". At this point, you can use the same pseudorange/light-cone cost-minimization calculations that are used in GPS receivers, or just take the two TDOAs and try to intersect the AOA lines.
These are the general principles.
(My comments are based on personal experience with geolocation of radio emitters in several success projects, but the same principles apply to sound.)


How to draw a vector between light and dark spots on an image I'm very new to Python and I need to draw vectors between dark and light spots on an image. Ive managed to convert it to greyscale but I'm struggling with identifying the dark and light spots. If I can identify the co-ordinates of the centre of the bright spots then I should be able to draw a vector on. The aim is to get it to work on an image with one pair of spots then extend this to multiple pairs.

Thanks to anyone who helps.

The keyword here is blob detection. Search OpenCV blob detection and I'm sure there are many examples.

One of the common methods is to use the Laplacian of Gaussian operator.
The Laplacian is a measure of the strength of the second-order derivatives. If you think of the blobs as peaks and valleys on a 3D surface, at the middle of the blobs the slope (derivative) of the surface is changing (second derivative) rapidly from negative to positive and vice-versa.

The Gaussian is there to average out this measurement over a wider area. Any area averaging kernel could be used, but the Gaussian is special* in that increasing the size of the kernel does not introduce any new features. Using multiple kernels creates a scale space, and allows you do scale invariant blob detection.

*I think the Poisson kernel can create a scale-space too with a slight rule relaxation.


Face detection in videos I have to detect faces in videos using machine learning methods like SVM,
I have to build it from scratch not using any inbuilt face detectors.My time limit is like 15hours from now.Can I get any valuable suggestions on how can I begin and proceed?
Thank you

If you only have 15 hours I would built it on a well known
easy to understand algorithm like Viola and Jones, which uses weak binary classifiers and combines them to detect faces. It is relativily fast algorithm that runs in realtime.

You can get more information on the algorithm in the following stack exchange Post


How is the sound converted to matrix in Matlab? I am using the audioread(filename,y,Fs) function of Matlab and it converts an audio file to a matrix. I want to know how is each value of that matrix generated. Is it possible to produce a specific sound like the sound of a guitar by giving your own values to the matrix and then sounding it using the sound() function?

The function audioread doesn't generate any values, it just reads audio samples stored in a file. If you want to generate the sound of a guitar, you need to look into sound synthesis, such as the Karplus-Strong method.


The complexity of zero padded DFTs Assume that I have L points and want to calculate its Fourier transform with N points (N is much larger than L, link N=256, while L=10), I am think if there is way to calculate it with lowest complexity.

At first, I try to calculate it by the definition Fourier transform, according to the definition, the complexity cost should be 2NL.

Then I try to calculate with fast Fourier transform (FFT). As far as I understand, I should make N points input (L points with N-L zeros) to produce N points output. In this case, the complexity should still be the NlogN+N, it is small than 2NL.

Thinking that L is much smaller that N, the result is unaccepted that it takes the same computational cost with N points input. I wonder do you guys have better idea?

The classic complexity measures are much less relevant today than they were in the days of fixed point ALU's taking 40 cycles to do a floating point multiply.

The answer to your question is maybe.  One can decompose an FFT in terms of smaller FFTs,

based on

  Rabiner, Lawrence R., and Bernard Gold. "Theory and application of
  digital signal processing." Englewood Cliffs, NJ, Prentice-Hall, Inc.,
  1975. 777 p. (1975).

there is a section that shows how to do a 1 dimensional dft in two dimensions.

see my answer at

Combine FFT's of shorter length than sample data to get spectrum of all data

and coppying the matlab code from there:

for a 96 point FFT, I formed a 32 X 3 matrix

clear all
M=3;
N=32;
x=linspace(0,10,M*N);
X=reshape(x,N,M).'; % read in as rows
Twiddle=zeros(size(X));
for i=1:M
for k=1:N
Twiddle(i,k)=exp(-1j*2*pi*(i-1)(k-1)/(NM));
end
end
X=fft(X) % fft on each column
X=X.*Twiddle;% element by element product
X=fft(X.').' ; %fft on each row
y=reshape(X,N*M,1); % read out as columns
figure(1)
plot(abs(y),'linewidth',2)
title('Composite DFT')
figure(2)
plot(abs(fft(x)),'linewidth',2)
title('Direct DFT')
figure(3)
plot(x,'linewidth',2)
title('time series')

You need to exploit the zero padding when you form the 2D Matrix. When you do the column FFT operation, simplify by noting the zero elements.

As a simple example, let's say we have 3 data points we want to zero pad for a 16 point DFT. A direct DFT is complexity $16^2$. We need to choose M and N, so that M*N=16. I want all the non zero data in the first row, so M=N=4 , so read in as rows is:
$$
X=\left| \begin{array}{cccc} x[0] & x[1] & x[2] & 0 \\
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
 0 & 0 & 0 & 0 \\
\end{array}\right|
$$
the next step is to do DFTs on the column, normally $ 4 \times 4^2$ complexity. The last column are zeros, and the DFT of a column is zero, so
our complexity is now $ 3\times 4^2$. Im reusing the same memory,
$$
\text{after column DFTs} \; \;X=\left| \begin{array}{cccc}
 x[0] & x[1] & x[2] & 0 \\
 x[0] & x[1] & x[2] & 0 \\
 x[0] & x[1] & x[2] & 0 \\
 x[0] & x[1] & x[2] & 0 \\
\end{array}\right|
$$
There was only one term in each column so I didn't really need to do any multiplies in this step and you don't either.
The Twiddle matrix is usually calculated offline in any FFT , and  forming it is not usually counted in the Complexity calculation.
In most DSP treatments, the top of the array is zero based, but in Matlab is one based.
so in MATLAB
$$ W_{16}^{ik}=
\exp(-\jmath 2\pi\frac{(i-1)(k-1)}{16})
$$
but in most DSP books
$$ W_{16}^{ik}=
\exp(-\jmath 2\pi\frac{(i)(k)}{16})
$$

In this case the product of the Twiddle matrix with $X$ is:
$$
\text{Twiddled X} X=\left| \begin{array}{cccc}
x[0] & x[1]1 & x[2]1 & 0 \\
x[0] & x[1]W_{16}^1 & x[2]W_{16}^2 & 0\\
x[0] & x[1]W_{16}^2 & x[2] W_{16}^4 & 0\\
x[0]& x[1]W_{16}^3 & x[2]W_{16}^6 & 0\\
\end{array}\right| \;
$$
given
$$
\text{Twiddle}\; X=\left| \begin{array}{cccc}
1 & 1 & 1 & 1 \\
1 & W_{16}^1 & W_{16}^2 & W_{16}^3 \\
1 & W_{16}^2 & W_{16}^4 & W_{16}^6\\
1 & W_{16}^3 & W_{16}^6 & W_{16}^{9} \\
\end{array}\right| \;
$$
at this step 9 Complex products

The rest should be obvious (and is left as an exercise)
Use the Matlab code as a template but do the short cuts when you can.


Examples of Causal and Complex Wavelets? I'm looking for wavelets that are both:

causal (a wavelet coefficient $W(j,k)$, where $j$ is scale and $k$ is translation, does not use data from $t > k$ to be estimated).
complex (such as the Morlet wavelet, to analyze phase differences between two time-series).

Do any exist?

You can always just build your own, for instance by taking the Morlet wavelet shift it to the right until, say, 99% of the area below the function is on the positive axis and then set everything on the negative axis to zero. But bear in mind that it's impossible to build a wavelet (or filter which is the same) that is causal and analytic, i.e., a filter that erases all negative frequencies, which you actually need to get that property you need (separation of phase and amplitude). However a pretty good approximation is always possible.


Phase of matched filter for chirp signals I am using the following system: Ultrasonic chirp signals with central frequency 100kHz are being transmitted using ultrasonic sensor. After the transmission through the medium(air/solid/liquid/combination) they are being received with hydrophone. The system is being controled by DSP and signal generation and acquisition is obtained via DAC and ADC. After reception we are performing demodulation and then complex matched filtering. Our aim is to detect the time of arrival between transmission and reception. So far I have been only observing the peak of absolute value of matched filter output as the moment that corresponds to time of arrival, but if I am correct the phase of the matched filter output can also be used.
My question is, should the phase of matched filter output be zero at the moment of the peak in the absolute value of matched filter output?

For a constant frequency signal, with no frequency offset between transmitter and receiver, the phase of a complex correlation, at the correlation peak, is the phase difference between the received signal and the reference signal used for the matched filter.  It need not be 0.  And that is for the actual correlation peak, which may have to be obtained by interpolation between the two samples that surround the peak.
I haven't worked out the math for a chirp.  I would not rely on the phase of the complex correlation being 0 at the correlation peak.


Beginner level confusion regarding symbols used in modulation For 4-QAM symbol sequence, each symbol can be represented by $k=\log_2(4) = 2$ bits. 4-QAM constellation can be generated using s = round(rand(1,T)*2 -1)
For 16-QAM, the symbol sequence takes 4 amplitude levels viz. $\{-3,-1,+1,+3\}$. Each symbol is represented by  $k=\log_2(16) = 4$ bits.
For 64-QAM, the symbol sequence takes 8 amplitude levels viz. $\{-7,-5,-3,-1,+1,+3,+5,+7\}$. Each symbol is represented by  $k=\log_2(64) = 6$ bits. 64-QAM constellation can be generated using : s = randi(8,1,T)*2 -9) +j*randi(8,1,T)*2 -9)) / sqrt(42)

Question: In general, for M-QAM modulation, what is the relation between $M$ and the type of symbols and the number of symbols. For ex: 64-QAM, why do we need 8 amplitude levels and why these symbols?

Please correct me where wrong.

An M-QAM modulation contains M different constellation symbols. Hence, in one symbol you can encode $\mu=\log_2(M)$ bits. (i.e. 6 bits for 64-QAM for example). Now, the bits are equally distributed along the real and the imaginary part of the constellation, i.e. you have $\frac{\mu}{2}$ bits for the real and imaginary part. To encode $\frac{\mu}{2}$ bits, you need $2^{\frac{\mu}{2}}$ different amplitude levels (i.e. 2^3=8 for 64QAM).

Another explanation for the relation of amplitude levels and M: Your M constellation points are arranged in a square. Hence, each side of the square contains $\sqrt{M}$ elements. These are your number of different amplitude levels.

Regarding why the amplitudes have the values of (-7, -5, -3, -1, 1, 3, 5, 7) for 64-QAM. These numbers occur due to two reasons:

the distance between adjacent symbols should be constant (with these numbers the distance is 2). Otherwise, different symbols have different "quality" or error rates. This would be unfair and also does not achieve the optimal distribution in terms of channel capacity.
the constellation should have zero mean (otherwise you would use power for transmitting the mean which does not contain information).

From these two criteria and the addition that you always have an even number of amplitude levels, you actually find that the amplitudes need to be symmetrically distributed around zero. Then, the simplest possibility is to use the values (-7, -5, -3, -1, 1, 3, 5, 7). However, also (-7/2, -5/2, -3/2, -1/2, 1/2, 3/2, 5/2, 7/2) yields a correct 64QAM constellation.


Estimate frequency and peak value of a signal's fundamental How can I find out frequency and peak value of a signal's fundamental?
Please do consider the following:

The peak value of some harmonics might be higher than that one of the fundamental
The frequency of the harmonic is not an integer (e.g. 3.45  Hz)
I have information about where to search for the frequency (3.45 +- 0.5 Hz).
I need to estimate frequency and peak value as accurate as possible

Any response or sample code (MATLAB) is highly appreciated.

There are a few Matlab implementations of various frequency estimators here.  However, they will probably be "tricked" by the higher harmonics.

That you know a range, perhaps the best technique is one of Eric Jacobsen's estimators. Have a look at his page.

In essence, Eric's estimators use the FFT bin at $X(k)$, and on either side of ($X(k+1), X(k-1)$) the peak.  If you choose your FFT length to match your known range (so that $k \rightarrow 3.45$ and $k+1 \rightarrow 3.50$ and $k-1 \rightarrow 3.4$, then you should be good to go.


How to express time delayed discrete signals Hi sorry this question might sound silly but,

a signal $$x[n]=[ 6 \: 4\: 0\: 2]$$ what is $$x[n-2]$$ then?

The problem with representing a discrete time just as a list like: $$x[n]=[6, 4, 0, 2]$$ is that you don't have a time reference.

The first index $n=0$ or the sample in $\textbf{bold}$ or $\underline{\mathrm{underline}}$ is usually taken as the start of the sequence.

For example, $$x[n]=[\mathbf{6}, 4, 0, 2]$$ or $$x[n]=[\underline{6}, 4, 0, 2]$$

So $x[n-2]$ is a delay of two samples, and $x[0]$ now happens two samples later: $$x[n]=[\underline{0}, 0, 6, 4, 0, 2]$$


channel gain calculation I want to calculate the SINR term in an up-link of wireless communication system but there is different models for simulation. I want to know what is the exact meaning of $|h|^2$ in this formula:

$\text{SINR} = \frac{p|h|^2}{\text{Interference+Noise}}$

Is this a Rayleigh fading envelop? If yes, what about path loss? Isn't any difference between near and far users? How can I simulate the channel gains correctly?
Thanks

In your definition of the SINR, $h$ is the channel coefficient. What it means is this: if we send a signal $x(t)$ and model our channel as stationary and frequency-flat then the signal we receive is $y(t) = h \cdot x(t)$ (plus noise) so that its power is $P_y = E\{|y|^2\} = |h|^2 \cdot P_x$. This is why $|h|^2$ appears in your SINR expression: the channel coefficient scales your received signal and therefore $|h|^2$ scales your received power.

Now, in practice, $h$ would contain all the effects that can happen to a signal: (free-space) path loss where $h$ is proportional to $1/d$ where $d$ is the distance (near/far users); small-scale fading, where $h$ varies rapdily so it's best described statistically (e.g., via a Complex Gaussian, which then leads to $|h|$ being Rayleigh distributed); shadowing which can be modeled as a slow variation of the variance of $|h|^2$, e.g., via a log-normal distribution; and more.

Whenever you start modeling $h$ randomly, your received SINR will become a random variable as well. This makes sense though: when you're in a fade, you have little power and as soon as you drop out of it, your power is back to being high. You can then look at statistical moments of the SINR or the quantities derived from it (such as the capacity, leading to things like ergodic and outage capacity, respectively.


Find cut-off frequency of a low-pass filter for a given output signal Suppose we have an ideal low-pass filter $H(e^{j\theta})$ with a cut-off frequency $\theta_c$ in the range of $0\leq \theta_c \leq \pi$.

I want to know the input signal $x[n]$, as well as the cut-off frequency $\theta_c$, which produces the following output signal $y[n]$.

To solve this I started off with transforming $y[n]$ to the frequency domain by calculating the DTFT:

$$
Y(e^{j\theta})=\sum_{n=-\infty}^{\infty}y[n]e^{-j\theta n}=1+e^{-j\theta}+2e^{-j2\theta}+e^{-3j\theta}+e^{-j4\theta}
$$

Now I can calculate $X(e^{j\theta})=\frac{Y(e^{j\theta})}{H(e^{j\theta})}$. This is only possible for $|\theta| \leq \theta_c$, since outside of this region $H(e^{j\theta}) = 0$ and $X(e^{j\theta})$ would approach infinity. But inside this region $H(e^{j\theta})$ is simply 1 and when transforming $X(e^{j\theta})$ back to $x[n]$ it is just the same as $y[n]$, a behaviour that is expected of a low-pass.

The thing I am not sure about is the cut-off frequency. If I simply plot $Y(e^{j\theta})$, I see that there are frequency components until $\pi$. So I could just say alright, $\theta_c=\pi$. But is it really like that? Because the low-pass is actually just dependent on $\theta$ and not all the multiples of it, which I find in the complex exponentials of $Y(e^{j\theta})$. So how do I calculate $\theta_c$?

set $\theta_c = \pi$ and set $x[n] = y[n]$ and you're done.  and you're not dividing by zero anywhere.


Extract graph data from map image I have a map (in image format) describing the power grid of a continent.

Power plants have flags next to them, indicating their type (hydro, gas, oil, etc.). The power lines have a color indicating their tensions (different voltages). The cities are marked as black dots.

Here's a piece of the map, courtesy of the ENTSO-E consortium

Full maps here. Things get much more complex in populated area.

As you can see, the power lines are often overlapping each other. Moreover they are not simple straight lines. Two cities (points) can be connected by a broken line, with no middle point in between.

I need to turn that map image into a graph (dataset). The broken lines may be straightened and have their full length marked as the weight of the arc. Or they may be handled by inserting middle points.

I'd like to automatize as much work as I can, but, since the map is so complex, I understand I'll probably need some manual tool to do the final pass.

What would you suggest?

I think finding the power-plants should be possible using an automatic tool, then corrections and addition of arcs could be done using an assisted manual tool.

EDIT: perhaps this is something that can be done with a plugin for GIS software? If programming is involved, my favorite languages are Java and Python. I know C++ too.

Start by getting rid of the background (Mathematica example):

  map = Import["http://i.stack.imgur.com/HYtbo.jpg"]

  fg = ImageAdjust[map, {1.0, 1.0}]

Then use a graph extraction tool like NEFI. If you want to do it all in Mathematica, that's possible too. You might want to add another preprocessing step to decompose the graph by edge type if you want to treat them separately, like you said.


Career in signal processing I am doing my MSc in electrical engineering, and there are not many courses in my university on signal processing. I want to develop my skills in signal processing such that I have a hold in one particular domain.

If anybody could help me and guide me on how to start developing skills it would be helpful.

Skills like which languages to be known. What are the basic things to be familiar with in basic signal processing and image processing? If there is any other domain of signal processing where signal processing engineers are handy, please let me know. Finally what books are helpful in developing the skills in DSP?

To be honest, most people end up doing something else other than what they directly study.

The career follows the job you get.

Get a good background in your departments core strengths. you really can’t anticipate what skills you will eventually develop with precision.  Soft skills are actually important.


Moving Average FIR filter coefficients I am trying to use a FIR filter making software to create a moving average filter and I'm a little confused on what I should put in the coefficients part? I read that the coefficients are always 1 for moving average so this is what I did - is that correct? (Picture below)
Other than the coefficients, I just wanted to double check (because I am very new to working with filters):
Moving average filters are symmetrical, have no L-th band filter, and are single rate filters right?

I read that the coefficients are always 1 for moving average so this is what I did - is that correct?

Technically no. For a moving average filter of length $N$ all the coefficients are $1/N$. In practice most people use $1$ for all coefficients and then scale the output with $1/N$ It gives the same result and it's typically more efficient, but if you were to look at the difference equation, the coefficients would show up as $1/N$

Moving average filters are symmetrical,

Depends a bit exactly how exactly you define symmetry. You can certainly time reverse it without changing anything.

have no L-th band filter,

I have no idea what you are mean by "L-th band filter"

and are single rate filters right?

The concepts of multi-rate and moving average are mostly orthogonal. You can use the together , individually or not all.


Duration of unknown rectangular pulse with additive white Gaussian noise Problem.

There is a discrete signal $f[i]$ (example below).

It is known, that $f[i]$ have a form of rectangular pulse with additive white Gaussian noise.

$f[i] = s[i] + n[i]$,
$s[i] = \alpha(\theta[i - i_{1}] - \theta[i - i_{2}]) + c$,
$i_{2} > i_{1}$,
$i_{1} > N$

Where
$\theta[i]$ is a Heaviside step function,
$n[i]$ is an additive white Gaussian noise,
$\alpha$ is a height of rectangular pulse,
$i_{1}$ is an index of the first sample of rectangular pulse,
$i_{2}$ is an index of the last sample of rectangular pulse,
$c$ is constant level of signal,
$N$ is adjustable parameter.

All parameters may have large range of values.
It is required to find value of $(i_{2} - i_{1})$ (duration of rectangular pulse in samples).

Possible solutions.

At the moment, I have tried two ways to solve this problem.

Low pass filter with threshold.

As first attempt, I have used simple scheme with low pass filter and threshold.
1. Apply FIR low pass filter with cutoff frequency equal to $0.05f_{sampling}$.
2. Estimate mean $m$ and dispersion $\sigma^{2}$ of filtered noise from first $N$ samples of the signal.
3. Set threshold $t = m + 3\sigma$.
4. Estimate $i_{1} = \min_{i}(f[i] > t)$.
5. Estimate $i_{2} = \max_{i}(f[i] > t)$.

Pros:
1. This algorithm is simple.
2. It is easy to write fast implementation.

Cons:
1. It is hard to estimate efficient value of filter's cutoff frequency. On the one hand low value can corrupt the form of short pulses. On the other hand large value decreases effect from filtration.
2. Algorithm isn't using all information, we have about the signal.

Regression analysis

As the second attempt, I have tried to approximate input sequence of samples with the function $s'[i]$.

$s'[i] = \alpha(\theta'[i - i_{1}] - \theta'[i - i_{2}]) + c$,
$\theta'[i] = \frac{1}{1+e^{-\frac{i}{\beta}}}$, where $\beta$ is a small parameter.

For approximation I have used the method of least squares with gradient descent for minimizing cost function.
1. Set initial values for $\alpha$, $c$, $i_{1}$, $i_{2}$.
2. Perform gradient descent.
3. Set threshold $t = c + 0.5\alpha$.
4. Estimate $i_{1} = \min_{i}(f[i] > t)$.
5. Estimate $i_{2} = \max_{i}(f[i] > t)$.

Pros:
1. This algorithm gives results with good precision.
2. It works for wide range of durations.

Cons:
1. It is very slow.

Question.

After all, I am not satisfied with the precision of the first algorithm and with the speed of the second one.
How would you solve this problem?
Is there any classical solution, that I failed to find?
Ideas, links, any feedback will be much appreciated.
Thank you.

You want a method that removes noise while preserving edges. This cannot be achieved well by linear filtering, as you noticed yourself. I know of two approaches that might work well for your problem. The first is median filtering, where samples inside a window are replaced by their median. The following plot shows the result of median filtering with a window length of 25 samples (in red):

The other, more complex, approach is Total variation denoising, which works very well for piecewise constant signals. There is a very good description of total variation denoising including Matlab code available: link.


Derive chernoff bound for $\mathrm{erfc}(x)$ I'm trying to derive the Chernoff bound $\mathrm{erfc}(x) \le \exp(-x^2)$, by first showing:

$$\mathrm{erfc}(x) = \frac{2}{\pi}\int_{0}^{\frac{\pi}{2}}\exp\left(-\frac{x^2}{\cos^2\theta}\right)d\theta$$

This equality should be derived by considering 2 independent, standard gaussian random variables $x_1, x_2$ and the region in the $x_1x_2$-plane where $|x_1|\le x$.

What I have so far, is come up with a diagram like this:

where I consider the the probability where $x_1, x_2$ falls in the square (4 $\times$ "blue square"). Then

\begin{align}
P(x_1<x,x_2<x) &= P(x_1<x)P(x_2<x)\\
&=\frac{1}{2\pi}\int_{0}^{x}\int_{0}^{x}e^{-\frac{x_1^2+x_2^2}{2}}dx_1dx_2\\
&=\frac{1}{2\pi}\int_{0}^{\frac{\pi}{2}}\int_{0}^{\frac{x}{\cos\theta}}e^{-\frac{r^2}{2}}rdrd\theta\\
&=\frac{1}{2\pi}\int_{0}^{\frac{\pi}{2}}\int_{0}^{\frac{x^2}{2\cos^2\theta}}e^{-v}dvd\theta\\
&=\frac{1}{2\pi}\int_{0}^{\frac{\pi}{2}}-e^{-v}\Big|_0^{\frac{x^2}{2\cos^2\theta}}d\theta\\
&=\frac{1}{2\pi}\int_{0}^{\frac{\pi}{2}}1-e^{-\frac{x^2}{2\cos^2\theta}}d\theta\\
&=\frac{1}{4} - \frac{1}{2\pi}\int_{0}^{\frac{\pi}{2}}e^{-\frac{x^2}{2\cos^2\theta}}d\theta
\end{align}
This is the probability for the blue square, so for the entire square, I have:
\begin{align}
4P(x_1<x,x_2<x)&=1-\frac{2}{\pi}\int_{0}^{\frac{\pi}{2}}e^{-\frac{x^2}{2\cos^2\theta}}d\theta
\end{align}
The probability of falling outside the square is:
$$1-\Big[1-\frac{2}{\pi}\int_{0}^{\frac{\pi}{2}}e^{-\frac{x^2}{2\cos^2\theta}}d\theta\Big] = \frac{2}{\pi}\int_{0}^{\frac{\pi}{2}}e^{-\frac{x^2}{2\cos^2\theta}}d\theta$$

I'm kind of stuck here. I think this is my $[Q(x)]^2$, but according to this paper (Eq 9), my $Q(x)$ is incorrect. How can I derive the $Q$-function in the paper using this gaussian RV derivation method?

ETA: As per Dilip Sarwate's suggestion, I considered the circle that circumscribes the squares, and I got the following by considering just 1 (upper-right) quadrant.

\begin{align}
P(\text{in the circle quadrant}) &= \frac{1}{2\pi}\int_{0}^{\sqrt{2}x}\int_{0}^{\sqrt{2x^2-x_2^2}}e^{-\frac{x_1^2+x_2^2}{2}}dx_1dx_2\\
&=\frac{1}{2\pi}\int_{0}^{\frac{\pi}{2}}\int_{0}^{\sqrt{2}x}e^{-\frac{r^2}{2}}rdrd\theta\\
&=\frac{1}{2\pi}\int_{0}^{\frac{\pi}{2}}\int_{0}^{x^2}e^{-v}dvd\theta\\
&=\frac{1}{2\pi}\int_{0}^{\frac{\pi}{2}}-e^{-v}\Big|_0^{x^2}d\theta\\
&=\frac{1}{2\pi}\int_{0}^{\frac{\pi}{2}}1-e^{-x^2}d\theta\\
&=\frac{1}{4}-\frac{1}{4}e^{-x^2}
\end{align}

\begin{align}
\text{P(in the blue square)} &< \text{P(in the circle quadrant)}\\
\frac{1}{4} - \frac{1}{2\pi}\int_{0}^{\frac{\pi}{2}}e^{-\frac{x^2}{2\cos^2\theta}}d\theta &< \frac{1}{4}-\frac{1}{4}e^{-x^2}\\
\frac{2}{\pi}\int_{0}^{\frac{\pi}{2}}e^{-\frac{x^2}{2\cos^2\theta}}d\theta &> e^{-x^2}
\end{align}
Let $v = \frac{x}{\sqrt{2}}$

\begin{align}
\frac{2}{\pi}\int_{0}^{\frac{\pi}{2}}e^{-\frac{v^2}{cos^2\theta}}d\theta &> e^{-2v^2}
\end{align}
I still don't quite understand how my lefthand expression became erfc(v), since it was calculated from $P(x_1<x,x_2<x)$, and how to eventually get $\mathrm{erfc}(x) \le \exp(-x^2)$.

I don't know why you are approaching the problem the way you are by going to a two-dimensional approach.

The Chernoff bound starts with the observation that for any random variable
$X$ and any $\lambda > 0$,
$$P\{X \geq a\} \leq E[\exp(\lambda(X-a))]\tag{1}$$
because $\exp(\lambda(x-a))\geq  \mathbf 1_{\{x\colon x \geq a\}}$.
Consequently,
$$P\{X \geq a\} \leq \min_{\lambda > 0}E[\exp(\lambda(X-a))].\tag{2}$$
For a standard Gaussian random variable, $E[\exp(\lambda(X-a))]$
can be found easily by using a method called "completing the square"
(in the exponent which as terms like $-\frac{1}{2}x^2+\lambda(x-a)$ in
it) and noting that the resulting integral is the integral
of a Gaussian density times a constant. Alternatively, you can
use the moment-generating function and read off the answer.
After minimization with respect to $\lambda$, the result is
$$P\{X \geq a\} = Q(x) \leq e^{-\frac{x^2}{2}}.\tag{3}$$
The two-dimensional approach leads to the slightly better result
$$Q(x) \leq \frac 12e^{-\frac{x^2}{2}}.\tag{4}$$
Note that you have proved already that for $x > 0$,
$$P\{X_1>x, X_2>x\} = Q^2(x) \leq \frac 14 e^{-{x^2}}$$
from which $(4)$ follows immediately, and $(4)$ is not the
Chernoff bound on $Q(x)$.


Histogram estimation (from sub-set of image pixels) Problem: How to find gray-scale image histogram without iterating all pixels and binning them.

Description:
Suppose that I have large set of images with no common scene and of different resolutions. For each image I would like to compute histogram but to save processing time, I want just an estimate rather than finding the full histogram.

One way which is more or less obvious is to take just sub-set of pixels in a clever way and use these to determine the estimate. The problem is then how to sample the image (which pixels to take) so that the estimate is really representative (similar to full histogram). What would be a good sampling technique to use? E.g. the simplest may be to take N samples at random positions. I have thought of some approaches but I'd be interested if there are any that were proven to give good estimate or, if there are none, in your ideas what could work.

There may be also different approaches to histogram estimation that I didn't think of other than (or complementing) taking sub-set of pixels. I'd be glad for pointing out anything that may improve estimate or processing speed.

I was searching for some time but there seems to be no scientific papers or other resources that would discuss this issue.

I think random sampling approach seems to be not effective, since the statistical population (pixel intensities) distribution in images is heavily localized. There might be more scientific approaches, but if you need only a rough near estimation I'd like to suggest a method. Here is how:

Blur the image with a low-pass mask (average out neighboring pixels) [Optional, if you are O.K with blurring computations]
Down sample (make it smaller like by a factor of 2) the blurred image
Compute histogram of the down-samples blurred version of the image.

Here is a code experiment:

I=imread('coins.png');
Mask=ones(5,5);
tic
H1=imhist(I);
toc
I_Averaged=imfilter(I,Mask); %OPTIONAL STEP
I_DownSampled=I(1:3:end,1:3:end);
tic
H2=imhist(I_DownSampled);
toc

The results:

Elapsed time is 0.003042 seconds. %For Original Hist
Elapsed time is 0.000365 seconds. %For DownSampled Hist

and the rough histogram and original image histogram are shown in below images:

1.Rough

  2.Original

Almost no difference, for 8 times faster computations.


What is normalized frequency I am working on DSP and finding difficulty to understand the term Normalized frequency often used with DFT & DTFT.

What is normalized frequency in DSP ? and how it is different from analog frequency ?

What is significance to normalize frequency in DSP ?

Why limit of normalized frequency is 2π ?

How FFT deals with normalized frequency ?

Normalized frequency is frequency in units of cycles/sample or radians/sample commonly used as the frequency axis for the representation of digital signals.
When the units are cycles/sample, the sampling rate is 1 (1 cycle per sample) and the unique digital signal in the first Nyquist zone resides from a sampling rate of -0.5 to +0.5 cycles per sample. This is the frequency equivalent of representing the time axis in units of samples instead of an actual time interval such as seconds.
When the units are radians/sample, the sampling rate is $2\pi$ ($2\pi$ radians per sample) and the unique digital signal in the first Nyquist zone resides from a sampling rate of $-\pi$ to $+\pi$.
How this comes about can be seen from the following expressions:
For an analog signal given as $$x(t)=\sin(2\pi F t)$$ where F is the analog frequency units in Hz,
When sampled at a sampling frequency of $F_s$ Hz, the sampling interval is $T_s=1/F_s$ so the signal after being sampled is given as:
$$x(nT_s)=\sin(2\pi F nT_s) = \sin\left(\frac{2\pi F}{F_s}n\right)$$
Where the units of normalized frequency, either $\frac{F}{F_s}$ in cycles/sample or $\frac{2\pi F}{F_s}$ in radians/sample is clearly shown.
This is illustrated below using $\Omega = 2\pi F$
Update: As @Fat32 points out in the comments, the units for sampling rate $F_s$ in the figure below is "samples/sec" (which is "Hz") in order for the normalized frequency to become radians/sample.

To visually see the concept of "radians/sample" (and most other DSP concepts dealing with frequency and time) it has helped me considerably to get away from viewing individual frequency tones as sines and/or cosines and instead view them as spinning phasors ($e^{j\omega t} = 1 \angle (\omega t)$) as depicted in the graphic below, which shows a complex phasor spinning at a rate of 2 Hz and it's associated cosine and sine (being the real and imaginary axis). Each point in a DFT is an individual frequency tone represented as a single rotating phasor in time. Such a tone in an analog system would continuously rotate (counter-clockwise if a positive frequency and clockwise if a negative frequency) at F rotations per second where F is the frequency in Hz, or cycles/second. Once sampled, the rotation will be at the same rate but will be in discrete samples where each sample is a constant angle in radians, and thus the frequency can be quantified as radians/sample representing the rate of rotation of the phasor.