How to understand the relationship between the filter time constant and the half-power cutoff frequency I learned that the time constant can be computed as $\frac{1}{2 \pi f_0}$, where $f_0$ is the half-power cutoff frequency of a high-pass filter.

However, I was wondering how the time constant and the half-power cutoff frequency can be related. For example, if the half-power cutoff frequency is 2 Hz, why must the time constant be $\frac{1}{2 \pi (2 \text{ Hz})}$ seconds? Why can't it be some arbitrary number, as long as at that time the output reaches $\frac{1}{e}$ of the starting value?

My knowledge of signal processing is limited to a very basic understanding of Fourier analysis and filtering, so a bit more illustration would be appreciated.

If you have a first order high-pass filter with transfer function

$$H(s)=\frac{s}{s+\omega_0}\tag{1}$$

then the squared magnitude of its frequency response is given by

$$|H(j\omega)|^2=\frac{\omega^2}{\omega^2+\omega_0^2}\tag{2}$$

Consequently, for $\omega=\omega_0$ we have

$$|H(j\omega_0)|^2=\frac12\tag{3}$$

and $\omega_0$ clearly is the $3$-dB (half-power) cut-off frequency.

In the time domain, we could look at the step response of that system. The step response corresponding to the transfer function $(1)$ is

$$s(t)=e^{-t/T}u(t),\quad T=\frac{1}{\omega_0}\tag{4}$$

where $u(t)$ is the unit step function. So the time constant $T$ and the $3$-dB cut-off frequency are inverses of each other. This is what you get from the math (Fourier transform).

I'll try to give you some intuition and insight into this result. In general a high-pass does not pass DC (i.e., a constant value), so after the step has been applied at the input, the output must decay to zero because the value of the step function  is constant after the step has occurred. The higher the cut-off frequency $\omega_0$, the faster the output must decay to zero, because the filter suppresses low frequencies better than a filter with a lower cut-off frequency. So a high value for $\omega_0$ implies a smaller time constant $T$. You can also look at it from the other side. Assume an extremely low cut-off frequency $\omega_0\approx 0$. This means that the filter hardly suppresses anything, i.e., it leaves the input almost unchanged. So at the output you will observe almost the same as at the input, i.e., a step. Only after a very long time will the output decay to zero (because DC is suppressed, after all). So a small value for $\omega_0$ corresponds to a large value of $T$.


How to gather similar objects in an image? I am working in a traffic tracking app, I am using a background subtraction approach to obtain information about the moving vehicles. After some morphological operations I obtain blobs of the moving objects. In some cases more than one blob belong to the same object, e.g.

How can I gather the isolated blobs to track only a single object? what features can allow me the fusion of these blobs in a single object? Thanks

First, I personally wouldn't recommend you using background subtraction in general for tracking applications.
There are way better methods like L1-minimzation of lasso problems, Kenalized Correlation Filters, Multiple Instance learnng, Struck, and whole bunch of stuff.

But if you must use background subtraction,, here are my thoughts:

You could obviously at the first frame apply some very basic background subtraction algorithm like using median filtering.
Now once you get these blobs, you can start generating a dictionary out of each blob.
That means you can now treat each blob as an individually part where you learn a dictionary for each. Usually in dictionary learning from first frame, people usually perturb the bounding box covering the blob a little (slight translations)(Maybe by gaussian sampled translations with zero mean) and then collect all these samples and vectorize them into your dictionary matrix.
You do this for all the other parts, or blobs if you will.

Now at the first frame you already have the bounding box governing the object ,which is very typical in any tracking problem. Also, you have the bounding box governing all these blobs individually.

What you do next in the next frame is to try track each of these parts independently say using simple cross correlation.
At last, you only need to learn one affine transformation that transforms the big box from the previous frame to cover all the blobs in the next frame.

In this approach you divided your main task of tracking from 1 bounding box to smaller ones on which you can use simple methods of tracking. Then you combine 1 bounding box by either learning an affine transformation or simply plotting a minimum bounding box that covers all the other parts.


What is the signal-noise model of laser induced fluorescence? The following image is an example of a fluorescence signal trace from the photomultiplier tube detector (PMT). Each spike represents one fluorescently labeled cell.

This signal model is definitely not additive white Gaussian, although the noise term is Gaussian. The signals are low variance but the point is that when the signal amplitude increases, the noise variance rapidly decreases. In the given example, above the intensity value of $0.6$ all three Gaussian signals are almost noise free. Around $0.4-0.5$ intensity of the noise is reduced but can be seen if zoomed. Below $0.4$ the noise is significantly apparent.

  How can this signal-noise relation be properly modeled? For example for additive noise we have $$x[n]=s[n]+w[n]$$ where $x$ is the received signal $s$ is the error free signal and $w$ is the Gaussian noise.

My aim is to simulate such signals in Matlab using a theoretical model. However, I was not able to see and/or find any model via google search or on some related books.

The physics behind this is due to detection of laser signal via photodiodes or PMTs. I was thinking that the noise source was probably the quantum noise but I am not sure because for the quantum noise, the noise variance does not decrease when the total number of photons $N$ increases. Instead it increases as $\sqrt N$.

What I tried was dividing the noise values by $k+s[n]^2$ before adding to the signal, where $k$ is some constant. I did this over additive Gaussian noise model. The results were okay. But I am not able to justify this theoretically. It is just my intuitive understanding.

Sources:

https://www.cut.ac.cy/faculties/fet/mem/research/research-units-and-laboratories/biolisys/research/imaging/cytometry/?languageId=1

https://en.wikipedia.org/wiki/Laser-induced_fluorescence

Second try!  This question is also at the Electrical Engineering stack, so I posted a short 'starter' answer there. A mod converted it to a comment and my tiny rep there prevented me making any further comments. This is what I posted there:

I suggest this link to cbeleites's answer to a question involving laser-excited Raman spectroscopy: https://stats.stackexchange.com/a/194247/247352 and this golden oldie paper: N.J. Dovichi, J.C. Martin, J.H. Jett, M. Trkula, R.A. Keller, "Laser-Induced Fluorescence of Flowing Samples as an Approach to Single-Molecule Detection in Liquids", Anal. Chem.,56 (1984) 348-354. There are thousands of laser-induced (aka, laser excited) fluorescence papers, but newer research papers tend not to dwell on the fundamentals.

My research career, prior to retirement, involved the use and study of laser-based spectrometries, including laser-induced fluorescence (LIF) of molecular and atomic species. So the reason for the above paragraph was that I simply cannot understand how the noise can decrease with increased number of photons. Likewise, the correlation of the signal and noise is puzzling. Hence the references above.

I have been trying to simulate what I think is happening and get this:

Maybe I am on the wrong track, but my model used 1) shot noise on the fluorescence peaks, 2) additive, white gaussian noise on the baseline and 3) some laser intensity fluctuations. It would help to know exactly how the PMT (or photodetector) output current leads to the observed fluorescence signals. Also, are there any unusual time constants that could result in unsuspected pulse stretching? The peaks in the OP's figure presumably have widths determined by the times it takes labelled cells to transit the region where the laser is focussed. This would be on the order of ms, so perhaps provides a clue about the correlation issue.

EDIT in response to OP queries. Without knowing any details, I assumed a fairly typical system: 1) CW laser excitation, 2) fluorophore-labeled cells, 3) a flow cytometer of some sort, 4) microscope objective, 5) PMT or other high gain photodetector, 6) preamplifier, 7) data acquisition, and so forth. As noted above, I assumed shot noise on the fluorescence emissions, white Gaussian noise (from the preamplifier) on the baseline, and about 0.5% Gaussian intensity fluctuations on the excitation laser. I ignored the variation in cell fluorophore loadings, because this has somewhat similar effect as laser intensity fluctuations. The simulation model (program, actually) is shown below:

Every simulation run (start = 0, end = 0.6, $\Delta t = 0.001$) results in a new temporal response plot. The numbers I used were just arrived at by trial and error, trying to get close to the OP's plot. Lowering the gains on the three peaks, via the triangular gain blocks that have 'x 1000' under them, results in lower SNRs. The next two plots show the results for gain = 300 and gain = 100.

I do not know if this is going to be of any help to the OP. Maybe it is over-simplified or I am missing something essential, but I did not have much to go on. I published my first laser induced (excited) molecular fluorescence paper in 1981 and my last in 2014, but I never encountered a molecular LIF system where the noise decreased with increasing number of emitted fluorescence photons. The correlation is also puzzling to me since I cannot think of how it might arise. Obviously, I am very curious to know what is going on!


What are the advantages, if any, of derivative sampling? In Five short stories about the cardinal series $[1]$, the author makes the following comment:

  Interestingly enough, Shannon goes on to mention that other sets of
  data can also be used to determine the band-limited signal--for
  example, the values of ƒ and its first derivative at every other
  sample point, the values of ƒ and its first and second derivatives at
  every third sample point, and so on.

The paper mentions some historical developments, but I'm curious what the "killer apps" are for derivative sampling.  Does it go by any other names?  Are there further generalizations of this approach?

A simple overview, or pointers to some references would be great.

--

J. R. Higgins, Five short stories about the cardinal series, Bull. Amer. Math. Soc. (N. S.) 12 (1985), no. 1, 45-89. http://bit.ly/plioNg

Papoulis introduced a generalization of the sampling theorem [1], of which derivative sampling approach is one case. The gist of the theorem, quoting from [2] is:

  In 1977, Papoulis introduced a powerful extension of Shannon’s sampling theory, showing that a band-limited signal could be reconstructed exactly from the samples of the response of $m$ linear shift-invariant systems sampled at $1/m$ the reconstruction rate.

Perhaps one reason why it's hard to search for the term is because Papoulis' generalized sampling theorem is mentioned more often than "derivative sampling". [2] is also a very good article which presents a broad overview of the sampling approaches at the time of publication. [3], also by the same author is an extension of [1] to the class of non-bandlimited functions.

As for applications, in a recent paper [4], the derivative sampling approach is used to design wideband fractional delay filters and the authors show that sampling the derivative results in smaller errors. From the abstract:

  In this paper, the design of wideband fractional delay filter is investigated. First, the reconstruction formula of derivative sampling method is applied to design wideband fractional delay filter by using index substitution and window method. ... Finally, numerical examples are demonstrated to show that the proposed method has smaller design error than the conventional fractional delay filter without sampling the derivative of signal.

While there certainly are more, I'll refrain from posting more references and application to keep it short (and avoid it turning into a list). A good point to start looking would be to check which papers have cited [1]-[3] and narrow down the list based on the abstract.

[1]: A. Papoulis, “Generalized sampling expansion,” IEEE Trans. Circuits and Systems, vol. 24, no. 11, pp. 652-654, 1977.

[2]: M. Unser, "Sampling - 50 years after Shannon," Proceedings of the IEEE, vol. 88, num. 4, p. 569-587, 2000

[3]: M. Unser and J. Zerubia, "A generalized sampling theory without band-limiting constraints," IEEE Trans. Circuits and Systems II, vol. 45, num. 8, p. 959–969, 1998

[4]: C-C Tseng and S-L Lee, "Design of Wideband Fractional Delay Filters Using Derivative Sampling Method", IEEE Trans. Circuits and Systems I, vol. 57, num. 8, p. 2087-2098, 2010


Is sparse signal in time domain is sparse in Frequency domain In OFDM system, if we have a multi-path environment, we supposed to have a sparse received signal in time domain, then  $FFT$ is performed to convert it into frequency domain.

My question, does the results of FFT in frequency domain must be sparse since its equivalent time signal is sparse too ?

thank you.

Your central question first:

  My question, does the results of FFT in frequency domain must be sparse since its equivalent time signal is sparse too ?

No. In fact, it does the opposite: signals sparse in one domain are typically dense in the other.

If you want so, that's happening for the same reason that we have Heisenberg's Uncertainty Principle: Operators like the DFT don't allow us to preserve correlation / concentration of energy when going from one domain to the other (or back).

Let me comment on your original statement:

  In OFDM system, if we have a multi-path environment, we supposed to have a sparse received signal in time domain, then FFT is performed to convert it into frequency domain.

That is false. The characterizing property of multipath environments is that they convolve the transmit signal with a nontrivial impulse response.

That converts even very sparse time-domain signals (e.g. a single dirac impulse) to dense ones. The only case where a channel would introduce sparsity would be if the the impulse response of the channel happened to be a matched filter to the transmit waveform.

But that can not be the case for an OFDM system, as that would directly contradict the principle of OFDM which is that you transport different bits of information on different frequency bins.

You might be considering a predistorting system, but then you'd mention that, and also, you'd know about the spectral properties of your signal (and hence could've negatively answered your original question yourself).

So, it feels like you've either taking something completely out of context, or you're misunderstanding something fundamentally.


What are some promising areas of research in digital signal processing? I intend to pursue a PhD program in digital signal processing (DSP). However, it seems nowadays that research has saturated in DSP and it's all about machine learning and deep learning. Are there are any areas in pure DSP that can still be explored for a PhD? If so, what are some research groups that focus on these areas?

Choosing a novel topic is a bold and wise move. Your question requires a bit of a bet, yet it is interesting to me, as I have not thought about this for a while. I will first drop two dual thoughts (they could be updated):

shed light on black-boxes: the present successes of machine learning and deep networks shall not make us forger that signals are our (DSP community) sensing door to a globally unknown world. Truly understanding what we can understand with data and models remains at our core, and we can be much better at that (to make it more interpretable, robust, model-based, explainable, reliable etc.)
save electrons in devices and computations: sensing and computing are becoming pervasive (IoT, etc.), and this requires materials (from the extraction of rare-earth elements for sensors and computing devices) and resources (form energy-harvesting sensing to bit-wise or qubit-related fast algorithms). Saving power (efficient filters, compression, fast algorithms, etc.) are in our DNA.

DSP is a blended field of knowledge, where people of mixed origins (computing, acquisition, information theory, physics, optimization, computer vision, statistics) meet. I hope that we DSP people can innovate here.
Additional stuff, as they come:

The Future of Signal Processing: conference program (2017)


How to statistically quantify the coverage of features in the data? I have an image. I find features on it using different methods. I want to quantify that which method performs the best. The creteria is how many features per unit area of the image are and how they are spread. More the features are spread, the better it is. How to visualize them.

So at the end the criteria is "how many features per unit area are found and how well are they spread all over the image"

It will be great if there are some suggestions that how to visualize the quantified results.

Thanks a lot.

To visualize the features, why don't you just do scatter plots? I assume that if you do feature extraction, you know where in the image the features are found. If different features have different "sizes", i.e. some features are point features, whereas others take up e.g. a 5x5 window (and are non-overlapping), you may want to adjust the size of your symbols for the scatter plot to visually improve the relative amount of coverage.

If you want to quantify the coverage, then you're looking at two independent measures

Density: To this end, you can simply divide the number of features by the image size.
Homogeneity: I've seen Ripley's K used to great success when looking at whether data is homogeneously distributed, or whether it is clustered.


Good Reference Problem to Test Filtering/Estimation Algorithms I am looking to figure out if a current filter algorithm I have built could be useful for some problems I am looking into at work. It isn't a Kalman filter, but is instead making estimations using a Neural Network and the latest $N$ time series measurements for some state.

I have tested it on a 2D particle tracking problem I put together with what I believe to be fairly large Gaussian noise, and it has performed decent as far as I can see. An example figure is below:

However, I don't know if this sample problem is challenging enough to prove this algorithm works sufficiently well and I also don't have any reference algorithm results to compare against.

Is there any sample problems seen in the literature or elsewhere that I could run to test this filtering algorithm against so I have a reference of what's good or not?

This answer might give you one example problem and a KF implementation to solve it.

The usual problems are kinematics problems, where you want to track position, velocity, and acceleration in one or more dimensions.

There are some nonlinear problems, too, like frequency / phase tracking but that might be overkill. See how it works on linear problems first.


2-sided Regions of convergence for Z transforms Given a z transform with one pole can you have a 2 sided Region of convergence or does 1 pole limit it to being only left or right sided? I know when you have two poles the 2 sided scenario is when a "ring" forms between the two poles in the z-space but when you have one pole does the case when the entire z space is in the RoC excluding the ring formed by the one pole count as the two sided case?

With one pole you have only two possible regions of convergence, either $|z|<r$  (where $r$ is the pole radius), which corresponds to a left-sided sequence, or the region $|z|>r$, corresponding to a right-sided sequence. I think the misunderstanding lies in the concept of 'the ring formed by the one pole', because there is no such ring.


Signal Reconstruction after fourier transform I'm working from an example posted here.  I understand the steps to acquire the fourier transform and can clearly see the spikes at normalized frequencies at 15 and 40 Hz from the 0-centered periodogram.  Knowing this, I believe that I can reconstruct a smoother version of the signal as:

$x_{\text{reconstructed}}(t)=\alpha_1 sin(30\pi t)+\alpha_2 sin(80\pi t)$.

I have two questions related to this reconstruction:

How do I obtain the coefficients $\alpha_1$ and $\alpha_2$ without an inverse fourier transform of the entire frequency domain data set?  Is there a more efficient way?
How could I have obtained the 15 & 40 Hz frequencies from the transformed data?  I know I can sort the transformed data to determine that these two frequencies had the highest two powers.  But if the data set were very large, this might be unfeasible.  Is there another way to determine the important frequencies?

Suppose you want to perform N points FFT, then the evenly spaced frequency vector is given by

  f= (0:NumUniquePts-1)*Fs/N;

where NumUniquePts = ceil((N+1)/2) is the number of unique points in f, and Fs is the sampling rate.

So if fftx=fft(x,N), then $\alpha_1$ = fftx[15/(Fs/N)+1] and $\alpha_2$ = fftx[40/(Fs/N)+1] (provided 15 and 40 can be divided by Fs/N).

See here for more.


Mathematical description of the ACF using fft2 I computed the acf of an image with the following code:

%# read in image
Img = imread('myImage.tif');

% get size of image
[N, M] = size(Img);

%# convert to double
I = double(Img);

%# subtract mean
I = I-mean(I(:));

%# normalize magnitude
I = I/sqrt(sum(I(:).^2));

%# compute acf
fft_I = fft2(I);
Acf = real(fftshift(ifft2(fft_I.*conj(fft_I))));

Now, i have to mathematically describe the procedure, and I am thinking that the FFT equation given in https://uk.mathworks.com/help/matlab/ref/fft2.html (scroll to the bottom) isn't enough?

Can someone please help me out here. Many thanks.

EDIT:
Please find the added mathematical description i have come up with and kindly help me verify its correctness.

When the signal is assumed to be ergodic, then its ACF can be computed using time averages which can also be computed using the following convolution of length $N$ sequence $x[n]$ by its conjugate symmetric version:

$$ \hat{r}_{xx}[m] = \frac{1}{N} ~~~x[m] ~~\star ~~ x^*[-m] ~~ $$

Taking the $2N-1$ point DFT of both sides yields the following:

$$ \hat{R}_{xx}[k] = \frac{1}{N} ~~~X[k] X^*[k] =  \frac{1}{N} |X[k]|^2  $$

In a MATLAB implementation, an fftshift is required to bring $r_{xx}[0]$ to the center position for convenience. Hence the ACF will be obtained by an $2N-1$ point inverse DFT/FFT as

rxk = (1/N)*real(fftshift( ifft( abs(fft(x,2N-1)).^2 ,2N-1) ) );

The 2D case follows very similarly to this 1D outline...


Overlap-windowing vs preserved information When windowing non-overlapped segments of data stream, part of the data is being "ignored" at the boundaries of the segments due to the tapering of window functions (non-rectangular). Some information is lost in further processing.

My question is: for a given window function, does overlapping the segments by $50 \%$ preserve $100\%$ of the data stream information ?

Depending what do you mean by preserve $100\%$ of data stream information. Usually $50\%$ of overlap is enough and you are not loosing too much information. Amount of data within your observation window is definitely sufficient for further analysis.

Although there are some windows, such as Hanning or modified Gaussian window which that are complementary. This allows you to perfectly reconstruct original signal from windowed segments.


What does the physical meaning of $\mathbf{U} \Sigma \ \mathbf{V}^{*} $ in the MIMO system? SVD is a method to cancel the interference in MIMO system,we often do the SVD of channel as below

\begin{align}
  \mathbf{H} &=
  \mathbf{U} \ \Sigma \ \mathbf{V}^{T} \\
%
\end{align}

So  $ \mathbf{U}$ means the beamforming direction of receiver? and $\mathbf{V}$ means the  beamforming direction of transmitter? and  the singular value in the diagonal of $\Sigma$ means the gain of channel ?Why can we use this to decrease or cancel the noise ?

Singular value decomposition (SVD) is a method of decomposing a matrix into two unitary matrices and a non-negative definite matrix of singular values.

In the context of MIMO communications, the $\mathbf{H}$ matrix defines the channel; it doesn't specify anything directly about the transmitter or the receiver.

This paper has the following diagram which shows one interpretation of the $\mathbf{U}$ and $\mathbf{V}$ matrices.  The SVD allows the channel to be decomposed into several parallel paths, with no crosstalk between them.


imperfect power control I want to simulate a direct sequence spread spectrum (DSSS) system and the receiver works with matched filter. I want the received power to be different for each user, in other words, SINR shouldn't be 0. Is there any idea how to do that in matlab?

I'm assuming this is the uplink channel since you want the received power to be different for each user.

A DSSS has a chip sequence that modulates the incoming signal. Basically the chips are just +1 and -1 and hence your TX signal power still is determined only by the incoming signal.

Generate random bits - A simple function randint.
Modulate them however you want - You can use qammod. - Check out the book "Fundamentals of Digital Communications" by Proakis - Chapter 3 for intricacies regarding energy Normalization for constellations.
Pass them through a Pulse Shaping filter of sufficient length. Check this link out for the same: http://www.mathworks.com/help/dsp/examples/pulse-shaping-filter-design.html
Multiply with the Chip Sequence. The way DSSS works is each user will have a Spreading code that is almost perfectly orthogonal to each other user. I think MATLAB can generate spreading sequences you might want to see: http://www.mathworks.com/help/comm/examples/spreading-sequences.html
Once you've generated the signals for all users, pass them each through a Flat fading Channel(you can consider different average received energy and generate the fade vector with variances appropriately - this will give you different SINRs) - You'd want your Chip period Tc to be very small compared to the DelaySpread of the channel and you can assume this which will hold the Flat Fading condition to be true. Multiply the fade vector for each user symbol by symbol with their corresponding signal points.
Receiver receives the Linearly combined signal and then you need to generate an equivalent length AWGN vector considering the average TX energy as the average TX energy you generated at the each user but the signals from each user suffered independent and different fades with different paths leading to different SINRs.
or AN ALTERNATIVE FOR (5,6): After spreading the sequences, alternatively you can multiply the whole signal for each user by a constant at the Transmitters before the fading channel leading to different energies(equivalent to different Path Losses for each signal) and linearly combine them all, Transmit and Receive on just one Fading Channel. Further you'd need to add AWGN considering the Average Signal energy of the signal you generated at the output of the Modulator(essentially the average energy of the Constellation).
The receiver now is like a Multi-User detector; you need to create copies of the received signal(as many number of users) and correlate with the Spreading sequences for each user.  You will essentially now remove the spreading code(Since synchronization is perfect) and you'll have the received signals for each user with different SINRs and so you can proceed to Matched Filtering using the same Filter you used for Pulse shaping at the Transmitter, Sampling, Demodulation using qamdemod and Detection using Maximal Likelihood(Minimum Distance).

You can refer to "Wireless Communications" by Andreas Goldsmith on further technical details in Chapter 3 for Fading and Chapter 13 and 15 for Multiuser Spread Spectrum Systems.

Note: Since you've titled this as Imperfect Power Control, I'd say the alternatives for 5,6 make more senes - the TXs are having bad knowledge about the channel and are not transmitting according to proper Powers.


1D filter for speed input and noisy position sensor I was tasked with designing a filter to smooth out a translating table. The system, we command a speed and a desired position, and measure it with two noisy position sensors. The commanded speed is also not always perfect, but it's usually within 10%. The slide can only move left and right.

My initial thought was to simply use a moving average (via circular buffer), and grow / shrink the array based on the speed. For high speeds, I would have a smaller array (10 or so measurements), and for lower speeds, I would have a large array (maybe 40 measurements). When stopped, I would have it grow to 100 or so. The reason for this is to decrease the phase offset at high speeds where resolution is not as important, and decrease the noise at lower speeds, when accuracy is more important. However, this doesn't really take into account the speed, so I feel like I'm throwing away useful information.

I'm figuring there has to be a better way. Any thoughts on how I would approach this problem more elegantly?

Just so we're on the same page, "high speed" is only about 3 inches per second, and I am sampling at 20Hz (but each time, I get two measurements). There is a ramp up/ ramp down portion of the translation such that whenever the table approaches the desired position, it slows to a crawl of 0.5 inches per second.

Thanks!

Sounds like an Extended Kalman Filter is the ideal candidate.

Bayes++ is a really simple drop-in library for implementing one, and the example here (https://github.com/Exadios/Bayes-/blob/master/PV/PV.cpp implements the question exactly (except it only takes in 1 observation)

Below is the code, in case it gets moved:

/*
 * Bayes++ the Bayesian Filtering Library
 * Copyright (c) 2002 Michael Stevens
 * See accompanying Bayes++.htm for terms and conditions of use.
 *
 * $Header: /cvsroot/bayesclasses/Bayes++/PV/PV.cpp,v 1.17.2.3 2005/04/07 06:39:22 mistevens Exp $
 */

/*
 * Example of using Bayesian Filter Class to solve a simple problem.
 *  The example implements a Position and Velocity Filter with a Position observation.
 *  The motion model is the so called IOU Integrated Ornstein-Uhlenbeck Process Ref[1]
 *    Velocity is Brownian with a trend towards zero proportional to the velocity
 *    Position is just Velocity integrated.
 *  This model has a well defined velocity and the mean squared speed is parameterised. Also
 *  the velocity correlation is parameterised.
 *
 * Two implementations are demonstrated
 *  1) A direct filter
 *  2) An indirect filter where the filter is preformed on error and state is estimated indirectly
 * Reference
 * [1] "Bayesian Multiple Target Tracking" Lawrence D Stone, Carl A Barlow, Thomas L Corwin
 */

#include "BayesFilter/UDFlt.hpp"
#include "BayesFilter/filters/indirect.hpp"
#include "Test/random.hpp"
#include <cmath>
#include <iostream>
#include <boost/numeric/ublas/io.hpp>

namespace
{
    using namespace Bayesian_filter;
    using namespace Bayesian_filter_matrix;

    // Choose Filtering Scheme to use
    typedef UD_scheme FilterScheme;

    // Square
    template <class scalar>
    inline scalar sqr(scalar x)
    {
        return x*x;
    }

    // Random numbers from Boost
    Bayesian_filter_test::Boost_random localRng;

    // Constant Dimensions
    const unsigned NX = 2;          // Filter State dimension   (Position, Velocity)

    // Filter Parameters
    // Prediction parameters for Integrated Ornstein-Uhlembeck Process
    const Float dt = 0.01;
    const Float V_NOISE = 0.1;  // Velocity noise, giving mean squared error bound
    const Float V_GAMMA = 1.;   // Velocity correlation, giving velocity change time constant
    // Filter's Initial state uncertainty: System state is unknown
    const Float i_P_NOISE = 1000.;
    const Float i_V_NOISE = 10.;
    // Noise on observing system state
    const Float OBS_INTERVAL = 0.10;
    const Float OBS_NOISE = 0.001;

}//namespace

/*
 * Prediction model
 * Linear state predict model
 */
class PVpredict : public Linear_predict_model
{
public:
    PVpredict();
};

PVpredict::PVpredict() : Linear_predict_model(NX, 1)
{
    // Position Velocity dependance
    const Float Fvv = exp(-dt*V_GAMMA);
    Fx(0,0) = 1.;
    Fx(0,1) = dt;
    Fx(1,0) = 0.;
    Fx(1,1) = Fvv;
    // Setup constant noise model: G is identity
    q[0] = dt*sqr((1-Fvv)*V_NOISE);
    G(0,0) = 0.;
    G(1,0) = 1.;
}

/*
 * Position Observation model
 * Linear observation is addative uncorrelated model
 */
class PVobserve : public Linrz_uncorrelated_observe_model
{
    mutable Vec z_pred;
public:
    PVobserve ();
    const Vec& h(const Vec& x) const
    {
        z_pred[0] = x[0];
        return z_pred;
    };
};

PVobserve::PVobserve () :
    Linrz_uncorrelated_observe_model(NX,1), z_pred(1)
{
    // Linear model
    Hx(0,0) = 1;
    Hx(0,1) = 0.;
    // Observation Noise variance
    Zv[0] = sqr(OBS_NOISE);
}

void initialise (Kalman_state_filter& kf, const Vec& initState)
/*
 * Initialise Kalman filter with an initial guess for the system state and fixed covariance
 */
{
    // Initialise state guess and covarince
    kf.X.clear();
    kf.X(0,0) = sqr(i_P_NOISE);
    kf.X(1,1) = sqr(i_V_NOISE);

    kf.init_kalman (initState, kf.X);
}

int main()
{
    // global setup
    std::cout.flags(std::ios::scientific); std::cout.precision(6);

    // Setup the test filters
    Vec x_true (NX);

    // True State to be observed
    x_true[0] = 1000.;  // Position
    x_true[1] = 1.0;    // Velocity

    std::cout << "Position Velocity" << std::endl;
    std::cout << "True Initial  " << x_true << std::endl;

    // Construct Prediction and Observation model and filter
    // Give the filter an initial guess of the system state
    PVpredict linearPredict;
    PVobserve linearObserve;
    Vec x_guess(NX);
    x_guess[0] = 900.;
    x_guess[1] = 1.5;
    std::cout << "Guess Initial " << x_guess << std::endl;

    // f1 Direct filter construct and initialize with initial state guess
    FilterScheme f1(NX,NX);
    initialise (f1, x_guess);

    // f2 Indirect filter construct and Initialize with initial state guess
    FilterScheme error_filter(NX,NX);
    Indirect_kalman_filter<FilterScheme> f2(error_filter);
    initialise (f2, x_guess);

    // Iterate the filter with test observations
    Vec u(1), z_true(1), z(1);
    Float time = 0.; Float obs_time = 0.;
    for (unsigned i = 0; i < 100; ++i)
    {
        // Predict true state using Normally distributed acceleration
        // This is a Guassian
        x_true = linearPredict.f(x_true);
        localRng.normal (u);        // normally distributed mean 0., stdDev for stationary IOU
        x_true[1] += u[0]* sqr(V_NOISE) / (2*V_GAMMA);

        // Predict filter with known pertubation
        f1.predict (linearPredict);
        f2.predict (linearPredict);
        time += dt;

        // Observation time
        if (obs_time <= time)
        {
            // True Observation
            z_true[0] = x_true[0];

            // Observation with addative noise
            localRng.normal (z, z_true[0], OBS_NOISE);  // normally distributed mean z_true[0], stdDev OBS_NOISE.

            // Filter observation
            f1.observe (linearObserve, z);
            f2.observe (linearObserve, z);

            obs_time += OBS_INTERVAL;
        }
    }

    // Update the filter to state and covariance are available
    f1.update ();
    f2.update ();

    // Print everything: filter state and covariance
    std::cout <<"True     " << x_true << std::endl;
    std::cout <<"Direct   " << f1.x << ',' << f1.X <<std::endl;
    std::cout <<"Indirect " << f2.x << ',' << f2.X << std::endl;;
    return 0;
}


Combined impulse responses I'm looking for a simple way to show that you can form a single impulse response that is the equivalent $M$ other impulse responses, ie:

$$
h\left(t\right)=h_0\left(t\right)\star h_1\left(t\right)...\star h_{M-1}\left(t\right)
$$

I know convolution is multiplication in the frequency domain, so I cooked up this "equation":

$$
h\left(t\right)=\mathcal F^{-1}\left\{\prod_{i=0}^{M-1} \mathcal F\left\{h_{i}\right\}\right\}(t)
$$

With $\mathcal F^{-1}$ being the inverse Fourier transform operator, and $\mathcal F$ being the Fourier transform operator. Expressing myself in math notation isn't really my strong side, what I intended to show was that the combined impulse response $h(t)$ equals the sequences multiplied together in the frequency domain.

Does it make any sense?
Or is there perhaps a simple way to show this?
What about circular convolution in the frequency domain, and infinite impulse responses in continuous time?

You can readily use the associativity of the convolution to explain this fact. Assume you send a signal $x(t)$ through a chain of convolutions $h_i(t)$:

$$ y(t) = (((x(t)*h_1(t))*h_2(t))*h_3(t)) $$

By associativity, this becomes

$$ y(t) = x(t) * ((h_1(t)*h_2(t))*h_3(t)) $$

I.e. you can see this operation as if sending your signal $x(t)$ through a single filter with impulse response $h(t)=h_1(t)*h_2(t)*h_3(t)$.


Randomized prefix code An input symbol is randomly chosen from a set of symbols with known associated probabilities and encoded using a prefix code that minimizes the expected codeword length. Is it possible to further reduce the total expected codeword length by randomly choosing the prefix code from a set of prefix codes with associated probabilities?

No, it is not possible. A proof follows.

Let $\mathrm{E}(L \mid 0)$ be the expected codeword length for the prefix code $0$ that minimizes it. The answer to the question is "no" if it is always true that:

$$\mathrm{E}(L \mid 0) \le \mathrm{E}(L),\tag{1}$$

where $\mathrm{E}(L)$ is the total expected codeword length using a randomly chosen prefix code $i$ from a set of $N$ prefix codes denoted by $i\in\{1, 2, 3, \ldots, N\}$. We can divide both sides of Eq. 1 by $\mathrm{E}(L \mid 0)$, which is positive, so we could equivalently show that:

$$\Leftrightarrow 1 \le \frac{\mathrm{E}(L)}{\mathrm{E}(L \mid 0)}.\tag{2}$$

Using the law of total expectation:

$$\frac{\mathrm{E}(L)}{\mathrm{E}(L \mid 0)} = \frac{\sum_{i=1}^{N}\mathrm{E}(L \mid i)P(i)}{\mathrm{E}(L \mid 0)} = \sum_{i=1}^{N}\frac{\mathrm{E}(L \mid i)}{\mathrm{E}(L \mid 0)}P(i),\tag{3}$$

where $\mathrm{E}(L \mid i)$ is the expected codeword length for a prefix code $i$ with an associated non-zero probability $P(i)$. Because the prefix code $0$ minimizes the expected codeword length, we have $\mathrm{E}(L \mid 0) \le \mathrm{E}(L \mid i)$ $\Rightarrow$ $\mathrm{E}(L \mid i)/\mathrm{E}(L \mid 0) \ge 1$ $\Rightarrow$ $\mathrm{E}(L \mid i)/\mathrm{E}(L \mid 0) - 1 \ge 0$ for all $i$. It follows that:

$$\frac{\mathrm{E}(L)}{\mathrm{E}(L \mid 0)} = \sum_{i=1}^N P(i) + \sum_{i=1}^{N}\left(\frac{\mathrm{E}(L \mid i)}{\mathrm{E}(L \mid 0)}-1\right)P(i) \ge \sum_{i=1}^N P(i)= 1,\tag{4}$$

where the second sum is never negative. This proves Eq. 2. Therefore, the answer to the question is "no".


Filtering a digital signal online in real-time using python I am currently trying to apply a bandpass filter to a signal in real-time. There are samples coming in with a constant sampling rate and i would like to calculate the corresponding bandpass filtered signal.

What would be the best way to do this? Do I have to filter the whole (or at least a huge bit) of the signal every time a few new samples came in or is there a way (like the sliding DFT) where it is possible to efficiently determine the new part of the filtered signal?

I would like to use a butterworth filter (for offline analysis I am currently using scipy's butter and lfilter). I know that this function can return a filter delay, but I don't know how to use it to get a constant signal.

Do I have to filter the whole (or at least a huge bit) of the signal every time a few new samples came in or is there a way (like the sliding DFT) where it is possible to efficiently determine the new part of the filtered signal?

Digital filters don't work like that -- basically, classical FIR or IIR can work on every single new sample. You should really read up on what these filters are, and how people model them.

  I would like to use a butterworth filter

Well, there's plenty implementations of that out there,

  I am currently using scipy's butter and lfilter

of which you already know one!

Now, a butterworth filter is a recursive thing, so to calculate the next part of your sampled signal, you'll need the last state. That's exactly the "filter delay state zi" that lfilter returns and can take in on the next call as zi parameter.

  but I don't know how to use it to get a constant signal.

I think you mean "to achieve continuous filtering".

Now, that being said, the point is that you're setting yourself up for writing your own streaming architecture. I wouldn't do that. Use an existing framework. For example, there's GNU Radio, which lets you define signal processing flow graphs in Python, and also is inherently multithreaded, uses highly optimized algorithm implementations, has a lot of in- and output facilities, and comes with a huge library of signal processing blocks, which can be written in Python or C++, if you happen to need to do that.

For example, a flow graph that takes in samples from a sound card, butterworth-filters them and writes them to a file is:

#!/usr/bin/env python2
# -*- coding: utf-8 -*-
##################################################
# GNU Radio Python Flow Graph
# Title: Butterworth Test
# Generated: Mon Feb  8 16:17:18 2016
##################################################

from gnuradio import audio
from gnuradio import blocks
from gnuradio import eng_notation
from gnuradio import filter
from gnuradio import gr
from gnuradio.eng_option import eng_option
from gnuradio.filter import firdes
from optparse import OptionParser

class butterworth_test(gr.top_block):

    def __init__(self):
        gr.top_block.__init__(self, "Butterworth Test")

        ##################################################
        # Variables
        ##################################################
        self.samp_rate = samp_rate = 48000

        ##################################################
        # Blocks
        ##################################################
        # taps from scipy.butter!
        self.iir_filter_xxx_0 = filter.iir_filter_ffd(([1.0952627450621233e-05, 0.00013143152940745496, 0.0007228734117410033, 0.0024095780391366808, 0.005421550588057537, 0.008674480940892064, 0.010120227764374086, 0.008674480940892081, 0.005421550588057554, 0.0024095780391366955, 0.0007228734117410089, 0.00013143152940745594, 1.0952627450621367e-05]), ([1.0, -4.4363862740719835, 10.215121830052535, -15.374408118154847, 16.57333784740102, -13.325056987818655, 8.133543488903097, -3.77641064765334, 1.3181452681671835, -0.3361758629961047, 0.05930166356243964, -0.0064815521348275, 0.00033130678123743994]), False)
        self.blocks_file_sink_0 = blocks.file_sink(gr.sizeof_float*1, "", False)
        self.blocks_file_sink_0.set_unbuffered(False)
        self.audio_source_0 = audio.source(samp_rate, "", True)

        ##################################################
        # Connections
        ##################################################
        self.connect((self.audio_source_0, 0), (self.iir_filter_xxx_0, 0))
        self.connect((self.iir_filter_xxx_0, 0), (self.blocks_file_sink_0, 0))

def main(top_block_cls=butterworth_test, options=None):

    tb = top_block_cls()
    tb.start()
    try:
        raw_input('Press Enter to quit: ')
    except EOFError:
        pass
    tb.stop()
    tb.wait()

if __name__ == '__main__':
    main()

Note that this code was auto-generated from a graphical flow graph that I just clicked together using the gnuradio-companion program:

If you want to know more about how to implement signal processing flow graphs in Python, go for  the GNU Radio Guided Tutorials.

EDIT: I liked @Fat32's answer pretty much! What he describes as a double buffering architecture is pretty close to what GNU Radio does:

An upstream block produces samples in sample chunks of arbitrary sizes, writes them to the output ring buffer (which is represented as an arrow in above picture), and notifies its downstream blocks that there's new data.

The downstream block gets notified, checks whether there's enough space in its output buffer to process the samples that are in his input ring buffer (which is the same as the output buffer of the upstream block), processes these. When finished, it informs the upstream block(s) that it has used up input ring buffer (which can then be reused by the upstream blocks as output), and the downstream blocks about new samples being available.

Now, GNU Radio being multi-threaded, the upstream block might already be producing samples again; in a normal GNU Radio application, nearly all blocks are "active" simultaneously and things scale pretty well on multi-CPU machines.

So, GNU Radio's main job is to give you this buffer infrastructure, the notifying and thread housekeeping, the clear signal processing block API and something to define how everything is connected, so you don't have to write what Fat32 describes in her/his post yourself! Notice that doing sample stream marshalling isn't that easy to do properly, and GNU Radio takes the hardness out of it and lets your concentrate on what you want to do: DSP.


Which of these features may not be needed or redundant? I am trying to extract some features of the accelerometer sensor and during the feature extraction each 5 second segment in the array is windowed and features are generated from these windows. Features extracted include:
•   Mean of the magnitude data
•   Maximum value.
•   Minimum value.
•   Standard Deviation value for each axis (x, y, and z).
•   Average Standard Deviation over 3 axes.
•   Root mean squared acceleration (RMS).
•   The sum of height of frequency component below 5 Hz.
•   Number of peaks in spectrum.

I am not sure if mixing all these is adding redundant values and affect the final results.

Thanks

Whether or not these features are redundant can only be concluded regarding your target variable that you want to deduce. Mathematically speaking, these features are not redundant, since you can not directly compute one from the other without additional information.

Thus each feature adds additional information.

For redundancy check in general, you can compute the mutual information or - often used - check for correlation.


Best Metric to Compare Sparsity of Vectors I solved the Basis Pursuit Denoising Problem looking for a sparse solution (I am in compressive sensing):
$$ {x}^{\ast} = \arg \min_{x} \left\{ \frac{1}{2} {\left\| A x - y \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} \right\} $$

for $ 100 $ different $ \lambda $ values with ADMM.

Now I have $100$ different $x^*$ and want to find "the best one" out of them. Moreover, I solved several optimization problems for CS under noise and noisefree cases (e.g. also Lasso, etc.)

What is the best measure to determine sparsity? Found the paper Comparing Measures of Sparsity by Niall Hurley and Scott Rickard, saying it's the Gini index. Do you agree? In my opinion sparsity in this case is not all, isn't it? It's also about solving the equation as good as possible, e.g. having a small norm on:

$$ {\left\| A {x}^{\ast} - y \right\|}_{1} $$

What do you recommend?

I am sorry I cannot comment your answer due to my low reputation. Gini and your suggested sparsity ratio ($l_1(x)/l_2(x)$) both give me the same value for $\lambda$. But

The problem I still see is that I cannot take into account how well the vector is solving the equation $Ax-y$. I would like to combine the residuum $l_1(A\hat{x}-y)$ and the sparsity $l_1(\hat{x})/l_2(\hat{x})$ into a new metric value. Do you know something from the literature?
Gini only works with positive values, but my measurement vector has lots of negative values. In my eyes, this is a heavy disadvantage because many sparse vectors can have negative values.


I want to get tilt sensor value for a given time t_0, but the tilt sensor data points are discrete. How do I interpolate between the discrete points? I have 3-dimensional readings from a tilt sensor (specifically these are rotational angles about X, Y, and Z axes) over time. Let's call these angles S. I want to infer S at a specific time t_0, but each reading is discrete and are at different time points across axes. Here is graphic1 to illustrate.

graphic1

I'm thinking the solution will be something like:

Construct a window in time centred on t_0, then smooth/interpolate between the discrete points within that time window to obtain the 3-dimensional value at t_0, like in graphic2.

graphic2

However, I'm not sure on the following:

How wide should the window be such that no information is lost? (maybe something about Nyquist theorem?)
What smoothing algorithm to use?
The readings themselves are noisy.

Am I thinking about this problem the right way? I'm open to other alternatives.

Also worth nothing that I have to write an implementation in TypeScript.

Assuming that the signal is prefiltered and uniformly sampled according to Shannon Nyquist the textbook solution is to upsample it  using some approximation of the ideal (frequency domain) brickwall lowpass filter, equivalent with a (windowed) sin(x)/x function in the time domain.

The problem is similar to showing a smooth (zoomed in) waveform of discretely sampled audio in an audio editor. Perhaps you can find the code you need in eg Audacity?

If the signal contains noise, you might be able to trade noise suppression for signal distortion by knowing the statistics of both and using a Wiener filter.

As this seems to be a control systems project, perhaps you can use a Kalman filter.


the best pre emphasis for MMSE equalizer The system consists of multiple low-pass filter and noise, $h_1$, $h_2$, $h_3$, ..., $n_1$, $n_2$, $n_3$, ...
At the transmitter, the power of signal is fixed; and at the receiver, MMSE equalizer is used to minimize the MSE.
My question is that if we can pre-emphasis the signal in the transmitter to improve the overall performance. If the distribution of signal and noise is known, what is the optimal pre-emphasis system?

welcome.
As MMSE estimators are unbiased, try reducing Cramér–Rao bound and maybe your MMSE equalizer can achieve a better performance. Note that in general the MMSE estimator is known to be asymptotically efficient.
Note also that you must specify criteria before say something is "optimal". Even so, it is usually not easy to prove the optimality.


A feature extraction method: Principal Component Analysis (PCA) I read some articles about PCA and I think the nice way to summarize it is:

  It is a way of identifying patterns in data, and expressing the data in such a way as to highlight their similarities and differences.

So in my code, we subtract the initial data from the mean and find principal components, then  we multiply the data by a variable called e as for our dimension reduction purposes.

What I don't understand here is that When we want to highlight the similarities and differences, it happens both about noise data, the data which is too big or too small, will be soften. What about the edges in our image? Will the same thing happen to those data? I mean, when we have a change of color from white to black and we don't want to lose such information, will PCA soften the data as well as noise? If yes, this is not what we want, is it? Is there a solution to solve this problem ?

Ok, it sounds like you are trying to do eigenfaces, right?

In that case, you have to think of your face images as points in a very high-dimensional space. For example, if your images are 32x32, then the space has 32 * 32 = 1024 dimensions.  Operating in so many dimensions is very difficult, because distances between points become almost meaningless.  With so many dimensions any point is very far away from any other point, even if the two points represent the same face.

That's why you want to reduce the dimensionality, and PCA is one way to do it. If you go from 1024 dimensions to, say, 10, then the distances will become more meaningful, which fits in with your idea of "highlight their similarities and differences".  However, PCA by itself will not tell you of any patterns in your data, beyond the directions of greatest variance.  PCA will project your data into a lower-dimensional space with decorrelated axes, but you have to do more work to actually discover the patterns in the data, such as train a classifier, or run a clustering algorithm.

Now, to understand what happens to the edges and transitions from light to dark, you have to display your eigenvectors as images and look at them. You will see that you don't really preserve the precise edge information, but you do keep course-level information about light and dark areas.


Difference between the two forms of input-output relationships of an LTI system $$y(t)=f(t)*h(t)\tag{1}$$

$$y(t)=H(s)e^{st}\tag{2}$$

$$H(s)=\int_{-\infty}^{\infty}h(t)e^{-st}dt\tag{3}$$

Let $f(t)$ in Eqn $(1)$ be $e^{st}$. In many worked out examples, I have found that the two equations/formulae give different responses. To be specific, the actual convolution (Eqn $(1)$) is giving the zero-state response, while the Eqn $(2)$ involving $H(s)$ is giving the forced response. Recognising that Forced response and Zero-State response are not equal (Z-S responses may involve characteristic modes of the system, while the Forced Response is the non-characteristic mode part of the Total response), I am confused what is the intuitive difference between the two equations? Which equation should be primarily used for finding the response to an input to a system? Note that in the examples, the inputs $e^{st}$ were multiplied by $u(t)$.

For an LTI system with impulse response $h(t)$ and transfer function $H(s)=\mathcal{L}\{h(t)\}$, an input signal $x(t)=e^{st}$ (note: no multiplication with the unit step function $u(t)$), results in an output

$$\begin{align}y(t)&=(x\star h)(t)\\&=\int_{-\infty}^{\infty}h(\tau)x(t-\tau)d\tau\\&=\int_{-\infty}^{\infty}h(\tau)e^{s(t-\tau)}d\tau\\&=e^{st}\int_{-\infty}^{\infty}h(\tau)e^{-s\tau}d\tau\\&=e^{st}H(s)\tag{1}\end{align}$$

This is a very basic result concerning eigenfunctions of LTI systems.

If the input signal is given by $x(t)=e^{st}u(t)$ then the output is not given by $(1)$, but it needs to be determined from the convolution integral:

$$\begin{align}y(t)&=\int_{-\infty}^{\infty}h(\tau)e^{s(t-\tau)}u(t-\tau)d\tau\\&=\int_{-\infty}^th(\tau)e^{s(t-\tau)}d\tau\\&=e^{st}\int_{-\infty}^th(\tau)e^{-s\tau}d\tau\end{align}$$

Note that in both cases you get the zero-state response because convolution with an impulse response cannot give you anything else.


Data transmission with SNR<0dB? I encounter a situation that it is very beneficial to limit the signal strength below noise floor.

I know that the Shannon–Hartley theorem suggest that it is possible:

I searched over the internet but cannot find any direct pointer.

What keywords should I look for? Are there any compiled resources on the topic?

Thanks

There may be some confusion about the definition of SNR.

In some applications you may be interested in reducing the ratio of the signal power spectral density to noise power spectral density.  "SNR" in this case refers to the ratio of the signal to noise power when you measure both over the same bandwidth.  In these cases, as Juancho pointed out, direct-sequence spread spectrum is a classical technique.  It deliberately spreads the signal power over a wider bandwidth.

The SNR used in the Shannon-Hartley Theorem refers to the SNR per channel usage.  To oversimplify the matter, it refers to the SNR per symbol after you demodulate the signal.  Spread spectrum does not decrease the SNR needed to achieve a given data rate; it just spreads the signal power in the SNR calculation over a wider bandwidth.  It is true that you can communicate with SNR lower than 0 dB: one typically does so in a practical environment using simple modulation (e.g., BPSK, QPSK) and a low code rate with a strong-performing code.  Whether one spreads the signal or not is a separate decision.


How do I start studying Gabor Filters ? I know very basic image processing. Please suggest me how do I and where do I start learning Gabor Filters ?

You know what I actually found really useful, Gabor's original paper, "Theory of Communication":

http://redwood.berkeley.edu/w/images/b/b6/Gabor.pdf

Of particular interest is how Gabor used banks of reeds, tuned to different frequencies, to visualise his filters. For me this led to a good conceptual understanding, especially of the uncertainty principle, before diving into the mathematics.


Are pulse shaping and windowing functions the same in concept? I read a digital communication problem and they interchange terms as windowing function and pulse shaping and refer to it as same entity. For example the following windowing function with Transition time $T_R$ used to smooth transitions between symbols is also refered to as pulse shaping

$w_T(t)= \left\{ \begin{array}{ll}
         \sin^2(\frac{\pi}{2}(0.5+\frac{t}{T_{TR}})) & \mbox{if $-T_{TR}/2 < t<T_{RT}/2$};\\
        1 & \mbox{if $T_{TR}/2 < t<T-T_{TR}/2$}\\
 \sin^2(\frac{\pi}{2}(0.5-\frac{t-T}{T_{TR}})) & \mbox{if $T-T_{TR}/2 < t<T+T_{TR}/2$};\\
\end{array} \right. $

My understanding is that a pulse shaping function limits signal bandwidth to transmission bandwidth.

How come the two are refered to as same?

Thanks

Windowing and pulse shaping are very similar, and are implemented in similar ways, but their purposes are different.

Windowing means multiplying a discrete signal by a pulse with certain spectral properties before calculating its Fourier transform. The purpose is reducing the effect of discontinuities at both ends of the discrete signal.
In QAM and related digital signals, the information is transmitted in the amplitudes of a sequence of pulses. Pulse shaping means choosing a pulse with the desired spectral and orthogonality properties.

Note that some digital communications textbooks assume (many times implicitly) a rectangular pulse shape, which is then converted to another pulse shape (using a filter, as in the example in your question). In my opinion, this formulation is unncecessarily complex -- it is better to explicity select  a pulse shape from the start.

Note also that pulse shaping may be understood as the convolution of a stream of deltas with a pulse shape; in that sense, pulse shaping is a sort of "dual" of windowing.


Checking the "cleanliness" of an image for Optical Character Recognition Is there some mechanism to say whether the image is "clean" for the purpose of image processing?

For example, in the following image, the topmost row is the cleanest and the bottom one is the least clean. As a result, the top-image gives accurate output; the bottom-one gives less accurate.

Examples of noisy and noise-less images

The purpose is to use the image for Optical Character Recognition.
And, I wanted to check beforehand, if the image is sufficiently clean or not.

Your images are gray-scaled. In order to perform the OCR you might want to first threshold it to apply OCR on the binary image.

The effectiveness of the thresholding will definitely affect the performance of the algorithm. A simple way is say Otsu's thresholding. In your case, you can think of this as character vs. non-character pixels. And the measure of the cleanliness can be treated as the amount of energy (such as class variance) of the non-character pixels.

Unfortunately the measure I described is post applying your (first step) of OCR. Hence, I am not sure if this helps you.


How to simplify an exponential functions in MATLAB? I have many defined functions as below and trying to simplify it with MATLAB simplify() function
$g(t) = 14.9je^{j31°} e^{j2t} - 14.9je^{-j31°} e^{-j2t}$
How can we simplify these kind of functions using MATLAB?
Thank you.

Symbolic math is included in MATLAB, but rarely used in signal processing. There are cases, nevertheless, when we wish to use symbolic math packages to simplify some expressions. One such alternative is also Wolfram-Alpha web application.
Anyway, apart from a misunderstanding, you indeed have a (simple) mathematical problem:
Given a complex number $z$ , you have
$$ z + z^* = 2 ~\mathcal{Re}\{ z \} $$
or a complex-valued function $f(t)$
$$ f(t) + f^*(t) = 2 ~\mathcal{Re}\{ f(t) \} $$
If you look at your function $g(t)$, you will see that it's a sum of a complex-valued function $f(t)$ and its conjugate $f^*(t)$, where
$$f(t) = 14.9 ~j~e^{j 31^o} ~ e^{j 2t} $$
and its conjugate
$$f^*(t) = - 14.9 ~j~e^{-j 31^o} ~ e^{-j 2t} $$
hence
$$g(t) = f(t)+f^*(t) = 2 \mathcal{Re}\{ f(t) \} $$
Then you can see that
$$\mathcal{Re}\{ f(t) \} = \mathcal{Re}\{ 14.9 ~j~e^{j 31^o} ~ e^{j 2t} \} =\mathcal{Re}\{ 14.9 ~e^{j (121^o +2t)} \} = 14.9 \cos(121^o +2t)$$
Then your $g(t)$ is:
$$ g(t) = 2 \times \mathcal{Re}\{f(t)\} = 2 \times 14.9 \cos(121^o +2t) $$


Sound reflections plot I have sound files that contain recordings of an impulse signal in a room.

I'd like to make such a plot from it. I basically want to compare early reflections.

How do I go about it? Is just taking the waveform, pulling the absolute values and plotting it in bar/impulse form the right way? Or is this just an idealized form? If this is all wrong, could anyone point me towards the correct way of comparing early reflections in a plot form?

This is an idealized case of the echogram that simply depicts the times of arrival and energy of reflections. You can see this type of plot in literature and as an output in modelling software such as CATT-A and ODEON.

If you have a real-life impulse response then you can obtain information about reflections by representing the IR in logarithmic scale, i.e.:


Linear systems: Square root of input product Hi guys i'm studying signals and systems, and my professor ask us if this signal is linear or not
$$y(t) = \big[x(t − 1)x(t + 1)\big]^{\frac 12}$$
the fact that is in the form of $x\cdot x$ told me that that's not linear but the square root confuses me can you explain to me if this is linear or not?


Derivative calculation I have a series of data (single array). If I take this data and plot it, I can see that there are multiple peaks. However If I zoom in to a section of data, I see that there is substantial noise.

I like to detect how many peaks in the data using little CPU time and energy. I thought of detecting the slopes of these peaks (could be sharp or smooth rising or falling) and look at the number of slopes to determine the number of peaks. (2 slopes positive and negative for every peak)

Any pointers how I can calculate the slopes in the presence of noise?

The code will go into an embedded system, memory is limited so preferably I like to implement something that doesn't require any significant data copy.

It depends on the type of noise and type of signal.  Show an example if you want a good answer.  But, that said, in general you probably want to low-pass filter the signal.  If I were you, I'd take a Fourier power spectrum to see if most of the noise is high frequency, and the signal I care about mostly in a lower range.  If they overlap, oh well that's life. I would have to think more about things.

One low-pass filter that's good for noisy signal in many cases is the Savitzky-Golay filter.  It is described in Numerical Recipes, and for Python there's a function in the Python Numpy Cookbook.  It is merely a convolution with a small kernel.   You pick the window size based on the width of the peaks or other features, wide enough to mush out the noise, but not wider than the features.  It can be small, say five points, or bigger like dozens, a hundred maybe.

You also pick a polynomial order - usually I use 2 or 4.  Order 2 is fine for when the window is small, < 10 points or it spans less than half a cycle or so (if your signal resembles a sine) while order 4 is better at matching distorted peak shapes, but likes to have around 9 or more points.  But a lot depends on the shape and frequency of the noise.

As other say in the comments, finding derivatives probably isn't the best strategy, but if you want find derivatives anyway, the Savitzky-Golay filter can do that - simultaneously smoothing and reporting the derivative instead of the signal.


Development of Hilbert transform relationship In the development of Hilbert transform relationships, Prof. Oppenheim has chosen
\begin{equation}
\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)\sum_{k=-\infty}^{\infty}\delta(\omega-\theta-2\pi{}k)d\theta
=X_R\left(e^{j\omega}\right)
\end{equation}
in his book "Discrete-Time Signal Processing, 2e" Chapter 11 pp807, Eq. 11-24. I failed to understand how he arrives at it.

Does anybody any idea how has derived, please let me know?

I tried some BUT IT DID NOT GIVE the above relationship:

1:

\begin{align*}
\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)
\sum_{k=-\infty}^{\infty}\delta(\omega-\theta-2\pi{}k)d\theta&=\frac{1}{2\pi}\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)2\pi\sum_{k=-\infty}^{\infty}
\delta(\omega-\theta-2\pi{}k)d\theta\\&
=\frac{1}{2\pi}\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)e^{j\theta{}n}d\theta\\&
=x_e(n)
\end{align*}

OR

2:

\begin{align*}
\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)
\sum_{k=-\infty}^{\infty}\delta(\omega-\theta-2\pi{}k)d\theta&=\sum_{k=-\infty}^{\infty}\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)
\delta(\omega-\theta-2\pi{}k)d\theta\\&
=\sum_{k=-\infty}^{\infty}X_R\left(e^{j(\omega-2\pi{}k)}\right)
\end{align*}

Your second try was already close. However, you treated the integral as if its limits were $-\infty$ and $\infty$. With the given integration limits there is only one delta impulse in the integration interval, i.e. only one index of the sum is relevant, all others are outside the integration interval. So if for some given value $\omega$ this index is $k=l$, then the result is

$$\int_{-\pi}^{\pi}X_R\left(e^{j\theta}\right)
\sum_{k=-\infty}^{\infty}\delta(\omega-\theta-2\pi{}k)d\theta=X_R(e^{j(\omega-2\pi l)})=X_R(e^{j\omega}\cdot e^{-j2\pi l})=X_R(e^{j\omega})$$


How to remove/reduce noise by using the Morlet wavelet? I'm trying to use Morlet wavelet to reduce noise in my signal.
I found formula for Morlet wavelet

$$\Psi(x) = \frac{1}{\sqrt{\pi\cdot \textrm{bandwidth}}} \cdot \exp\left(i \cdot 2\pi \cdot \textrm{centerFreq} \cdot x\right)\cdot\exp\left(\frac{x^2}{\textrm{bandwidth}}\right)$$

(sb told me I shoud use formula for complex Morlet wavelet)

I also got CWT formula

$$C = \int f(t) \cdot \frac{1}{\sqrt{\textrm{scale}}} \cdot \Psi\left(\frac{t-\textrm{shift}}{\textrm{scale}}\right)dt$$

I convert both formulas to java code

public double morletRe(double x)
{
return (1/Math.sqrt(Math.PI*fb))*Math.cos(2*Math.PI*fc*x)* Math.exp(x*x/fb);
}

public double[] cwt(double[] data,double scale, double position)
{
double[] newData = new double[data.length];
    for (int t = 0; t<data.length; t++ )
    {
    newData[t] = data[t]*(morletRe(t-position)/scale)[0])/Math.sqrt(scale);
    }
return newData;
}

But in the results I get a strange chart which is completely not related to my data,
is anyone can point me what I'm doing wrong?
regards.

I'm not sure of the reliability of your source... I worked with wavelets on my master's thesis and the basic wavelet transform is a convolution of the wavelet an your signal.
In most practical cases, it's best to go through Fast Fourier Transform (FFT), as it's mathematically equivalent.

With WT being the wavelet transform, and FFT() being the Fast Fourier Transform of what's inside the brackets, iFFT() is the inverse-FFT. You would compute the WT as such:

FFT(WT) = FFT(signal) * FFT(wavelet)
WT = iFFT(WT)

Now, for the wavelet formula itself, I believe there has to be a specific scaling factor to get a reliable result. Here's what we used:

$$\Psi = exp\Big(-\frac{ik_0x}{\lambda} - \frac{x^2}{2\lambda^2}\Big)$$

Where $\lambda$ is the wavelength (or scale) and $k_0 = 5.336$ (source: J.F. Kirby, Computers & Geosciences 31 (2005) 846–864)
I did not include the normalization factor because you just have to divide your wavelet by it's highest value, and subtract its average, and you get the same result with less headhache.
Off course, you can interchange $1/\lambda$ by $2\pi*frequency$ in the formula.

In order to qualify for the "Morlet" name, I think you have to use the same wavelength in both places in the formula. Else, you over-sample or under-sample the frequency part of your signal. $k_0$ is the best ratio between the two, as explained by Kirby in his article.


Is a wavelet-based correlation measure worth any additional computational overhead? I have used both correlation and coherence as measures of correlation between signals.  I was thinking that a time-frequency approach would give me the best of these worlds.

My question is whether this extra data adds enough to the overall picture of the signal to justify the increased computational cost associated with doing the wavelet transforms as part of the calculation?

Reference: an ArXiv paper: "A cross-correlation technique in wavelet domain for detection of stochastic gravitational waves" by S.Klimenko, G.Mitselmakher, A.Sazonov

First off, you should use whichever tool is appropriate for the job.  Correlation vs coherence vs wavelet-based correlation are all different things, so this question is kind of like asking "Which is better? Screwdrivers or hammers?"  It depends on what you're trying to do, and whether you care about similarity in time, frequency spectra, or both.

Second, I have only a minimal understanding of wavelets, but your assumption that wavelets require more computation might be wrong.  The Fast Fourier Transform takes $O(n\log n)$ operations, while the Fast Wavelet Transform takes $O(n)$.  So the wavelet method might actually require less computation, depending on whether you can use the information you get out of it.

Empirically, producing n outputs from n real inputs, the multi-level wavelet transform in PyWavelets becomes faster than NumPy's FFT when n is greater than about 4096.

However

It's Python, and the two implementations could be very differently efficient.  I don't even know if wavedec() would be considered a Fast Wavelet Transform.  They use the abbreviation DWT in their documentation.  Are Haar DWT and FWT the same thing?
The time varies depending on the wavelet used.  Meyer wavelet takes 6 times as long as Daubechies to produce the same amount of data.
I still don't understand how the FWT tiles the time-frequency plane, or if producing n outputs is sufficient to get the same kind of similarity measurement as an n-point circular cross-correlation using FFTs. (Technically it's a time-scale plane, not time-frequency, but I think they're the same for the complex Morlet wavelet?)  FWT is a "critical sampling" of the plane, and produces the same amount of data as FFT, so it seems fair to compare them.

The main point is that computation time is at least roughly similar for both, so I don't think you should worry about it when deciding which to use.


What acoustic properties of sounds enable us to distinguish the owner of the voice? Different people reading the same sentences, we can know from the voice who is speaking.

What makes the voice distinguishable?

Some possible perceivable differences might include:

differences in pitch frequency (bass vs. soprano voice, etc.)
differences in overtone series (formants, head resonances, accent)
variations in modulation of pitch over time (vibrato, tremolo, accent)
differences in voiced to unvoiced ratios (raspiness?)
variations in typical durations of each vowel or consonant
variations in cadence modulation (accent?)
differences in choice of phonemes for each word or phrase (accent?)


Where has this complex component come from? I'm trying to do an exercise where I have to link the images below to eachother. I'm fairly sure R1 and R4 go together but I'm not sure where the 0.25+/-j0.43 has come from.
If my reasoning is correct the 2.5 is the frequency and the 1 is the DC offset of the signal (please correct me if I'm wrong) and the equations with cosines have been found with Euler.

Thanks

At a quick glance, it looks like you're right. The impulse with an area of $1$ at zero frequency respresents a DC offset of 1. The $0.25 - j0.43$ complex value represents the amplitude and phase of the cosine at 2.5 Hz.

Specifically, the two frequency domain impulses at +/- 2.5 Hz correspond to time-domain complex exponentials:

$$
\left(0.25 - j0.43\right) \delta\left(f - 2.5\text{ Hz}\right) + \left(0.25 + j0.43\right) \delta\left(f + 2.5\text{ Hz}\right) \Leftrightarrow \\ \left(0.25 - j0.43\right)e^{j2\pi 2.5 t} + \left(0.25 + j0.43\right)e^{-j2\pi 2.5 t}
$$

(This can be shown from a table of common Fourier transform pairs, like this one)

Next, let's rewrite the $(0.25 \pm j0.43)$ factors in polar form:

$$
(0.25 \pm j0.43) = 0.50 e^{\pm j2\pi0.17}
$$

The right-hand side of the equation then becomes:

$$
0.50 e^{j2\pi0.17}e^{j2\pi 2.5 t} + 0.50 e^{-j2\pi0.17}e^{-j2\pi 2.5 t}
$$

$$
0.50 e^{j2\pi (2.5 t+0.17)} + 0.50 e^{-j2\pi (2.5 t+0.17)}
$$

We can then apply a form of Euler's identity:

$$
\cos(x) = \frac{e^{jx}+e^{-jx}}{2}
$$

$$
0.50 e^{j2\pi (2.5 t+0.17)} + 0.50 e^{-j2\pi (2.5 t+0.17)} = 1.0 \cos\left(2\pi (2.5t+0.17)\right)
$$

So, your signal should have a DC component of 1, plus a cosine at frequency 2.5 Hz with amplitude of approximately 1 and a phase offset of 0.17 cycles. From inspection, it appears that figure R1 is close to that form.


Compression algorithms specific to complex signals I am looking for (lossy or lossless) compression algorithms dedicated to complex signals. The latter could be composite data (like the left and right for stereo audio), a Fourier transformation or an intermediate step of a complex processing, an Hilbert pair of generic  signals, or complex measurements like some NMR (Nuclear magnetic resonance ) or mass spectrometry data.

A little background: I have been working on real signal and image coding, and mildly on video and geological mesh compression. I have practiced adaptive lossless and lossy compression (RLS, LPC, ADPCM), Huffman and arithmetic coding, transforms (DCT, wavelets), etc. Using standard compression tools separately on each of the real and the imaginary parts is not the main aspect of the question. I am interested in compression-related ideas specific to complex data, for instance:

sampling of the complex plane,
joint quantization of module and phase,
optimal quantization of 2D complex values: Loyd-Max is "optimal" for 1D data. I remember that 2D optimal quantization is generally more complicated. Are there 2D binnings dedicated to  complex, or for instance Gauss integers?
entropy coding methods, arranging along complex "features" (angle, magnitude),
entropy calculation for complex data,
use of specific statistical tools (e.g. from Statistical Signal Processing of Complex-Valued Data: The Theory of Improper and Noncircular Signals)
integer transformations for Gauss integers.

I have not found much tools. My scarce references are:

An experimental audio codec based on warped linear prediction of complex valued signals, 1997
Compression of Complex-Valued SAR Images, 1999  (Download)
Lossless Signal Processing with Complex Mersenne Transforms, 2003

I would be interested in more authoritative references. A compression file format dedicated to complex data would be a plus.

Complex signals are a special case of multidimensonal signals (where the dimension is two).
A lossy approach tackling compression of multidimensional signals is vector quantization. A very good resource is the book:

"Vector Quantization and Signal Compression",

co-authored by Robert M. Gray.

Vector Qquantization is a classic lossy source coding technique that is used to quantize $M$-dimensional source vectors into a set of $N$-dimensional code vectors which is referred to as the codebook, and $N<M$. A scalar quantizer is a special case with $M = 1$.

A complex signal in general can be seen as two correlated sources. Of course, if you don't want to assume any correlation between the real and imaginary part, then the answer is trivial since the two sources can be compressed independently and separately (using the standard compression techniques of real signals).

If you are interested in information-theoretic aspects of lossy/lossless compression of correlated sources, Wyner-Ziv and Slepian-Wolf coding are the keywords. These schemes are refereed to as distributed source coding schemes. However, since a complex signal can be modeled by two correlated sources, they can be considered here as well.

Finally, for analytic signals, since the imaginary/real part can be perfectly reconstructed from the other one by using the hilbert transform, the compression problem boils down to compression of a single real signal, for which there are plenty of standard schemes.

In response to the updated question:

Loyd algorithm generalized for vector quentization is referred to as LBG Algorithm (abbreviation of the name of the authors of the following paper)

Y. Linde, A. Buzo and R. Gray, "An Algorithm for Vector Quantizer Design," in IEEE Transactions on Communications, vol. 28, no. 1, pp. 84-95, Jan 1980.

It was again extended in 1989 for entropy coded vector quantization:

P. A. Chou, T. Lookabaugh, R. M. Gray, “Entropy constrained vector quantization,” IEEE Transactions on Signal Processing, vol. 37, no. 1, pp. 3142, January 1989.

I can recommend visiting www.data-compression.com/vq for a descent explanation of the scheme. It works in a simple way:

Divide the vector space.
Update the partitions.
Repeat the procedure for the new partitions until the quantization error is acceptable.

At the end, all data points that belong to a partition are quantized to the coordination of the center of partitions. The centers are indexed such that to represent all quantized values only the index is stored/sent.

Animated illustration from the same website for 2D Gaussian source:


Confused with convolutions in scipy I tried to write my own circular convolution function in python using the fact that for two signals $f$ and $g$ we have

$$ \widehat{(f * g)} = f \cdot g $$

So I tried this

from scipy import array, zeros, signal
from scipy.fftpack import fft, ifft, convolve
def conv(f, g):
  # transform f and g to frequency domain
  F = fft(f)
  G = fft(g)
  # multiply entry-wise
  C = F * G
  # transfer C to time domain
  c = ifft(C)

  return c

However, for f = array([1,2,3,4]) and g = array([5,4,3,2]) I get

conv(f,g) = [ 34.+0.j  32.+0.j  34.+0.j  40.+0.j]
convolve.convolve(f,g) = [  1.48219694e-322   2.07507571e-322   3.06320700e-322   3.26083326e-322]
signal.convolve(f,g,'same') = [14 26 40 29]

The first function is mine, the second one comes from the fft pack and the third one is from the scipy signal package. All return different results. How does this happen and what is the "correct" one?

The circular convolution I want is

$$(f*g)[n] = \sum_{p=0}^{N-1} f[p]g[n-p] $$
where $f[p] = f_{p \mod N}$. Can somebody help me?

Edit: I added a direct calculation for $(f*g)$

def direct(f,g):
    r = zeros(len(f))
    for k in range(len(f)):
        for p in range(len(f)):
            r[k] = r[k] + (f[p] * g[k-p])
    return r

which returns the same result as my own implementation above. So my implementation seems correct. However, why do convolve.convolve and signal.convolve return other results?

I do not know what convolve.convolve does but the output of signal.convolve is the linear convolution (as opposed to circular convolution) of the two sequences. And since you supplied 'same', the output has the same size as the first input vector, and the samples are taken from the center of the complete output.

The circular convolution is indeed

34   32   34   40

and the linear convolution is

5   14   26   40   29   18    8

from which you obtained samples 2 to 5 (by using the mode 'same').


Sampling $x(t)=\cos(4\pi t)+\cos(2\pi t)$ Imagine that we sample the signal $x(t)=\cos(4\pi t)+\cos(2\pi t)$ with a certain sample frequency $f_s$ and we obtain $x[n]$. Now, by ideal interpolation, we get $y(t)$ from $x[n]$.

How can we know the used sample frequency $f_s$ by looking at $y(t)?$

As a particular case, which $f_s$ verifies that $y(t)=2\cos(2\pi t)?$

How can we know the used sample frequency fs by looking at y(t)?

You can't.

The sampling theorem states that any sampling higher than twice the highest signal frequency allows for perfect reconstruction. In you case any sampling frequency higher than 4 would result in the same $y(t)$.


Bandwidth of M-ary modulation schemes I am confused with regard to how bandwidth increases with increase in number of bits transmitted for given duration.

1) Usually say in PCM system, when we increase the number of bits (n) for a given sampling frequency (fs), the bitrate increases and minimum transmission bandwidth thus increases.

Rb = nFs  and Bandwidth = Rb/2.

2)While in M-ary schemes,when the number of bits is increased the bandwidth decreases because in time domain  the signal is expanding,so obviously in frequency domain the overall width should reduce.

I am not able to understand the context in which either is true!

The required bandwidth depends on the symbol rate (or baud rate). If you have a binary system (i.e. two symbols) then the bit rate and the symbol rate are equal. If you use $M>2$ symbols, you can transmit more than 1 bit per symbol and your symbol rate, hence also the required bandwidth, is decreased. The problem with increasing the number of symbols is that you need more energy to keep the average distance between the symbols (i.e., the robustness against noise) the same. As pointed out by @vaz in a comment, the latter is only true for passband PAM systems, such as QAM or PSK.


Are there any good texts on EEG analysis geared towards programmers I just have a quick question. I'm a second year computer science major with a background in C# and C++. I'm interested in studying neuroscience after graduation, and I've been researching some subjects on my own in my free time. I'm thinking of purchasing a Neurobit Optima 4 and writing some applications for it (specifically I want to convert unspoken speech into commands an application can understand and execute), but while I'm still in the research stage, I was wondering if anyone knows about any text books or websites that discuss the analysis of EEG data or cover topics such as Discrete Wavelet Transforms that offer programming examples?

Thanks in advance to anyone who can help

Your application is not something that can be done very easily. The project would be well-suited for a PhD topic to be done over 4-5 years.

If you really want to work towards it, learn

1) Signal processing. Knowing just one technique (like the DWT) will not get you anywhere. Get a thick book on the basics.

2) Machine learning (supervised and unsupervised). Your application will need classifiers.

3) Basics of EEG. Again, get a thick book. Focus on continuous-data analysis.

4) Start reading publications in your field. Try sciencedirect.com, ieeexplore.org, google scholar, etc.

5) Offer to work as a research/teaching assistant with neuroscientists at your local university.


Is white noise WSS by nature or not? I want to know what is the difference between white noise and WSS white noise. is there any difference between them or they're equal?
and what about white Gaussian Noise?

White noise is not "WSS by nature" whatever you mean by that phrase but it can be treated as a (zero-mean) WSS process insofar as its effects in linear systems are concerned.
For example, standard linear system theory ways when the input to an LTI system is an ordinary WSS process $\{X(t)\}$ with autocorrelation function $R_X(\tau)$, then the output of the LTI system is a WSS process
$\{Y(t)\}$ with autocorrelation function $R_Y(\tau)$ given by
$$R_Y = h\star \tilde{h} \star R_X \tag{1}$$
where $h(t)$ is the impulse response of the LTI system and $\tilde{h}(t) = h(-t)$ is the time-reversed impulse response of the LTI system. The power spectral densities are related as
$$S_Y(f) = |H(f)|^2S_X(f)\tag{2}$$
where $H(f)$ is the transfer function of the LTI system. If $\{X(t)\}$ is a white noise process with autocorrelation function $K\delta(\tau)$ and we pretend that $(1)$ and $(2)$ are still applicable, we get that $\{Y(t)\}$ is a zero-mean WSS process with autocorrelation function $R_Y = K\cdot h\star \tilde{h}$ and power spectral density $S_Y(f) = K\cdot |H(f)|^2$. Of course, mathematicians would laugh at this calculation but physical experiments using the naturally occurring thermal noise in electrical circuits as a stand-in for a white noise process show that these results are pretty close to reality. As engineers, we seek equations that match the universe as we observe it (physicists seek universes that match their equations while mathematicians don't care) and so we go blithely on our way treating white noise as a WSS process in linear systems and everything works out OK. The troubles start when we start treating white noise as a WSS process in nonlinear systems and the world comes crashing down about our ears and we need to start paying attention to what the math people are saying.
White noise is referred as white Gaussian noise if we pretend or claim or insist that $\{Y(t)\}$ is a Gaussian process which means that not only are all the random variables $Y(t)$ Gaussian random variables, but every finite set $\{Y(t_1), Y(t_2), \cdots, Y(t_n)\}, n \geq 2,$ of random variables is a set of jointly Gaussian random variables. Standard random process theory says that when a Gaussian process is passed through an LTI system, the output is a Gaussian process but this fact does not allow us to reverse-engineer the result and claim that all the $X(t)$'s are also Gaussian random variables.


What would the behavior be of an FFT without the phase (imaginary part) If I have a sound recording with a sinusodial frequency with a constant amplitude and I do an FFT on the time domain of this recording to get the frequency back:

If I were to ignore the imaginary part and just watched the amplitude (real part) would there be anything different about it? Would it also be constant?

I'm trying to do an FFT and I notice that the frequency spike oscillates up and down as if it were a sine wave, I have ignored the imaginary part on my FFT.

The real part isn't the amplitude.

The imaginary part is the sine wave part of the FFT.  Thus exactly periodic- in-aperture sine waves (with a phase of zero at the window start) won't show up on the real part, and the FFT real part will see a zero signal.  Only sinusoids with a periodic cosine wave (with a phase of zero at the window start) component will show up at all.

The "beating" you are seeing is your signal alternating between being in and out of phase compared to a cosine basis vector.


Clarification Regarding Changing the Variable Assume I have that :
$$
X(j\omega)=\begin{cases}e^{-2\omega}&\text{if $|\omega|<2\pi$}\\0&\text{if otherwise}\end{cases}
$$
Say, I want to make this change substitution $X(j(\omega-\Omega))$. Would this imply :
$$
X(j(\omega-\Omega))=\begin{cases}e^{-2(\omega-j\Omega)}&\text{if $-j(2\pi+\omega)<\Omega<-j(\omega-2\pi)$}\\0&\text{if otherwise}\end{cases}\;\;???
$$

The implication seems a little bit wrong, but that is okay.
When dealing with the functions in the frequency domain, let’s discuss it over the function $X(j\omega)$, when shifting or scaling the function on the frequency axis, the operations are done on the angular frequency $\omega$ or the frequency $f$.
On the frequency axis, there is no imaginary number, that is, the frequency axis is full of real numbers. Therefore, calling a point on the frequency axis as e.g. $j2$ Hz is incorrect.
When I do calculations on the paper, I prefer using $\omega$ instead of $j\omega$ when defining the variable of the Fourier transform of a function e.g. $X(\omega)$. However, in the case of transferring a Fourier transform into the Laplace transform under specific conditions, $j\omega$ notation aids in preventing confusion of where to put the $s$ variable as in this process, we do assignment $s = j\omega$.
I've changed the given function a little bit i.e. I've taken the absolute value of the angular frequency $\omega$ in the power of the exponential in order to illustrate what happens better.
Okay, along with these knowledge, let’s focus on your function. The amplitude spectrum of the function $X(j\omega)$ is given below:

As it is noted, the mathematical explanation of the function is;
$$X(j\omega)=\begin{cases}e^{-2|\omega|},&\text{if $|\omega|<2\pi$}\\0,&\text{otherwise}\end{cases}$$
So, let’s shift the function $X(j\omega)$ by $\Omega$ to the right on the frequency axis as I’ve assumed $\Omega>0$. This shifting process, which is a frequency shifting, is done in the same way as time shifting. The resultant shifted function $X(j(\omega - \Omega))$ will be in the figure below:

So, the resultant mathematical expression becomes;
$$X(j(\omega - \Omega))=\begin{cases}e^{-2(|{\omega - \Omega}|)},&\text{if $|{\omega - \Omega}|<2\pi$}\\0,&\text{otherwise}\end{cases}$$
The expression $|{\omega - \Omega}|<2\pi$ can be represented comprehensively as ${{-2\pi}+{\Omega}}<{\omega}<{{2\pi}+{\Omega}}$.


Combining two images I am working with speckle noise reduction in ultrasound images. I have used wavelet transform for removing the noise  by eliminating certain frequencies in order to eliminate any existing noise. Since in an image HH,LH,HL components contains most of the noise. We can eliminate noise by eliminating those components.

  step 1 : read an input image.
  step 2 : apply discrete wavelet transform.
  step 3 : Eliminate LH, HL, HH  by making these components zero.
  step 4 : Take inverse wavelet transform(let us consider this as output 1).
  step 5 : read the same input image and again take discrete wavelet transform.
  step 6 : eliminate LL component by making zero and keep LH,HL & HH because it contains edge information. (trying to remove noise from the edge information & add it back to Output 1)
  step 7 : take inverse wavelet transform .
  step 8 : Applying a Hybrid Filter (Median- Weiner) to the above image
  step 9 : We get the filtered image ( let us consider this as output 2)

I have represented the same in the below figure.

Now, I need to combine output 1 and output 2 . I just added both the outputs but am not getting proper output . Please suggest me the technique which I should use to combine these two output images. Kindly tell me whether it is possible or not and also tell me whether this concept is correct or wrong.

You can implement discrete wavelet transforms on output1 and output2. Use LL from output1 and LH, HL, and HH to form a new block, then implement inverse wavelet transform.


Nyquist criterion for zero isi I am working with the Nyquist criterion for zero ISI when I find something that makes me think that I misunderstood this criterion.

What Nyquist says is that to avoid ISI in the sampling we must ensure that the sample frecuency is less than 2 times the bandwith of the signal, so if the channel filter is something like a rise cosine the ISI in the sampling will be $0$.

At the same time, we have (another) Nyquist criterion which says that to be able to perfect reconstruct a signal, the sample frequency must be at least $2$ times the bandwith of the signal.

So... taken into account that the case of $F_s=2B$ is not implementable, how can we assure no ISI at the same time that we have perfect reconstruction? Aren't this two criteria contradictory? Am I missing something?

Thanks!!!

My answer is related to this question
Why root raised cosine filter can eliminate intersymbol interference (ISI) ?

A classical system with Nyquist pulse $p(t)$ is:

The equivalent baseband signal has bandwidth limited in $[-1/2T,1/2T]$ and a sampling rate $F_s > 1/T$ is enough to avoid aliasing.

We normally sample at rate $t=kT$. I wonder what would happen if suddenly I sample faster than the normal sampling rate, for example at $t = kT/2$ instead of $t = kT$ ?

In the frequency domain, the composite filter $g(t) = p(t) \star p^*(-t)$ has Fourier transform $G(f)$ satisfies $\sum_k G(f-k/T) = T, \forall f$. Double the sampling rate $\sum_k G(f-2k/T)$ is not constant any more. This is true at least for typical pulses such as sinc or raised cosine.
In the time domain, thank Dan Boschen for the visualization below, the sample at impair $k$ of $kT/2$ is contributions of two adjacent symbols.

Thus, is it true that for free-ISI-Nyquist-pulse-shaped signal, the largest sampling rate is the value of nominal (Nyquist) bandwidth $1/T$ ?

One can say that we still have correct data by taking only one per two samples (downsampling). What if I sample at rate $F_s > 1/T$ that satisfies Shannon-Nyquist sampling theorem to avoid aliasing ? We still have ISI.

OK another one can say that the criterion $F_s > 1/T$ is in the sense that the entire information after being sampled is intact, just take a baseband filter to get back the continous time spectrum and then sample it at the correct sampling rate at $t=kT$. Have I answered my questions by myself ?


What type of sensors have complex Gaussian noise? What type of sensors have complex Gaussian noise?

Do strain gauges also exhibit such noise?

No physical system suffers from, or produces, complex Gaussian noise. Quadrature systems can be modelled as being complex, with the in-phase branch corresponding to a real value, and the quadrature branch corresponding to an imaginary value. If each branch presents Gaussian (thermal) noise, then the noise can be modelled as complex Gaussian, too. However, this is just a modelling technique: the actual physical noise is never complex.

To answer your first question: sensors that operate in quadrature, and are modelled using the complex envelope, can be modelled as having complex noise.


Introductory book on time-frequency analysis? I'm looking for an introductory book to time-frequency analysis. The book should be practical in nature and not mathematics heavy. Suggestions?

I recommend this book:

Understanding Digital Signal Processing, Richard G. Lyons

This book explains the basics concepts of digital signal processing, which includes time-frequency analysis, in a very intuitive way.


Reconstruction using sinc The signal that I produce above. What is the reason for it to slide to the right? In oversampling at Nyquist rate can I make like below picture ? do you think the signal i produced at nyquist rate is correct, why so on the right? I was expecting to draw as below picture.
(Blue is the my signal, Red ones sampling and Yellow is the reconstruction using sinc.)
If needed I can share my all code.
Thank you all.

After seeing the code, I think the mistake may be

sinc((t(i)-n*Ts1)/Ts1)

Change it to

sinc((t(i)-(n-1)*Ts1)/Ts1)

An extra Ts1 delay was being applied to sinc. The index $n$ starts from 1 for MATLAB but it should start from 0.


Partitioning the energy of a signal between its components In physics, the total energy of a system is the summed energy of its component subsystems, and energy must be conserved, ie, adding the energy of the components gives the total energy.
In signal analysis (discrete time signals), one would expect that the energy of a signal would be defined in such a way that energy would again be conserved.  This would make it consistent with normal usage in physics.
There are many types of component analysis. In general, there is no requirement that the components are zero-mean or orthogonal, and in data-dependent component extractions, the components are not even predictable, but in all cases it is expected that they linearly sum to the original signal values.
This is incompatible with the usual definition of signal energy as the sum of the squared component values, since the component energies do not sum to the original signal energy.
This is confusing.  Is there another approach to defining the component energies that conserves energy but does not destroy the requirement that the linear sum of components must be equal to the original signal?

Say you decompose the signal $x(t)$ into two components: $$s(t) = x(t) + y(t).$$ As you say, there are myriad ways to do this. The energy of $s(t)$ is $$E_s = \int_{-\infty}^\infty s^2(t) dt.$$
If we calculate it using the decomposition, we obtain:
\begin{align}E_s &= \int_{-\infty}^\infty (x(t) + y(t))^2(t) dt \\
                 &= \int_{-\infty}^\infty x^2(t) dt + \int_{-\infty}^\infty y^2(t) dt + 2\int_{-\infty}^\infty x(t)y(t) dt
\end{align}
Clearly, if you require that $E_s = \int_{-\infty}^\infty x^2(t) dt + \int_{-\infty}^\infty y^2(t) dt$, then it is necessary that $$\int_{-\infty}^\infty x(t)y(t) dt = 0$$ or, in other words, $x(t)$ and $y(t)$ must be orthogonal.
By the way, I have never seen it claimed that all existing decompositions have the energy-conserving property you are looking for, not even in the context of physics rather than engineering.