Histogram Equalization - mapping a certain portion of pixels to pure black and white How can I make sure to map, for example, 7% of the pixels to pure black (0) and pure white (255) in a grayscale image?

Compute the image histogram, then its running integral (i.e. an array T with as many entries as there are gray levels, such as T(k) = sum_{i=0, 1..k} (histogram(i)).

The last element of the integral is, of course, the area of the image.
For the lower threshold, start from the left and find the level at which T exceeds your desired fraction. Analogously, for the high threshold, start from the right and find the level at which T goes below your high fraction.


How to "draw" the function/wave to send symbols using QAM? Say you have the following 8QAM constellation diagram:

Say, you want to send the sequence:

001 110 101 000 001

How would you draw the function/wave containing the symbols?

My problem in understanding this is:

I know you use amplitude and phase shifting for defining symbols. Each symbol has a complex function like this:

$$g(t) = A \times e^{i(ft + \phi)}$$

($A$ = amplitude, $f$ = frequency, $\phi$ = phase)

You have one function for each of the possible symbols. But how do you draw the sequence of symbols (or even one symbol)? Isn't $e^{ix} = \mathrm{cos}(x) + i \times \mathrm{sin}(x)$ ? Do I draw the real part or the imaginary part?

Thanks!

The baseband QAM signal is complex, and the only way to draw it is by doing two drawings, one for the in-phase (real) component, and one for the quadrature (imaginary) component.

The passband QAM signal, though, is real, and it is a pulse-shaped carrier whose amplitude and phase depend on each symbol. Myself, I would "draw" it using Matlab or some other numerical program. I don't see the point in trying to draw it by hand.

For example, in Matlab, assuming QPSK and transmitting three symbols with a rectangular pulse shape, at a rate of one symbol per second and a carrier of 5 Hz:

t = 0:0.01:2.99;
p = ones(1,100);         % pulse shape
s = [1+1i, -1-1i, 1-1i]; % symbols to transmit
u(1:100:201) = s;
qam_bb = conv(u,p);      % pulse shaping to generate baseband signal
subplot(3,1,1);
plot(t,real(qam_bb)); title('Baseband In-Phase'); axis([0 3 -1.3 1.3]);
subplot(3,1,2);
plot(t,imag(qam_bb)); title('Baseband Quadrature'); axis([0 3 -1.3 1.3]);
qpsk = real(qam_bb.*exp(2*pi*5*t*1i));  % upconvert baseband signal
subplot(3,1,3);
plot(t,qpsk); title ('Passband QPSK');

produces this plot:


Memory requirements for FFT I have a system with limited memory. So basically the question is this: If I have (say) 24KB of data space and 4096 data samples of single precision floats (16KB) can I do an FFT on the data? More generally, what is the minimum data space required to do an FFT?

From the documentation of the CMSIS-DSP FFT function:

  The FFT functions operate in-place. That is, the array holding the
  input data will also be used to hold the corresponding result. The
  input data is complex and contains 2*fftLen interleaved values as
  shown below.

It is possible to implement the FFT in such a way, that it does not require additional buffers but uses the memory of the input array for its output. This method is usually referred to as "in-place" calculation. From the above citation I think that the CMSIS-DSP library does exactly that (though we can not know for sure without the source code).

So in this case, yes, a memory size of 16 KB is sufficient.

That an in-place calculation is possible can be understood from the butterfly diagram, see for example here.


Does oversampling improve processing gain? I am used to Time-Bandwidth product describing the processing gain or improvement in SNR for a spread spectrum signal using a matched filter.

One could view the bandwidth $B$ as the sample rate (complex) and the time $T$ is the integration time and thus there are $TB$ samples.  Thus, the matched filter produces a gain of $10\log_{10}N$, with $N$ being the number of samples.

So with all of this in mind, suppose a bandpass filter removes all noise power except for the signal bandwidth, thus the noise power is $kTB$.  The signal power is the amplitude $A$ squared, $A^2$.  If we over sample the signal by some factor $M$, then I want to say that the processing gain is now $MTB$ because there are $M$ more samples.

I also have worked some theory of this as well as a MATLAB script and all is indicating that it is true.  Am I missing something?  Why not over sample signals to provide more processing gain using a longer matched filter?

You need to define what you mean by SNR and processing gain, because Marcus Muller and Qasim Chaudhari can both be considered correct depending on your definitions of SNR and processing gain. A processing gain is usually taken to be the SNR at the output of a system divided by the SNR at the input of a system. If you consider the SNR to be the signal power within the band of interest divided by the noise power within the band of interest, then there is no processing gain as you increase the sampling rate. This assumes an ideal (sinc) filter. In fact, without a perfect matched filter and perfect stopband attenuation (equal to zero), there is a processing loss caused by the out-of-band noise aliasing back into the spectrum after downsampling to the symbol rate. Although in all practical cases, this is very small and not really worth talking about. In other words, there is a fundamental limit to the SNR you can achieve by oversampling in this case. Without oversampling, it is likely that you will have poor filtering, and out-of-band noise may leak into the spectrum of interest (degrading it).

Now consider the SNR to be defined as the signal power divided by the total noise power contained in the samples (i.e. across the larger Nyquist bandwidth). In this case, you can reduce the noise level by filtering the excess noise that was added out-of-band. In other words, you will get a processing gain that is equal to the sampling rate divided by the true bandwidth of your signal (one over the baud interval for single carrier pulse shaped systems). This is often called a processing gain in ADC literature for example. Notice that if in your original samples, the matched filter design is good enough to reject out-of-band noise to negligible levels then there is no benefit in letting more noise into the system simply to reject it in order to obtain a high processing gain. It is probably worth noting that oversampling does result in a higher order matched filter which typically would result in a better stopband attenuation (making the excess noise less pronounced). This effect is very subtle and will never allow you to achieve a better SNR than you could by using a perfect matched filter (stopband attenuation = 0) at a lower sampling rate.


Periodicity property of DFT in time and frequency domain I am trying to understand the periodicity of the DFT. How can this property (both in time and frequency domain) be used and can be helpful while developing on a DSP?

It would be good to see some source code or pseudo code, in MATLAB preferably, where this property is exploited/demonstrated.

If the expression that defines the DFT is evaluated for all integers
  $k$ instead of just for $k = 0, \dots, N-1$ , then the resulting
  infinite sequence is a periodic extension of the DFT, periodic with
  period $N$.

  The periodicity can be shown directly from the definition:

  $$ X_{k+N} \stackrel{\mathrm{def}}{=} \ \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} (k+N) n} = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} k n} \underbrace{e^{-2 \pi i n}}_{1} = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} k n} = X_k.$$

  Similarly, it can be shown that the IDFT formula leads to a periodic
  extension.

  Source: Wikipedia.

While the periodic property of DFT isn't very widely utilized, it often causes aliasing problems.

Source: http://www.dspguide.com/ch10/3.htm


Objeckt tracking - filter training (MOSSE) For a project I'd like to set up an easy to implement tracker. I read several papers and chose to go for the mosse (minimum output of squared error) tracker from Bolme et al.

The  formula to obtain the mosse filter is (equation 5 in the paper):

H* = sum(Gi x Fi*)/(sum(Fi x Fi*) +e)

Where Gi is the fourier transform of the desired output, a 2D gaussian shaped peak centered on the target in the training image.
Fi is the fourier transform of the training image.
x denotes elementwise multiplication.
the asterisk indicates that it is the conjugate of the fourier transform.
e = 0.1, is a regularisation parameter that induces noise for enhanced noise tolerance.

To generate training images I took the first frame and applied small affine transformations to it.

But somewhere I made a mistake. Since my filter looks nothing like the one in the paper.

mine looks like this.
on the left the input image and on the right the obtained filter.

I even checked if the transformation to the training images and the desire output are faulty. However, they look fine to me.

Unfortunately I have no background in signal processing or linear algebra.
So I would be really grateful if you could help me.

Below I added the Matlab code for my failed filter.

% get variables to store the numerator and denominator
sum_numerator = zeros(size(filter));
sum_denominator = zeros(size(filter));

for a = 1:numberOfTrainingImages
    Fi = fft2(image);
    Gi = fft2(desiredOutput);

    % sum_numerator is the correlation between input and desired output
    sum_numerator = sum_numerator+ ( Gi .* conj(Fi));

    % sum_denominator is the energy spectrum of the input
    sum_denominator = sum_denominator + ( Fi .* conj(Fi));

    % Generate a training image and a according desired output
    [image, desiredOutput] = generateTrainingImage (OriginalImage);
end

% elementwise division to get the filter
H = sum_numerator ./ sum_denominator;

The solution to the problem was, to rearrange the quadrants of the filter. After rearranging the quadrants, my filter looks like this:

A lot more like the one in the paper. However, it still lacks in tracking precission.   To tweak this I played arround with the 2D gaussian. And obtained the filter below. Still not optimal but it's an improvement.

This is the formula for the 2D gaussian
$$
g_{i}(x,y) = e^{-\frac{(x-x_{i})^2+(y-y_{i})^2}{\sigma ^{2}}}
$$

In the paper I mentioned in the question, Bolme et al. used $\sigma = 2$. Instead of this I obtained $\sigma$ now dependent on the size of my template

$$
\sigma = \sqrt{{imgSize_{x} \cdot {imgSize_{y}}}} \cdot \frac{1}{16}
$$

like used by Galoogahi et al. (2015). Correlation Filters with Limited Boundaries


How does a correlation receiver ensure which transmitted data does a user see? Hello! I understand that using correlation, one can set up a system wherein multiple users can transmit at the same time using the same channel. And theoretically, I understand that if the correlation is equal to 0, then one can restrict which one will be able to receive which kind of data. For example, If User 1 transmits data X and User 2 transmits data Y, add them, and send them in the same channel, then User 3 can receive data X provided that it knows the correlation. And another user cannot receive data X.
What I don't quite understand is the mathematics behind it. How do I prove that data X transmitted by User 1 gets received by User 3, and not another User. What equations can I use, or how does this use the autocorrelation function $$ R_{xy}\left[ k \right] =\sum_{n=-\infty }^{+\infty }x\left[ n-k \right]y*\left[ n \right] $$
Thanks!

Correlation is a sum of products. With a "correlation function" (such as autocorrelation function, cross-correlation function, etc) we are actually doing multiple correlations, one for each possible change in some variable of interest. In this case, the OP has shown a cross-correlation function since it is the correlation between $x$ and $y$ with the variable of interest being a delay index $k$. It is also very important to point out that it is a complex conjugate multiplication (the OP did not raise the * so may not be clear so rewriting here):
$$R_{xy}[k] = \sum_{n=-\infty}^{+\infty} x[n-k]y^*[n]$$
In this case the variable of interest is the delay index $k$ where the correlation (sum of products) is being computed for each possible delay index $k$: we delay one waveform relative to another and repeat the sum of products computation (lots of processing!), and here we would be looking for either a delay index $k$ where the correlation is strongest to the waveform we are looking for (so if $y[n] = x[n] + noise$ for example where "noise" could be other users), and we would desire that other users be uncorrelated so not have any correlation response.
We could also use this property for timing and synchronization by using codes that have an auto-correlation property where they are very strong when the codes align and nearly zero for every other offset. Auto-correlation is given as:
$$R_{yy}[k] = \sum_{n=-\infty}^{+\infty} y[n-k]y^*[n]$$
A simple example of this is a 11-bit barker code given as 1 0 1 1 0 1 1 1 0 0 0
If you replace all zeros with -1 in the above sequence, you will see that the autocorrelation is 11 when the code aligns and -1 for every other rotational shift (by doing a sum of products which is correlation).
In terms of using this with multiple users (CDMA), consider that each user has a unique code, and very importantly we'll assume that the power from each user that a receiver sees is about the same (this is very important for CDMA to work properly and thus a reverse power control mechanism is often used, or as it the case of Satcom such as GPS, the transmitters are so far away that all receivers receive about the same power level from each of the satellites). If each user has a unique code and those codes for each user are known to be uncorrelated, we can use cross correlation of correlating the received signal given by $r[n]$ with the code for a given user given by $y[n]$ to extract the data for that user (one way to do this such as with spread spectrum is to modulate the code sequence so that the entire sequence is inverted according to the data and thus the data is at a much lower rate than the code). The math is to just use the cross-correlation function, or more likely use a delay lock loop to align the local code with the code position in the received signal and just keep computing the correlation: if our local code is known to be aligned from a prior acquisition operation, then we simply multiply the received signal by the local code and accumulate over the sequence (symbol) period. This is an integrate and dump receiver. Often there will be additional integrate and dump channels to maintain tracking (early-late).


Do you know this phenomenon in DFT? Please let me ask you about a phenomenon in DFT. Below DFT program outputs excellent results for low sample rates, but for high results are very bad.

Specifically: for sample rate 50 samples/sec, segment 1000 samples, existing frequency 36.7634 Hz, grid from 35.2364 Hz to 40.7896 divided by 40, outputs maximum magnitude squared 1.22629e+15 corresponding to testing frequency 36.7635.

But for sample rate 50,000 samples/sec, segment also 1000 samples, and same existing frequency and grid, outputs maximum magnitude squared 8.46881e+14 corresponding to testing frequency 38.2907 Hz. Error is 5 significant digits.

Are you aware of this phenomenon? If yes please explain me that.

// Grid is created by variables lowTestFreq and highTestFreq
// both of double type, and integer type M.

#include <iostream>
#include "math.h"

int main()
{
//Unchanged variables.
constexpr double pi = 3.141592653589793; //pure number.
int n = 0; // sample.
double xn = 0; //pure number.
double Re = 0; //pure number.
double Im = 0; //pure number.
double xn_cos = 0; //pure number.
double xn_sin = 0; //pure number.
double PowerSpectrum =0; //pure number. Name given by my dear dsp Teacher Mr. Lyons.

//Sampling frequency and samples of segment.
int samplFreq = 50000; // sample/second.
int N = 1000; //sample.

// 1st existing sinusoidal.
double ampl_1 = 70000; //pure number.
double existFreq_1 = 36.7634; // Hz.
double phase_1 = 0.67 * pi; // rad.
double unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

// 2nd existing sinusoidal.
double ampl_2 = 60000; // pure number.
double existFreq_2 = 4505.75; // Hz.
double phase_2 = -0.37 * pi; // rad.
double unitAngle_2 = (existFreq_2 / samplFreq) * (2 * pi); // rad/sample.

// Testing grid.
double testFreq = 0; // Hz.
double testUnitAngle = 0; // rad/sample.
double testAngle = 0; // rad.
double lowTestFreq = 35.2364; //Hz.
double highTestFreq = 40.7896; //Hz.
int M = 40; // pure number.
int m = 0; // pure number.
double testFreqStep = (highTestFreq - lowTestFreq) / M; // Hz.

for (m = 0; m <= M; ++m)
{
testFreq = lowTestFreq + m * testFreqStep; // Hz.
testUnitAngle = (testFreq / samplFreq) * (2 * pi); // rad/sample.

Re = 0; Im = 0;

for (n = 1; n <= N; ++n)
{
xn = ampl_1 * sin(n * unitAngle_1 + phase_1)
   + ampl_2 * sin(n * unitAngle_2 + phase_2);

testAngle = n * testUnitAngle;
xn_cos = xn *  cos(testAngle);
xn_sin = xn * -sin(testAngle);

Re += xn_cos;
Im += xn_sin;
}
 std::cout << "Test Freq " << testFreq; // Hz.
 PowerSpectrum = Re * Re + Im * Im; //pure number.
 std::cout << " MagnSquared " << PowerSpectrum << std::endl;
}
return 0;
}

If you are sampling 50 samples/sec your nyquist frequency is only 25Hz, so your test signal is being aliased. You should expect different results if your sampling rate is less than double the highest frequency in your signal.


Why the plots of spline and cubic interpolation are exactly same? I am trying to watch difference between cubic interpolation and spline interpolation using matlab plot but i am getting same plots in both cases using interp1 command
My code is below. In below code ,if i use spline in place of cubic, i still exactly get same plot,why?
clc
clear all
close all
x = 0:pi/4:2*pi;
v = sin(x);
xq = 0:pi/16:2*pi;
vq2 = interp1(x,v,xq,'cubic');
plot(x,v,'o',xq,vq2,':.');

The best way to understand cubic splines (under whatever name or how specified) is to derive and code them yourself.
What you are looking for is a parametric equation with these properties:
$$
\begin{aligned}
f(0) &= y[n]    = y_0\\
f(1) &= y[n+1]  = y_1\\
f(2) &= y[n+2]  = y_2\\
f(3) &= y[n+2]  = y_3\\
\end{aligned}
$$
Your candidate function is a cubic polynomial, of course.
$$ f(t) = a t^3 + b t^2 + c t + d $$
Plug 'em in.  You get four equations with four unknowns ($a,b,c,d$)
Once you have solved the equation, use it to interpolate between the middle two points.
Extra credit:
Do the vector equivalent in N-space with the conditions given as two points in space and the velocity vectors at each point with one unit of time to get there.
Then code your solutions.  Print the results so you can see some precision.  And compare your answer to Matlab's.
Then post your results here.


Morphology - Why does the concept of opening and closing images exist when we already have erosion and dilation? If we consider the closing operation of an image, we first dilate the image and then erode it. So, first we expand the foreground (white pixels) and and then we take away some part of the foreground. The result of this particular operation is similar to dilation, but less destructive of the original shape boundary.

However, we can combine erosion and dilation in different ways to control the extent of increase or decrease of the shape boundary. Then, why do we specifically define closing and opening operations?

As you said closing is completely defined by a dilatation followed by an erosion with the same structuring element (SE) (respectively opening is an erosion followed by a dilation with the same structuring element).

The key point is that is must be the same structuring element and this is much more specific than an erosion followed by a dilatation with different random SE.

Closing and opening operations are defined just for convenience, because giving an image and a SE the result can easily be predicted and interpreted by the user.

  The result of this particular operation is similar to dilation [...]

Not really, dilatation is more like a thickening operation while closing fills holes and smooth borders accordingly to the SE (opening with a disk rounds angular borders for instance).

Try erosion, dilatation, opening and closing on simple shapes (square, disk, triangle, hollow shapes, set of close shapes) with simple SE (share, disk, lines) to better visualize the effect of each of them. you can also try erosion followed by dilatation with different SE to see the difference with an opening.

To sum up opening and closing are defined for convenience only, as it is a specific case of erosion and dilatation composition that have noteworthy properties.


DWT with more precise scale With a signal x of length $2^{16}$,  the DWT (computed for example with pyWavelet's wavedec) has only 16 or 17 rows (the wavelet scale is dyadic).
Is there a DWT with more scales? (ex : 1024 rows instead of 16!)

You will need to use a continues wavelet transform (cwt) where you can specify the resolution of the scales (as done in the image you posted mathworks.fr/help/releases/R2013b/wavelet/ref/cwt.gif). The DWT only has $J$ scales by definition for $2^J$ signals.


Finding audio signal time delay: cross-correlation vs. impulse response? I need to find a time delay between the two audio signals, one comming from one set of speakers, and the other being the source signal that is initialy sent to speakers. I send some noise through those speakers and they come to microphone with a time delay compared to a source signal.

Cross-correlation method
As I understand, there is a cross-correlation method where I send both recorded signals to frequency domain via Discrete Fourier Transform (DFT), then multiply the frequencies, and finaly go back to time domain via Inverse DFT. There should be a peak in time domain telling me what the delay is (let's disregard the accuracy for now.)

IR method
There is also an Impulse Response method, where I also take two signals to freq domain via DFT, then divide the two to get a system transfer function, and finaly do an Inverse DFT of the Transfer function; this will give me the Impulse Response, which should also have a strong peak at the delay point.

My question is:
Am I right in the short descriptions above? The main difference between the two methods seems to be that frequencies are multiplied in the cross-correletion method, and dividied in the IR method? Have I missed something important?

Your description of the cross-correlation method is not entirely correct. What you describe is actually a convolution. For a cross-correlation you would perform an FFT on both signals and then multiply one with the complex conjugate of the other, and finally perform an IFFT.


Fourier Series Proof I want to ask Question about the Fourier series in continuous time domain while reading a book signals and systems Alan Oppenheim. I have confusion in understanding the statement on page 189 of its 2nd edition.
To derive alternative form of Fourier series, we first rearrange the summation in Eq.3.25 as
$$x(t) = \sum^{+\infty}_{k\ =\ 1} [a_k e^{jk\omega_0t}+a_{-k} e^{-jk\omega_0t}]$$

where Eq. 3.25$$x(t) = \sum^{+\infty}_{k\ =\ -\infty} a_k e^{jk\omega_0t}$$
I want to know two things
1) How summation from -infinity to +infinity changes to 1 to +infinity

2) How we get this term $$a_{-k} e^{-jk\omega_0t}$$ in the equation.

You just need to split the sum into the positive and the negative indices (without forgetting the index $k=0$):

$$x(t)=a_0+\sum_{k=1}^{\infty}a_ke^{jk\omega_0 t}+\sum_{k=-\infty}^{-1}a_ke^{jk\omega_0 t}\tag{1}$$

By simply changing the sign of the index $k$ in the second sum in $(1)$ you get

$$\begin{align}x(t)&=a_0+\sum_{k=1}^{\infty}a_ke^{jk\omega_0 t}+\sum_{k=1}^{\infty}a_{-k}e^{-jk\omega_0 t}\\&=
a_0+\sum_{k=1}^{\infty}\left[a_ke^{jk\omega_0 t}+a_{-k}e^{-jk\omega_0 t}\right]\end{align}\tag{2}$$


Baseline correction for continuous EEG signals I'm new to the topic so forgive the trivial question.

I am trying to understand data preprocessing for EEG signals. As I understand, after recording continuous EEG signals, they should be band-pass filtered and segmented into epochs (so the window of time that is relevant to the experiment). After that, these epochs should be baseline corrected. Can anyone explain what this concept of baseline correction means?

A baseline is generally an unwanted artifact that alter interpretation or computation on signals. The concept is not totally well-defined, and several different names coexist, with similar meaning:  drift, continuum, trend, background.

It is often described as a large-scale, slowly-varying reference level, that should be removed to compare, or adjust the signal with respect to a globally flat or  zero-reference. It can be seen as a generalized (non-constant) offset or bias  correction. When removed, it is generally simpler to further process the signal, quantify parameters, extract features.

In some domains, this smooth reference may also  include random noise.

Similar problems are answered here:

How to remove or filter the drift problem in measured Strain signal?
Help with automated baseline removal


Segmentation of small artifacts - Image Processing I am working on a task of segmenting small artifacts from a set of images like below:

I am looking for a generic approach, but unsure of what to search?

I tried looking at density flow analysis, but don't know if that fit in.

I am unable to come up with a single pattern in all the images. For example I look out for "loops", but few images have loop plus more lines, like in (d) above.

One property I am certain about, the percentage of in these artifacts will be very less compared to the whole image.

Can you help me with:

Any algorithm you know which will solve the problem?
Any research paper, where a similar problem is solved. It is okay if that is not the topic of research paper, as this kind of work is usually present in per-processing stages.

If this is really what your input looks like then you should be able to handle this with the standard morphological image processing operations.

Assuming foreground as white and background as black, it basically comes down to a morphological open with a big enough structure element to elide the tails. The only trouble is that doing so will not work given that the input has holes so large the lines in the body of the shape itself would also get nuked by the open operation. One way to handle this is to perform the open on the figure without holes and then put the holes back in.

For example, in C++ using OpenCV:

void Test()
{

    Mat src_bw;
    cv::threshold(
        imread("some_file.png", CV_LOAD_IMAGE_GRAYSCALE),
        src_bw, 0, 255,
        CV_THRESH_BINARY | CV_THRESH_OTSU
    );

    imshow("src_bw", src_bw);

    // perform contour extraction using RETR_CCOMP i.e. we are assuming input doesn't have nested islands/holes just simple holes.
    std::vector<std::vector<Point> > contours;
    std::vector<Vec4i> hierarchy;
    findContours(src_bw.clone(), contours, hierarchy, RETR_CCOMP, CV_CHAIN_APPROX_SIMPLE, Point(0, 0));

    // iterate over the contours, draw the top-level contours into "blobs"
    Mat blobs = Mat::zeros(src_bw.size(), CV_8U);
    for (int i = 0; i < contours.size(); i++)
        if (hierarchy[i][3] == -1)
            drawContours(blobs, contours, i, 255, CV_FILLED, 8, hierarchy, 0, Point());

    imshow("blobs", blobs);

    // do an 'open' on blobs which will remove the artifacts.
    auto shape = MORPH_ELLIPSE;
    Mat strel = getStructuringElement(MORPH_ELLIPSE, Size(5,5));
    Mat blobs_opened;
    morphologyEx(blobs, blobs_opened, MORPH_OPEN, strel);
    imshow("blobs_opened", blobs_opened);

    // draw the holes back into the output
    Mat output = blobs_opened.clone();
    for (int i = 0; i < contours.size(); i++)
        if (hierarchy[i][3] != -1) {
            drawContours(output, contours, i, 0, CV_FILLED, 8, hierarchy, 0, Point());
        }

    imshow("output", output);
}

yielding the following:


What is the duration of the spreading sequence I have a question regarding the spreading sequence duration using Walsh matrix. I've read a paper related to that, they consider the signal $S = R(S) + I(S)$, where $R(S), I(S)$ are the real and imaginary parts of the signal $S$ resulted from the modulation $QAM$.

what is done, briefly, a Walsh-Hadmard matrix is built with dimension $4$x$4$, then the real are imaginary parts are multiplied by the spreading code taken from the Walsh matrix and then transmitted into receiver. So the resulted signal before transmitting is $S = W_1R(S) + W_2I(S)$ where $W_1, W_2$ are two different spreading code taken from the Walsh matrix.

So, without talking about why that process was done, My question is related into the signal after spreading itself:

Will the duration of transmitted signal "after spreading" is similar to duration of signal $S$ before spreading?I think yes, because the duration of chip of spreading sequence is shorter by $n$ times of duration of signal $S$, it's similar to this question here.. If that right, what's the downside of using spreading sequence? I mean when using it, what's the disadvantages we face compared with any other technique. Finally, does using spreading sequence as above process mean that we are using CDMA? ..

thank you

When you multiply the real and imaginary parts by $\mathbf{W}_1$ and $\mathbf{W}_2$, respectively, the resulting signal is no longer $S$. It will be a matrix.

Regarding your questions:

1- Yes, the signal after spreading has the same duration.

2- One of the disadvantages of Walsh-Hadamard codes is that they don't have clear auto-correlation peak, which makes synchronization at the receiver side challenging. Also, they must have a dimension of $2^k$, i.e., you cannot have arbitrary dimension to support an arbitrary number of users. In your example, to make the system efficient, you need 4 users, otherwise, one user will occupy 4 times the minimum required bandwidth, which isn't spectrally efficient.

3- Yes. This is called direct sequence (DS) CDMA.


Plotting the spectrum of a continuous sinusoid function I am just learning about spectrums and I am a little confused. I have given a simple example and where I get stuck.

  Plot the spectrum of  $x(t) = 5\cos\left(800\pi t + \frac{3 \pi}{7}\right)$.

So I use reverse-Euler's to write this as
$$x(t) = \frac{5 e^{i \frac{3 \pi}{7}}}{2} e^{i 800\pi t} + \frac{5 e^{-i \frac{3 \pi}{7}}}{2} e^{-i 800\pi t}$$

What is the spectral amplitude in each case? I am unclear if the frequency of $800 \pi t$ has amplitude $5$ or $\frac{5 e^{-i \frac{3 \pi}{7}}}{2}$.
If it is the second do I make two seperate plots for the complex and real parts, or do I just plot $\bigg\lvert \frac{5 e^{-i \frac{3 \pi}{7}}}{2}\bigg\rvert $?

Given your signal $x(t) = 5\cos(800\pi t + \frac{3 \pi}{7})$. And considering your correct decomposition of it into complex exponentials:

$$x(t) = \frac{5 e^{i \frac{3 \pi}{7}}}{2} e^{i 800\pi t} + \frac{5 e^{-i \frac{3 \pi}{7}}}{2} e^{-i 800\pi t}$$

You next apply the continuous-Fourier transform properties to get:

$$X(\Omega) =  5\pi \left( e^{i \frac{3 \pi}{7}} \delta(\Omega - 800\pi) + e^{-i \frac{3 \pi}{7}} \delta(\Omega + 800\pi) \right)$$

Now your Fourier transform has two impulses with weights of $5\pi e^{i \frac{3 \pi}{7}}$ and $5\pi e^{-i \frac{3 \pi}{7}}$ respectively.

The spectral amplitudes are infinite, by definition of the impulse function, but they have the weight respectively. And for the magnitude you would take the absolute value of the amplitudes, which is therefore infinite. Yet again you can represent the magnitudes of those impulses by their weights' absolute value.


Inter Symbol Interference (ISI) / Equalization in Digital Signal Processing I am studying the OFDM system and before to start I have decided to study the main concepts of modulation. And I have spent yet a lot of time to understand the concept of an equalization to reduce ISI.

I have found so many lectures and ppts, which provide and explain that. But it doesn't make sense for me.

What I understand, that the equalizer is a digital filter, that should be applied after receiver part. Is it right? And it reduces the ISI. A zero-forcing equalizer reduces the ISI complete.

But I can't find more detailed information about how the equalizer reduces the ISI. What is a principle of this type of the digital filter/ equalizer?

What is mean in this picture?

In OFDM DFT and IDFT work as an equalizer, sampling...

One more question: What is a role of a pulse shaper?

If the channel just adds noise and does not cause any distortion (as in the first figure), you can design the transmit filter $g(t)$ such that, when combined with its matched filter at the receiver, the Nyquist criterion for zero intersymbol interference (ISI) is satisfied. If the channel adds linear amplitude and phase distortion (as in the second figure), then generally the Nyquist criterion will not be satisfied and you will observe ISI, which will increase the error probability.

A linear equalizer is designed to (at least partly) compensate the effect of the channel $h(t)$. There are different optimality criteria, such as the zero-forcing criterion you mentioned, which eliminates ISI completely, but which does not take into account the noise. A zero-forcing equalizer basically inverts the channel response, if this is possible at all (with a causal and stable filter). If the noise is not negligible, it is usually better to use a different criterion such as the minimum mean square error (MMSE) criterion, which jointly minimizes the effects of ISI and noise.

Since many channel responses are unknown and/or time-varying, i.e., the impulse response $h(t)$ changes over time, equalizers must be adaptive. In practice, most adaptive equalizers use a sort of least mean squares (LMS) adaptation algorithm because it is relatively simple and robust.

This topic is quite broad and complex, so I would recommend that you read a good introductory text on digital communications, instead of gathering unrelated pieces of information from the internet.

Also take a look at this related question and its answer.


IIR filter and phase shift I heard that an IIR filter has delay and phase shift effect, but not clearly what are they. For example, 1-pole Lowpass has 90 degree phase shift, 2-pole Lowpass has 180 degree? So does it mean two 2-pole Lowpass filters has 360 degree phase shift? Beside the HighPass, LowPass, AllPass filter, other IIR filter like LowShelf, HighShelf, Peak... has phase shift too?
When I apply an IIR filter using biquad transfer function, does it has any phase shift diffirent if I use direct form I (4 delay registers) vs direct form II (2 delay registers)?
If IIR filter has phase shift, the phase shift occurred at the filter's target frequency, or all the frequencies? How can I phase shift (all frequencies) 90 degree clockwise & counterclockwise just using AllPass filter?
I know some of my questions are noob, but I'm just study dsp myself, and there're somethings that I really confused about. Sorry for my bad English too :)

I believe you have received no answers to your questions because your first three questions are so strange, so puzzling, that no one where to start to even begin to compose answers.  See:

http://www.bing.com/videos/search?q=digital%20filter%20phase%20response&qs=n&form=QBVR&pq=digital%20filter%20phase%20response&sc=0-0&sp=-1&sk=#view=detail&mid=3DC56A834B9C701B5FE93DC56A834B9C701B5FE9

and

http://www.bing.com/videos/search?q=digital%20filter%20phase%20response&qs=n&form=QBVR&pq=digital%20filter%20phase%20response&sc=0-0&sp=-1&sk=#view=detail&mid=85921946745DEDD05D7385921946745DEDD05D73

Good Luck.


Meaning of Rect and Train of Rect Spectra The Fourier transform of $x(t)=\operatorname{rect}(t)$ is $X(f)=\operatorname{sinc}(f)$

The Fourier transform of a periodic train of rectangular pulses $x(t)=\sum\limits_{n=-\infty}^{\infty}\operatorname{rect}(t-n)$ is a sampled sinc:

Picture taken from here.
From a math point of view it is fine because:

The periodic train of rects is a periodic signal that can be expressed as a discrete sum of sine waves according to Fourier Series Expansion.

The periodic train of rects is the convolution of a rect and a train of dirac deltas. The product of their Fourier Transforms is a sampled sinc since they are a sinc and again train of dirac deltas in frequency domain.

So, my question is a about the physical meaning of this. A single rect has a continuous spectrum, that means that it has a certain amout of power (and energy) at each possible frequency (from $-\infty$ to $+\infty$). Why, if I repeat the rect infinite times, do I remove some frequencies from its spectrum and leave only a discrete set of frequencies? Why does repeating periodically a signal means removing some frequency intervals from the single pulse spectrum?
In time domain, frequency may be seen as the speed of variation of the signal. Well, it is the same for a single rect and a train of rects (if we suppose the same fall and rising times and same rect duration).

Why does repeating periodically a signal means removing some frequency intervals from the single pulse spectrum?

Because your signal is periodic. A periodic signal only contains frequencies that are integer multiples of the signal's repetition frequency. In other words, if the signal is periodic all of it's constituent parts must be periodic as well.

In time domain, frequency may be seen as the speed of variation of the signal.

Not really. "speed of variation" would be the first derivate of the signal. It's somewhat related to higher frequencies since the Fourier Transform of a differentiator is proportional to the frequency itself but it's not the same thing. If you look at a square wave you also need high frequency to keep the flat parts nice and flat.

So, my question is a about the physical meaning of this.

That's always a tricky proposition. One possible answer is "there is no physical meaning since it's physically impossible to have a signal that's infinite in time". Sine waves are a useful mathematical concept but they DO NOT EXIST in reality. Of course, for practical purposes you can always generate something that's "close enough" for your specific application.
Another explanation
Let's look at what happens if you repeat the rect just once. In the time domain you can represent this as a convolution with two pulses, i.e.
$$h(t) = \delta(t) + \delta(t-T_0) $$
So you add a delayed copy to the original signal. In the frequency domain this looks like
$$H(\omega) = 1 + e^{-j\omega T_0}$$
The delay turns into a phase shift so you add add a phase shifted copy of the spectrum the original one. If the phase shift at a certain frequency is 0 degree, they will just add. But if the phase shift at a different frequency is 180 degrees they will actually cancel. So some frequencies will be amplified others will be cancelled that's called a "comb filter" a graph of the amplitude with these notches looks like the teeth of a comb.
The more you repeat the original rectangle, the more "combing" you get and eventually all frequencies are gone except the ones that line up exactly with the period of the repetition.


Determine the Given system is linear or nonlinear Does the following define a linear or nonlinear system?

$$ y(n) -4 y(n)y(2n)=x(n) $$

Using the definition of linear system as one which can be put in the form:

$y_0 + a_1 y_1 + \dots + a_n y_n = b_0 x_0 + \dots + b_m x_m$

because of the multiplicative term $y(n)y(2n)$ we can conclude that the system is nonlinear.


FM: Do rational ratios of carrier to modulator signals really produce harmonic spectra? I've read that in FM, all rational ratios of the carrier frequency to the modulator frequency produce harmonic spectrum in the FM signal. Link to what I read.

I don't understand why this is. If for example carrier frequency is 10Hz and modulator frequency is 5Hz, then the spectrum will include the frequencies 15Hz, 25Hz..., which are not integer multiples of 10Hz.

And for some ratios, the frequencies in the sequences come in "pairs",  the ratio 10:4 gives the sequence 2,4,8,10,14,16..., which doesn't look like a harmonic series at all.

What is the proof for this?

That paper (1973) is about more natural, simple and rich audio synthesis using FM (which later became Yamaha's very populer FM synthesis audio chips).

In that respect, the carrier frequency is taken very low compared to the modulating tone frequency, unlike in communication applications. Hence the the side lobes cross the zero frequency boundary and are reflected back into the positive band using simple trigonometry.

Another alternative to see this is to use positive and negative frequency version (exponential) of Fourier series analysis, to see that negative frequency side lobes penetrate into the positive (vice versa) band when the carrier frequency is insufficiently low.

Note that the resulting Bessel spectra of the single tone modulated FM signal predicts that the side lobes about the carrier will be placed at the frequencies of:

$$f_n = f_c \pm n f_m $$
where $f_c$ is th carrier frqeuency, and $f_m$ is the single tone modulator. As you can see, if $f_c$ is sufficiently low, then for some $n$, the component at $f_n = f_c - n f_m$ will be less than $0$ Hertz and will be reflected into the positive frequency $f_r = n f_m - f_c$.

So the FM side lobes that reflect from negative frequencies and the already existing ones do mix up to create the final audio spectrum which can be harmonic or inharmonic depending on the carrier frequency and modulating tone frequency ratios.

This was a genuine way of producing harmonically rich audio synthesis in legacy analog synthesizers.


Number of antenna to perform MUSIC algorithm When the MUSIC algorithm is used for direction of arrival (DOA) estimation it can identify up M-1 signals, where M is the number of array elements. MUSIC can also be used for spectral frequency estimation. Does MUSIC require more than one receiver (or more than one receive antenna) for frequency estimation of a single signal?

In the field of array signal processing, the signal model is
$$
\mathbf{x}(n) = [x_0(n), x_1(n), \ldots, x_{M-1}(n)]^T
$$
where $M$ denotes the number of sensor and $n$ is the discrete time.
While in the field of spectral estimation, the signal model is
$$
\mathbf{x}(n) = [x(n), x(n-1), \ldots, x(n-M+1)]^T
$$
where $M$ is the number of samples of the observed signal, and $x(n-m)$ represents the previous signal before $m$ samples.
The autocorrelation matrix is
$$
\mathbf{R}_{xx} = E\{ \mathbf{x}(n) \mathbf{x}^H(n) \}
$$
and it is usually replaced by its estimation $\mathbf{\hat{R}}_{xx}$ in practice.


Inverse Fourier Transform Dirac impulse with scaled argument Currently, I am dealing with the sampling problems and I don't understand how to calculate inverse Fourier transform of a scaling impulse function

$\textrm{IFT}\{\delta(\Omega T)\} = ?$, $T$ is sampling period but we can see it as a constant.

It comes from this equation since I have to find the $y(t)$ of reconstructing signal

$$Y_r(j\Omega) = j\omega_0\pi\delta(\Omega T - \omega_0) - j\omega_0\pi\delta (\Omega T + \omega_0) $$

The solution for $\textrm{IFT}\{\delta(\Omega T)\}$ is $\frac{1}{2\pi T}$, but I don't know how to get them.

I got stuck when trying to expand and calculate it by integral of inverse FT.

This is just a matter of variable substitution in the integral. The inverse Fourier transform of $\delta(\Omega T)$ is

$$\mathcal{F}^{-1}\big\{\delta(\Omega T)\big\}=\frac{1}{2\pi}\int_{-\infty}^{\infty}\delta(\Omega T)\,e^{j\Omega t}\,d\Omega\tag{1}$$

Substituting $\alpha=\Omega T$ gives $d\alpha=d\Omega\cdot T$, and $(1)$ becomes

$$\mathcal{F}^{-1}\big\{\delta(\Omega T)\big\}=\frac{1}{2\pi T}\int_{-\infty}^{\infty}\delta(\alpha)\,e^{j\alpha t/T}\,d\alpha=\frac{1}{2\pi T}\tag{2}$$

The last equation is a consequence of the sifting property of the Dirac delta impulse.


ifft of signal without phase information I have a chemical spectrum with amplitude on the Y axis and wavenumbers on the X axis, sampling interval I think this is $1.5^{10}$ hertz. I would like to reconstruct the amplitude/time representation of the original signal by performing a ifft of the signal however I have no phase information. Since the information is unchanging (constant) over time I don't think that the phase information is necessary for reconstruction of the original signal. Can I set the phase information to zero and still reconstruct the signal in this case and how do I incorporate the sampling interval in the matlab ifft routine for accurate frequency calculation?

This is what I think I should be using

sig = abs(ifftshift(ifft(complex(data.YData))));

Thanks

Dave

This is tricky. The phase information impacts the time domain form very much. For example white noise and a delta impulse have exactly the same amplitude spectrum, the only difference is in the phase.

Without any further assumptions or extra knowledge this can't be done. For example if you know that the signal is stationary, you could just try a random phase. If it's the impulse response of a causal system, you could try a minimum phase.


Understanding the DCT I have been trying to understand the definition of the DCT and found a really good website describing its definition.

Here is the link DCT.

What I don't understand is the part of the DFT. Why is the summation of the DFT multiplied by $\frac{1}{\sqrt{2N}}$?

Thank you.

Given $\{x_n, 0 \leq n \leq N - 1\}$, the conventional DFT is defined as
$$X_k = \sum_n x_n e^{-j2\pi kn/N}$$
The Parseval theorem gives
$$\sum_{n=0}^{N-1} |x_n|^2 = \frac{1}{N}\sum_{k=0}^{N-1} |X_k|^2$$
The DFT in your link is called Normalized DFT (NDFT):
$$Y_k = \frac{1}{\sqrt{N}} X_k = \frac{1}{\sqrt{N}} \sum_n x_n e^{-j2\pi kn/N}$$
that
$$\sum_{n=0}^{N-1} |x_n|^2 = \sum_{k=0}^{N-1} |Y_k|^2$$
Let me cite the idea behind this Normalized DFT (please follow the link if you desire more details)

While this definition is much cleaner from a geometric signal theory point of view, it is rarely used in practice since it requires slightly more computation than the typical definition. However, note that the only difference between the forward and inverse transforms in this case is the sign of the exponent in the kernel. Advantages of the NDFT over the DFT in fixed-point implementations (Appendix G) are discussed in Appendix A.
It can be said that only the NDFT provides a proper change of coordinates from the time-domain (shifted impulse basis signals) to the frequency-domain (DFT sinusoid basis signals). That is, only the NDFT is a pure rotation in $\mathbb{C}^N$ , preserving both orthogonality and the unit-norm property of the basis functions. The DFT, in contrast, preserves orthogonality, but the norms of the basis functions grow to $\sqrt{N}$ . Therefore, in the present context, the DFT coefficients can be considered denormalized frequency-domain coordinates.


IDTFT of convolution in the frequency domain I have tried everything. If you actually know how to solve this could you provide a hint?
$$ e^{-2j\Omega}\frac{ \sin\left( \frac{7\Omega}{2}\right)}{ \sin\left( \frac{\Omega}{2} \right)}\star \frac{\sin\left( \frac{10 \Omega}{2} \right)}{\sin\left( \frac{\Omega}{2} \right) }$$
Ideally I would like to find the Fourier of every "fraction" in separate and then use properties:
$x(n - n_o) \rightarrow e^{-jn_0\Omega}X(\omega)$ so I don't mind for $$ e^{-2j\Omega}$$
but I have 2 problems:

I cannot use $\displaystyle \frac{\sin\left(\left(n+\frac 12\right)\Omega\right)}{\sin\left(\frac \Omega 2\right)}$ for $(n+1/2) = 10/2$ because $n \in Z$
In DTFT , in my book there is no property like in continous time to transform convolution in $\Omega$ domain to multiplication in time domain so I don't know what to here as well.

update:
After some comments and help from people who answered :
I am going to try to do it as juch $\frac{sin(10\Omega/2)}{sin(\Omega/2)}= \frac{sin(10\Omega/2)}{sin(\Omega/2)}e^{-j\Omega(10-1)/2}e^{j\Omega(10-1)/2}=\Big[\frac{sin(10\Omega/2)}{sin(\Omega/2)}e^{-j\Omega(10-1)/2}\Big]e^{j9\Omega/2}$
I am to take advantage of the property : $\Big[\frac{sin(10\Omega/2)}{sin(\Omega/2)}e^{-j\Omega(10-1)/2}\Big]e^{j9\Omega/2} \rightarrow 2\pi F^{-1}{\Big[\frac{sin(10\Omega/2)}{sin(\Omega/2)}e^{-j\Omega(10-1)/2}\Big]} * F^{-1}[e^{j9\Omega/2}]$
The result is :
$F^{-1}[e^{j9\Omega/2}] =$
$\frac{1}{2\pi}int_{\pi}^{\pi}e^{j9\Omega/2}e^{j\Omega n}d\Omega = \frac{1}{2\pi}\frac{e^{j\Omega(9/2 +n)}}{j(9/2+n)}\Big|_{-\pi}^{\pi}=\frac{4(-1)^n}{2\pi(n+9)}$ ( i think)
and $F^{-1}[e^{j9\Omega/2}]=1$ for $n \in [0,9]$ and 0 anywhere else.
Now we need to compute the convolution of those 2:
the result should be non - zero only when $n \in [0,9]$ so:
$F^{-1}\Big[\Big[\frac{sin(10\Omega/2)}{sin(\Omega/2)}\Big]e^{-j\Omega(10-1)/2}\Big] = \begin{cases} \frac{4(-1)^n}{(n+9)} & n \in [0,9] \\
                         0 & else \end{cases}$

It seems to me that this exercise is meant to combines  the basic properties of the DTFT (gathered here: Table of DTFT Properties). Matt described the product/convolution property. You also get the time-shift/complex modulation. I suspect (thought I did not do the calculations) that the issue with the factor $10$ could be tackled with a change in variable: $10\Omega = 5\times (2\Omega)$, and the use of the time-scaling (expansion) property (see Time scaling of discrete-time sequences and the DTFT):
$$ S(c\Omega) \Longleftrightarrow \left\{
                \begin{array}{ll}
                  x[n/c] \textrm{ if } n/c \textrm{ is an integer } \\
                  0 \textrm{ otherwise.}
                \end{array}
              \right.$$
Other hints hinge around Dirichlet kernels:
$$D_N(x) =\frac{\sin\left(\left(N +1/2\right) x \right)}{\sin(x/2)}$$
They are also called asinc or psinc (aliased or periodic cardinal sine or sinc), and related to finite-support discrete windows. If $*$ is the Convolution sign, the resolution may use the product/convolution properties of Fourier.
And more generally (at the bottom of the page on Dirichlet kernels), you have the identity:
$$\sum_{n=0}^{N-1} e^{jn\Omega} = e^{j(N-1)\Omega/2}\frac{\sin(N \, \Omega/2)}{\sin(\Omega/2)}\,,$$
which relates the DTFT of a discrete time window $w_{[0,N-1]}$ (from index $n=0$ to $n=N-1$) to ratios of sines with a phase correcting term. You can check details at Discrete-time Fourier transform of a window function.


Fixed point conversions Have following questions about Fixed point conversion format.

Suppose I want to change Fixed Point format as

X1(m1,n1) (→) X2(m2,n2),m1+n1=m2+n2[Word length same],  n2>n1[ Higher Fractionlength]

Suppose  X1=00001100 and suppose I do X1<<(2^3)

Q1)  If I interpret X1 as pure Signed integer, then X1= 12 and above operation means
 Multiply by 2^(n2-n1)=8 so that X2= 01100000 = 96. (Decimal point is irrelevant here?). Is the interpretation correct?
  or
 Assuming Fixed point  notation, above means
 X1(m1,n1) (→) X2(m1-K,n1+k),m1+n1=m2+n2[Word length same],  n2>n1[ Higher Fractionlength]  and k=n2-n1
  For given example, this means X1(3,5)=000.011 00(0.375)(→) X2(0,8)=  .01100000=0.375
  So, the value remained same.

Q2) Can we always interpret this as decimal point being fixed and bits moved left?Is this a logical shift or arithmatic shift?

Q3) Ideally, we would like to preserve resolution and range covered. Under what condition would this be guaranteed for the above shift? for example, consider following conversions:

Consider  X1(3,5)=101.011 00(5.375)(→) X2(0,8) is 101.01100<<3 = .01100000 =0.375.So, we have error.
   Q4) Is this rounding or truncation and why?

Consider  X1(m1n1) = 000...0 . abc......z  , [m1 = all zeroes],  X1 << k = X1(m1-k,m1+k) where k<=m1
            would  not change the value, as bits to the left and right  of decimal point remain same

Q5) So, if fixed point fractional notation is assumed( Positive fraction, with magnitude <1), Left shifts smaller than integer word length would preserve the fractional value?

Suppose I have fixed point number X1(m1,n1) with Word length m1+n1. Suppose  I want to express
it in X2(m2,n2) format so that Word length is reduced ,Fraction length is increased, Integer Word length reduced    but the value remains the same.For example:
X1(3,5) = 000.01100(0.375) -> X2(1,6)
Steps for above could be:

Left shift by n2-n1 bits so that 000 011 00 <<1 = 000 110 00

Keep Least Significant m2+n2 bits so that X2 = 0.011000(0.375) or Keep Most Significant m2+n2 bits so that X2 = 0.001100(0.1875)

Q6) If I was writing a routine that would automatically do conversion, )(say this was Sum result that needed to be converted to output (m2.n2) format automatically), how are what decides whether I keep LS or MS bits or doing a left or right shift?
Suppose I have fixed point number X1(m1,n1) with Word length m1+n1. Suppose  I want to express  it in X2(m2,n2) format so that Word length is increased ,Fraction length is increased but the value remains the same.For example:
 X1(3,5) = 100.01100(-3.625) -> X2(3,13)
and  X1(3.5) = 100.01100(-3.635) -> X2(1,11)

Q7) How can I go about that?  What would be steps if result required increase in Word length and Integer Word Length?
Suppose I need to convert Fixed Point number as per
X1(m+1,n) (→) X2(m,n)[ Word length reduced, Fraction length same, integer Word length reduced]

One of the references I found says do following :
  (x1)&(0xff..e) and keep MS m2+n2 bits
  With X1(1,4) = 0.1011(0.875) this becomes 0.101 = 0.75

Other says: Do left shift by 1 and Keep Ms m2+n2 bits
With X1(1,4) = 1.011(0.875) this becomes 1.011(0.625)

Obviously, i don't perverse the value.
Q8) Is there a correct way to do it and the logic behind it?

I know these are lot of questions but it will help me clear lot of confusing stuff.
Thanks
sedy

As mentioned in the comments, your question is pretty confusing and confused, so I might not be able to answer all your smaller questions, but I'll try to shed some light on the topic. This should help you understand what's going on, and hopefully this understanding will help you answer all further questions by yourself.

In the following I assume two's complement representation. Note that if $M$ is the word length, and $n$ is the number of fractional bits, a number $x$ is represented as

$$x=K\cdot 2^{-n}\tag{1}$$

where $K$ is an integer is the range $[-2^{M-1},2^{M-1}-1]$. Now if you want to change the number of fractional bits to $m$, you have

$$x=L\cdot 2^{-m}\tag{2}$$

Combining (1) and (2) gives

$$L=K\cdot 2^{m-n}\tag{3}$$

If $m>n$, $K$ needs to be left-shifted by $m-n$ bits to obtain $L$. This means that you lose range (i.e. the maximum number that can be represented becomes smaller), so there is a possibility of overflow. If $m<n$, $L$ is obtained from $K$ by right-shifting $n-m$ bits, but for two's complement you need to replicate the sign bit, i.e. if the original sign bit was set, the new right-shifted number must also have the sign bit set. For $m<n$ you can lose precision because the resolution is reduced.

An example makes things easier to see. Let $M=6$, i.e. the integer $K$ is in the range $[-32,31]$. Let the number of fractional bits be $n=3$. E.g.

$$x=\texttt{101.011}=(-2^{5}+2^{3}+2^{1}+2^{0})\cdot 2^{-3}=-21\cdot 2^{-3}=-2.625$$

If we want to change the number of fractional bits to $m=2$, you need to do an arithmetic right-shift, resulting in

$$y=\texttt{1101.01}=-11\cdot 2^{-2}=-2.75\neq x$$

i.e. precision is lost. If instead we want $m=4$ fractional bits we need to do a left-shift resulting in overflow:

$$y=\texttt{01.0110}=22\cdot 2^{-4}=1.375\neq x$$

This is clear because by increasing resolution we have decreased the range of possible numbers to $[-32,31]\cdot 2^{-4}=[-2,1.9375]$, which makes it impossible to represent the original value $x=-2.625$.

There are of course examples where you would neither lose precision nor would you get overflow. E.g. (again with $M=6$ and $n=3$)

$$x=\texttt{001.010}=10\cdot 2^{-3}=1.25$$

For $m=2$ we get

$$y=\texttt{0001.01}=5\cdot 2^{-2}=1.25 = x$$

i.e. no loss in precision, and for $m=4$ we get

$$y=\texttt{01.0100}=20\cdot 2^{-4}=1.25=x$$

i.e. no overflow.

I hope that this addressed at least some of your questions and that things have become a bit clearer. If you haven't done it yet, you should check out the basics here, here, and here.


FFT-based coarse carrier recovery for QPSK An FFT-based coarse frequency acquisition block is being tested in GNU Radio. The algorithm used can be found on the GOES satellite user manual and is summarized in the figure below.

The algorithm works by:

Carrier detection by absolute-law, square-law or fourth-power law. Square law detector involves squaring the IQ signal while fourth-power law uses the fourth-power of the IQ signal.
Finding the FFT
Averaging successive FFT sequences
Finding the peak by looking at 8 adjacent bins with the highest magnitude.

I have been able to implement the algorithm in GNU Radio, at least for BPSK. The algorithm is able to recover the carrier at all offset scenarios (offset < 10% symbol rate, offset ~ symbol rate, offset > symbol rate).

According to the GOES manual, the same algorithm should work for OQPSK, and I assume by extension, QPSK. However, the carrier derived by the algorithm for the case of QPSK is pretty bad, especially in low-SNR scenarios.

I have two questions

Are there any literature sources as to how the frequency detectors (square, fourth-power) work? I have basically implemented the algorithms but I do not fully understand their mathematical derivations.
Is there a more effective frequency detector for QPSK apart from square/fourth-power/absolute laws that could remedy the situation?

Regards,
Moses.

EDIT

Disclaimer: Using a fouth-law as proposed by @DanBoschen resulted in expected BER performance. For all intents and purposes, the answer is already accepted. But, I have observed some behavior of the nature of FFT relating to this question. I figured there wouldn't be any need to open a new question since we are already here.

In the flowgraph, the sample rate used is 8Mhz. The fourth-law-derived carrier can be tracked from -4Mhz to +4Mhz. There seems to be an "FFT ambiguity" when the offset frequency is above $\frac{1}{8}$-th the sample rate i.e anything above 1Mhz. Taking, for example, an offset of 1.5Mhz. The fourth-law detector will create a carrier at 6Mhz, which will be wrapped to -2Mhz in the second Nyquist zone. Dividing this by 4 results into an offset of -0.5Mhz, instead of 1.5Mhz, as shown in the figure below (derived carrier in green).

We can solve this problem by simply using a sample rate which is 8x the maximum offset for the fouth-law or 4x the maximum offset for the square-law. But, of course, that's not an ideal solution. I was wondering if there are better ways.

A square operation creates an unmodulated tone for a BPSK signal at 2x the carrier frequency (a pure tone for the case that the signal was unfiltered or rectangular pulses with perfect phase and amplitude balance in the BPSK modulation, and typically a stronger carrier with weaker sidebands in the more common filtered or pulse-shaped cases). For QPSK signals (including OQPSK) a fourth law operation is required which produces the same at 4x the carrier frequency.

This can be explained by knowing that multiplying signals causes their phases to add. So a signal multiplied by itself (square) that is modulated 0° to 180° would then become 0° to 360° which is the same as 0° (unmodulated). For QPSK the states are 0°, 90°, 180° and 270° where it can be seen that multiplying any of those phases by 4 results in 0° using the same modulo 360 operation when adding phase.

It is typical when using this approach for carrier recovery to also use a PLL to track the recovered 4x carrier in order to clean up the residual (much lower) sidebands resulting from the finite bandwidth modulated signal and then the PLL-filtered signal is frequency divided to complete the carrier recovery.

I tend not to use this approach given the higher sampling frequency required and the comparative complexity to other methods when implementing all digital systems but this would make for a relatively simple analog approach using analog frequency multipliers. I detail another approach for carrier recovery specific to BPSK, QPSK and QAM that would be in my opinion much more efficient for a digital radio implementation in this post: High modulation index PSK - carrier recovery

This is also an interesting approach by fred harris for a frequency lock loop that I believe would work for QPSK carrier acquisition with additional subsequent phase tracking: How does this FLL work?


Steady state value of a complex convolution I am trying few problems on the introductory part of DSP. One of the problem asks to calculate the steady state response of a system with impulse response $h[n] = (\frac{j}{2})^{n} u[n] $ to an input signal $x[n]=\cos(\pi n) u[n]$. However, I have no idea how to handle this complex $j$ in the convolution. I want someone to confirm if my solution is correct or not. If it is incorrect please tell me where I am getting wrong.

I start with calculating the Z transform of $h[n]$.
$$\begin{align}
H(z) &= \sum_{n=0}^{\infty}h[n] \, z^{-n} \\
 &= \sum_{n=0}^{\infty}\left(\frac{j}{2}\right)^{n} \, z^{-n} \\
 &= \frac{1}{1-\frac{j}{2}z^{-1}} \\
\end{align}$$
such that $|z| > \frac{1}{2}$.

Now the z transform of $x[n]$ is ,
$$\begin{align}
X(z) &= \sum_{n=0}^{\infty}\cos(\pi n)z^{-n} \\
 &= \sum_{n=0}^{\infty}\Big[\frac{e^{j\pi n}+e^{-j\pi n}}{2}\Big]z^{-n} \\
 &= \frac{1}{2} \left(\frac{1}{1-e^{j\pi}z^{-1}}+\frac{1}{1-e^{-j\pi}z^{-1}}\right) \\
 &=\frac{1}{1+z^{-1}} \\
\end{align}$$

The convolution output in the z domain is $Y(z)$.
$$Y(z) = H(z)X(z) = \frac{2}{(2-jz^{-1})(1+z^{-1})}$$
Applying the final value theorem, I get, steady state output $y_{ss}$ as
$$y_{ss} = \lim_{z \to 0} \, zY(z) = 0$$.

I definitely think that I am doing something wrong over here. Please help. I don't know how to handle this $j$ in the convolution operation like this, please help me in that too.

You can't apply the final value theorem here. First of all, the correct form of the final value theorem for the $\mathcal{Z}$-transform is

$$\lim_{n\to\infty}y[n]=\lim_{z\to 1}(z-1)Y(z)\tag{1}$$

Eq. $(1)$ assumes that the poles of $(z-1)Y(z)$ are inside the unit circle, which is not the case in your example. In your case, the limit $(1)$ does not exist.

What you have to do is find the frequency response of the system, which is just the transfer function evaluated on the unit circle $z=e^{j\omega}$ (if it converges on the unit circle). Note that the region of convergence (ROC) you've found does not make sense because $|z|$ must be a real number. The requirement for the sum to converge is

$$\left|\frac{j}{2z}\right|<1\quad\Longrightarrow\quad |z|>\frac12\tag{2}$$

Since $H(z)$ converges on the unit circle $|z|=1$, the frequency response is given by

$$H(e^{j\omega})=\frac{1}{1-\frac{j}{2}e^{-j\omega}}\tag{3}$$

An LTI system's response to a sinusoidal input $x[n]=\cos(\omega_0n)$ is

$$y[n]=|H(e^{j\omega_0})|\,\cos[\omega_0n+\phi(\omega_0)]\tag{4}$$

where $\phi(\omega)$ is the system's phase response:

$$H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{5}$$

So in the limit, after all transients have died out, the system's (steady-state) response to a switched sinusoidal input $x[n]=\cos(\omega_0n)u[n]$ is given by $(4)$. What is left now is to determine the magnitude and phase of $H(e^{j\omega})$ given by $(3)$ and plug these values into Eq. $(4)$.


What is the state-of-the-art of using signal-processing for "cleaning-up" videos? TWO PART QUESTION - We've all seen the movies/TV shows where the police/feds/spies use computer software to take a grainy photo and do a "clean-up" to see a better picture and more details. I assume the concept is executed by some sort of uber-power pixel-smoothing or anti-aliasing type of algorithms to fill in the blanks based on deductive processing.

PART 1: How real is this technology in the public/commercial software world? I'm not asking about any speculation on alleged secret gov software or such, I just want to know where we actually are with this concept today? How much is fully automated vs human-assisted.

PART 2: Assuming there actually is some reality with this technology for photographs the second part of this question is how (if at all) has this been applied to videos? Again the issue of fully automated vs human assisted is of interest here.

At the heart of this post is the ultimate question of how viable is today's software for being able to take an old VHS or DVD recording and process the frames to create a new HD-resolution remaster. Considering that doing this would mean cleaning up tens-of-thousands of frames for even a simple wedding video I am not expecting this technology to be fast of course.

NOTE: Per in-topic and meta discussions I went ahead and cross-posted this in two other suggested SE forums to acquire their viewpoints and expertise on this matter. So far (it is still a little "early in the day" so to speak) I have received some pretty interesting information in 2 of the 3.

When this is all done (when good answers have been selected) I would appreciate a way to merge these for the benefit of all three SE communities:

Computer Graphics: https://computergraphics.stackexchange.com/questions/1880
Signal Processing: What is the state-of-the-art of using signal-processing for "cleaning-up" videos?
Video Production: https://video.stackexchange.com/questions/17363

Entertainment industry greatly exaggerates what video enhancement can do.

Spatial and temporal noise reduction

Noise reduction by Neat Video and MSU Denoiser Filter, which are somewhere at (according to a comparison from 2007) or near state-of-the-art for hard noise, takes advantage of that adjacent frames usually are of the same scene but with independent noise. So noise can be filtered both spatially and temporally. But that kind of processing is in my experience just cosmetic. It does not really reveal things that a human would not already see through the noise, it just makes the video look nicer. Both filters have an option to adjust settings automatically, but because noise reduction seems always a compromise between noise reduction and introduction of artifacts, the final "strength" of noise reduction should be adjusted to taste manually.

Motion estimation and superresolution

MSU Super-Resolution Filter reconstructs a new video at 2x horizontal and vertical resolution. It follows a similar paradigm as the above two filters, but also estimates motion between frames. The original pixelization (and other artifacts that do not follow the motion such as blocking and sensor Moiré) can be reduced by compensating for motion and then filtering temporally. This improves resolution and makes detailed things like pixel-deprived text and faces more clear. Resolution can be increased only if the pixelization is not properly anti-aliased. Artificial pixelization of faces, if not done carefully, could be reduced by the same or similar processing. This is a kind of enhancement that the human visual system is not accustomed to do natively, and can actually reveal hidden useful information.


Output of system defined by differential equation I don't fully understand how the output of a system can be derived from the system's differential equation and a given input.

For example:

$$y(0-) = 1 $$
$$y'(0-) = -2$$
$$ u(t) :   \text{Heaviside function} $$
$$ S:  y''(t) +2y'(t) +y(t) = u''(t) -2u'(t) + u(t)$$

Can someone show me a good way to find $y(t)$, given $u(t)$ and $S$?

What you want to do is take the Laplace transform of both sides of S.  The first and second derivative terms will be replaced by expressions of $Y(s)$ and the conditions at $t=0-$, so you can solve the equation for $Y(s)$ and then do the inverse Laplace transform, or more commonly, use partial fraction expansion to go back to the time domain.

See numbers 35 and 36 in the table of Laplace transforms here.


Is there use of $\pi/8$-DQPSK as shown in MATLAB's documentation? I was reading MATLAB's documentation on dpskmod. They gave an example of an $\pi/8$-DQPSK. I felt that this is strange as I have seen $\pi/4$-DQSPK and $\pi/8$-D8PSK, but I have not seen $\pi/8$-DQPSK.

Attached is a screenshot of the documentation showing $\pi/8$-DQPSK with 16 clusters in the constellation:

My question is: Is $\pi/8$-DQPSK actually used in industry or is this just a demonstration of the flexibility of MATLAB's functions?

I believe the correct name is pi/8 D8PSK. The ETSI TETRA Release 2 uses pi/8 D8PSK (digital trunked mobile standard for Professional Mobile Radio and Public Access Mobile Radio). The purpose of pi/4 DQPSK and pi/8 D8PSK is to have lower peak to average power ratios by not allowing the trajectories to go through the origin. Lower peak to average ratios means you can drive power amplifiers in the transmitter further into saturation, resulting in higher power efficiency (more transmit power from available DC power).   pi/8 provides for higher data rates at the expense of the more SNR that would be needed to get the same error rate performance.

More info on TETRA:

(see page 40 for details on pi/8 modulation which is identical to what the OP has shown):
https://www.etsi.org/deliver/etsi_tr/102500_102599/102580/01.01.01_60/tr_102580v010101p.pdf

Motorola TEDS:
https://www.motorolasolutions.com/en_xu/products/tetra/teds.html


A Practical Method for Correcting Phase Errors in a Transceiver System Comprising of OFDM + DUC and DDC I am back with a question related to my OFDM Transceiver System that I am designing where I use DUC and DDC. I could get rid of the high BER problem. I realized that the problem was due to the phase error occuring between the transmitter and the receiver side as I saw the constellation diagram at the receiver rotating constantly. So to get rid of this, I obtained the phase difference between the transmitter and the receiver and I rotated the complex signal with the obtained phase error. However, I feel this is only an ideal method to get rid of the phase error. I read about the phase error compensation techniques but those are mostly when there is a noise in the channel. But in my case, I still haven't added any noise in the channel. I suspect the phase error is either because of the DUC and the DDC, or because of the pulse shaping filters that I used after OFDM modulator.

So I was wondering if anyone could tell me why is this phase error occuring even without a noisy channel, and what can I do to mitigate the phase error with some estimation.

I have attached the screenshot of my system along with the waveforms obtained. Below are the specs -

Random Integer Generator Specs

Set Size - 2

Initial Seed - 37

Sample Time - 1/1e6

Samples per frame - 20

OFDM Modulator Specs

FFT length - 32

Number of guard bands - [6;5]

Insert DC Null - Checked

Cyclic Prefix Length - 16

Number of OFDM symbols - 1

Number of transmitted antennas - 1

Raised Cosine Transmit Filter Specs

Rolloff factor - 0.2

Filter span in symbols - 8

Output samplers per symbol - 4

Linear amplitude filter gain - 1

DIgital Up- Converter Specs

Interpolation factor - [2 8]

Minimum order filter design - Yes

Two-sided bandwitdth - 2.4 e6

Passband ripple - 0.05 dB

Stopband ripple - 90 dB

Center Frequency  - 30e6

Input sample rate - 9.6e6

Digital Down-Converter Specs

Decimation factor - [8 2]

Minimum order filter design - Yes

Two-sided bandwitdth - 2.4 e6

Passband ripple - 0.05 dB

Stopband ripple - 90 dB

Center Frequency  - 30e6

Input sample rate - 153.6e6

Thank you in advance!

I don't think it would be reasonable to expect the phase to be exactly the same. The raised cosine filter for example has a linear phase response. Also don't forget that the auto-correlation for the cyclic prefix detection might find the frame starting at any of the 16 samples of the cyclic prefix, inducing an undetermined phase shift.

Do you have pilot sub-carriers? One simple way is to analyse their phase (which you know a-priori were transmitted with 0 (or some other known) phase) and rotate the sub-carriers in between the pilots (linearly interpolating between the phase of the pilots).


Relation between poles and cut off frequency of a filter I noticed that for a low pass Butterworth filter, as the cut off frequency reduces, the poles get closer to the unit circle. Is it applicable to all filters and why is this so?

As the normalized cut off frequency reduces from $\pi/2$ to $0$ (close to $0$- really low cut off frequency), the poles get closer to unit circle (the distance between the poles and unit circle reduces monotonically). Similar behavior is observed when we increase the cut off frequency from $\pi/2$ to $\pi$ (close to $\pi$- really high cut off frequency), the poles get closer to unit circle again- at $\pi/2$ they seem to be farthest away from the unit circle.

I would really appreciate if anybody could help me understand why this is so intuitively and also the math behind this (if there is any). Thanks in advance.

The Butterworth filter is in fact a continuous-time filter, and what you are seeing is very probably a bilinear transformation applied to this filter in order to get a discrete-time filter.

Software such as Octave or Matlab do this when you ask for a digital Butterworth filter.

The original Butterworth filter, in continuous time, $H(s)$ has all poles equally distributed around a semicircle in the left of the $s$ plane. The radius of this circle is the cutoff frequency (in $rad/s$).

The relationship between frequencies in discrete-time ($\theta$) and continuous-time ($\omega$) given by the bilinear transformation is: $\omega = \frac{2}{T} \tan (\theta/2)$ (where $T$ is an arbitrary sampling period used during design).

You can imagine the left $s$ semiplane as mapped into the unit circle in $z$, or the imaginary axis mapped to the unit circumference.


Periodicity of transfer function of FIR filter proof (Parks and Burrus, Digital Filter Design) In Digital Filter Design by Parks and Burrus, p. 19.

The transfer function of an FIR filter is given by the $\mathcal Z$-transform of $h(n)$ as:

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

(where $h$ is the filter)

The frequency response of a filter is defined as

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

where $\omega$ is frequency in $\textrm{rad/sec}$.

Then the text proceeds to show that $H(\omega)$ is periodic with period $2\pi$:

\begin{align}
H(\omega + 2 \pi)&= \sum_{n=0}^{N-1} h(n) e^{-j(\omega+2\pi)n}\\
&= \sum_{n=0}^{N-1} h(n) e^{- j \omega n} \color{red}{e^{-j2\pi n}}\\
&=H(\omega)
\end{align}
Could someone clarify how is that equal to $H(\omega)$ when there's the extra term $\color{red}{e^{-j2 \pi n}}$?


Windowing a signal for FFT I am actually working on a project involving the cross-correlation of two signals. My cross-correlation code is working fine, but I suspect that I should apply a window to my data set in order to avoid spectral leakage.

As I am not AT ALL specialist of signal processing, I have a few questions :

My goal is to precisely track the maximum value of the cross correlation. Which type of window should I chose ?
Applying the window before zero-padding the data set is the good way to do it, right ?
Not really related to signal processing as a theory, but are the Accord Framework windowing functions efficient ?

Thanks a lot !

I suspect that you suspect that you need to window because there is a huge literature on windowing, but all things FFT, aren’t about window.

At a minimum, we use a boxcar window so it makes sense to consider other windows but if the noise in your application is flat over the band the signal occupies, the boxcar is very likely the appropriate window.

If the noise isn’t flat, a whitening filter is  typically a better idea.

Spectral leakage is a problem when you have strong interference tones and the leakage masks a weaker set of tones that are your signal of interest.  A signal of interest leaking on itself is not usually a problem.  If you do have strong interference, an adaptive algorithm will provide an effective optimal window.

If accuracy is your main concern you really should be asking yourself what the bounds on accuracy will be and that will be governed more by FFT size, how many averages and what kind of averaging you do, which also leads to issues like how much your received signal changes.

The subject can be complicated.  There are literally thousands of papers in the IEEE library that touch on the subject.  You can dedicate yourself to thousands of hours of study, or you can just try a few windows on some data and pick one.

In many situations, a customer will not be interested in details, but some customers will and you can find yourself being asked questions about issues that you didn’t know existed.

I suggest you should know what a Cramer-Rao bound is and what it is for your problem, and Pearson Correlation.

So some answers.

If you minimally averaging, a window smooths the cross correlation. If you are minimally averaging, don’t throw away data.  There are too many factors to pick a specific type. You need to experiment.

If you window, do so prior to zero padding.


How to find correlation/cross-correlation of two signals in real time? Considering there are two signals and the signals are real time in nature. Visibly when one signal does up the other signal does down and vice versa. Sometimes both the signals moved towards the same direction and sometimes moves away from each other.
How to determine the juncture when both the signals change their direction.

Correlation is a multiply and accumulate process. So to do what you describe, simply multiply the two signals and accumulate (integrate) the resulting output. The integration time is set based on your knowledge of how long the waveforms could be correlated; a longer integration time will give you more processing gain against background noise but will also effect how rapidly you can observe changes in the correlation. If the integration time is longer than the interval over which the correlation changes the resulting correlation will be averaged out (therefore not observable).

Here is a specific approach that would minimize processing required while giving you a new result after every incoming sample, when we can assume synchronous inputs:

Multiply the two inputs, and place the result in a "sliding FIFO buffer" (FIFO = First In First Out). The buffer is of length N where N is the number of samples in our "integration time". This is only for purpose of maintaining what the resulting multiplication was N samples ago, as there is no need to do operations on every element in the buffer on each update or to keep a buffer of the input samples.

For each input, multiply the two and add the resulting product p to the buffer, which causes the last sample in the buffer to "fall off".  Add the result to an accumulator register r while subtracting the product from N samples ago (now in the last position in the buffer):

$r(n) = p(n) - p(n-N) + r(n-1)$

Where:

$r(n)$ is the current value to store in the accumulator

$p(n)$ is the product of incoming samples $x(n)$ and $y(n)$

$r(n-1)$ is the previous value that was in the accumulator before this update

$p(n-N)$ is the value of the product from N samples ago (get this from the FIFO buffer).

Thus at any given time r(n) represents the (unnormalized) correlation of x and y over N samples.

Just to note, and not important if you don't understand this part: This is a CIC (Cascade Integrator Comb) implementation of a moving average filter. A block diagram showing this process is shown below:

The above will give us a result that is not normalized. If a true normalized cross correlation was needed, we would need to scale the result by the number of samples as well as the standard deviation of x and y. Note that same approach can be used to determine the variance for both x and y: Use the same block diagram with x as both inputs to determine the variance of x (since the variance is the sum of the squares divided by N), and similarly y as both inputs to determine the variance of y. In most cases we would want to avoid this if we can be content with a result that is proportional to the correlation given the increased processing requiried due to the need to divide and do square roots (and three of the above block diagrams). If you have to go down this path, look into the Cordic Rotator for computing the square root if you can spare doing multiple operations between each incoming sample, and use an N that is a power of 2 so that the divide by N is a simple bit shift.


Corresponding frequency domain filter of spatial domain filter in image processing? My question is:
Suppose that we have low pass spatial domain filter that averages 4-connected neighbors of that pixel and it doesn't consider its pixel in averaging. Find its corresponding filter in frequency domain and show that it is a low pass filter.

I know that I must take forward Fourier transform from spatial domain filter and its output is frequency domain filter, but I don't have any idea what is spatial filter is? I think that it is $w(x,y)$ such that $-1 \le x \le 1 $ and $-1 \le y \le 1$ and $w(-1,0) = w(0,-1) = w(0,1) = w(1,0) = 1$ and $w(-1,-1) = w(0,0) = w(1,-1) = w(-1,1) = w(1,1) = 0$ but I don't know how to take Fourier transform from this?

How do I take the Fourier transform of my 3×3 filter kernel?

So, the $M\times N$-DFT is very well-defined to be (aside from the usual ambiguity about the factor up front)

$$
X[k,l] = \frac{1}{\sqrt{MN}}\sum_{n=0}^{N-1} \left(\sum_{m=0}^{M-1}x[m,n]\cdot e^{-jmk\frac{2\pi}{M}}\right)e^{-jnl\frac{2\pi}{N}}
$$

In your $M=N=3$ case, that becomes

$$
X[k,l] = \frac{1}{3}\sum_{n=0}^{2} \left(\sum_{m=0}^{2}x[m,n]\cdot e^{-jmk\frac{2\pi}{3}}\right)e^{-jnl\frac{2\pi}{3}} \tag1
$$

Furthermore,

$$ x = \begin{bmatrix}
0 & 1& 0 \\
1& 0 & 1\\
0 & 1&0 \\
\end{bmatrix} $$

so that all there's a lot of zeros.

Simply expand the sum from $(1)$, filling in the values of $x$.

$$\begin{array}
\,3X[k,l]&=& &\left(e^{-jk\frac{2\pi}3}\right)e^{j0} & \text{first row}\\
&& +& \left(e^{-j0}+e^{-jk\frac{4\pi}3}\right)e^{-jl\frac{2\pi}3}& \text{second row}\\
&&+ &\left(e^{-jk\frac{2\pi}3}\right)e^{-jl\frac{4\pi}3}&\text{third row}\\
&=&&e^{-jk\frac{2\pi}3}\\
&&+&e^{-jl\frac{2\pi}3} + e^{-j\frac{2\pi}3\left(l-2k\right)}\\
&&+&e^{-j\frac{2\pi}3\left(k-2l\right)}
\end{array}$$

So,

$$X[k,l]=\frac13\left(e^{-jk\frac{2\pi}3}+e^{-jl\frac{2\pi}3} + e^{-j\frac{2\pi}3\left(l-2k\right)}+e^{-j\frac{2\pi}3\left(k-2l\right)}\right)$$

Set in all values for $k$ and $l$, and get your $3\times3$ DFT output .


Reference Code for Positive Basis Pursuit Denoising I am trying to reconstruct a positive sparse signal using compressed sensing (friedlanders code), I cannot find a way to impose the positivity constraint for this implementation. I have seen some papers discussing positive signals case but cannot find any code implementations, can anyone please point out to any relevant code.

I assume you're after the following optimization problem:

$$\begin{align*}
\arg \min_{x} \; & {\left\| x \right\|}_{1} \\
\text{subject to} \; & A x = b \\
& x \succeq 0
\end{align*}$$

This is pretty simple problem if we pay attention fo the fast that given $ x \succeq 0 $ then $ {\left\| x \right\|}_{1} = \boldsymbol{1}^{T} x $.

This means the above problem is equivalent of:

$$\begin{align*}
\arg \min_{x} \; & \boldsymbol{1}^{T} x \\
\text{subject to} \; & A x = b \\
& x \succeq 0
\end{align*}$$

Now, this is a Linear Programming problem which can be easily solved using MATLAB or any other LP solver.

For instance, in MATLAB it will be something like that:

numRows = 10;
numCols = 100;

mA = randn([numRows, numCols]);
vB = randn([numRows, 1]);

vE = ones([numCols, 1]); %<! Vector of ones to sum vX

vX = linprog(vE, [], [], mA, vB, zeros([numCols, 1]), inf([numCols, 1]));


Change pitch of a wav file without changing the Sample rate? I'm creating a Python library I call PAL, (Python Audio Library) and I need some help. What ways are there to change the pitch of a wav file without changing the sampling rate? I have complete access to the data contained in said file as numbers so if there is an equation or method that would work please let me know. To recap

Change the pitch
Don't touch the sampling rate
I can mathematically change the samples freely

If you're looking for a library solution, there's a few, for instance librubberband, binaries available in many Linux distros.

If you need a wheel to reinvent, well, approaches may vary, from granular stretching (more applicable to time stretch, but these two effects always come hand in hand), to FFT uppitch (probably somewhat similar in the frequency domain), I doubt there's a standard consensus on how to implement this effect.


What the entropy equations mean? In wavelet packet image compression, different types of entropy methods can be used, like Shannon and log-energy.

Shannon entropy uses this equation $\mathrm{ent}= -\sum (x^2 \times \log(x^2))$, where log-energy equation is $\mathrm{ent}=\sum(\log x^2)$, what this equation means?
which technique can achieve better results for images (in terms of PSNR)? and how can I determine the appropriate technique?

Entropy can be used as a criterion to decide whether a wavelet packet branch should be further decomposed, or to prune a wavelet packet tree. This carries similarities with segmentation techniques like split-and-merge.

Let us take one sign convention for entropy. One hopes that entropy diminishing, locally, across some levels, meas that this part of the image is simplified by
 one wavelet decomposition, and when the entropy starts to increase (or vice-versa), further decomposition levels become useless to simplification. Following the entropy gradient is a way to find some "optimal" wavelet packet splitting.

But first,  optimal for what purpose is unclear from your question so far. Denoising, compression? Entropy can provide a nice guess to predict compression performance with an actual coder, but should be applied to binned/quantized wavelet coefficients, and may be limited unless you use higher-order/cross-entropy to take into account advanced vector coders (as compared to scalar ones).

Second, traditional wavelets packets for image are implemented separably, and notoriously, separable schemes and packets lack in efficient representation for  piecewise regular images, and induce increased aliasing, that may yield mildly efficient decompositions.

Codes implementing such decompositions are presented in Matlab, see for instance: wpdec2: Wavelet packet decomposition 2-D, and the reference paper Coifman, R.R.; M.V. Wickerhauser (1992), “Entropy-based algorithms for best basis selection,” IEEE Trans. on Inf. Theory, vol. 38, 2, pp. 713–718.


Smoothing signal / detecting bumps in a data stream (EDIT: This question follows from Extracting Binary Magnetic-Strip Card Data from raw WAV)

Here is my signal (top line) and a basic IIR filter applied (bottom line)

(EDIT: my task is to break the signal into binary 0 (frequency F) and binary 1 (frequency 2F) -- that's why it is called F2F.  So I need to process it in such a way that guarantees no false peaks.  While the screenshot makes it look trivial, there is a potential problem of getting a double peak, and also of getting false positives in the trough between real peaks.)

My question is, what methods are available for smoothing this signal?  Is IIR my best bet?

I can see at least three possibilities:

IIR y[n] = 0.9*y[n-1] + 0.1*x[n]  where y[x] =0 when x < 0
Moving / windowed average -- place a Bell curve with area 1.0 over the surrounding say w=10 samples each side and integrate bellSmooth(x) = integral[x-w,x+w] { bell(k).samp(k) }dk
Determine the expected frequency and FFT / remove higher order bins / reverse FFT

I may have answered my own question, but probably this is incomplete and I'm sure I am using the wrong terminology.  Also I can't really predict the pros and cons.   The last method is less attractive as it requires knowledge of the basic signal frequency.  But then so does the second method; I need to choose an appropriate window length.

Are there any other methods?

Effects of Averaging

Using a moving average filter will smooth out the irregularities in the signal. The noise becomes E/N where N is the length of the moving average filter. The side effect of using a MA is that the signal peaks become wider and shallower.

In addition, the frequency content of the signal will change. A moving average filter in the time domain is the same thing as convolving the frequency domain signal by a sinc function everything gets smudged out.

Peak Detection Algorithm
Peak detection is a common problem in 9/10 engineering problems. (not really, but a TON depend on them)

Typically this is what is done:

Median Thresholding

1) Look for all peaks in your signal. (i.e., a point that is larger than the two
   adjacent points
2) take this list of points and for each one of them compute:
   med_threshold = median(Peak,Width) + constantThresholmedian   where median is the
   median value of the data centered at "Peak" with Width being the number of
   points to look at.
       a) The Width(usually written as Lambda in literature) and constantThreshold
          (usually written as C) are determined by trial and error and using the ROC
          curve (Acronym below)
3) if the peak's magnitude is above this threshold accept it as a true peak.
   Else: Discard it, its a false peak
4) Generate a Receiver Operating Characteristic Curve(ROC) to how well the algorithm
   is performing.

Here is an example:

suppose we have the signal X = [ 0 0 0 0 1 3 **9** 2 1 1 **2** 1 1 ]
1) 9 and 2 are both potential peaks
2) Lets use a window of 5 and  a threshold =2
so at 9 we have [1 3 9 1 2] -> [1 1 2 3 9]  so Median(9,5) = 2
9 > 2 +2, therefor its a peak
Lets take a look at 2: [ 1 1 2 1 1] -> [1 1 1 1 2 ] Median(2,5) = 1
2 < 1+2, therefor it is NOT a peak.

Determining Frequency

Now that You've effectively found the time localization of the peak try to find their frequency:

1) Use the locations of the peaks to generate a pulse train
      a) this means create sum(Dirac_delta[t-L(n)]) where L(n) is the nth time that
      you've localized through median thresholding
2) Apply FFT Algorithm
3) Look for largest peak.

Alternate Frequency Estimation

1) Think of this like a beat in a piece of music (I learned about thresholding by
   researching Onset Detection.
2) Compute the average time distance between detected peaks.
3) now call your results BPM or PPM (pulses per minute)

Additional Avenues of Research

While you may be satisfied by the peaked signal as is, there are algorithms that are applied to a whole different beast of problems called Onset Detection.

Onset Detection is a big area in Music Information Retrieval Research. Its used to determine when a note is being played.

If you think of your tape head signal as a highly sampled signal, you can apply many of the algorithms that you would find in this paper:

http://www.elec.qmul.ac.uk/people/juan/Documents/Bello-TSAP-2005.pdf


About image compression vs. DCT I have recently started reading about image processing, my knowledge in image processing is still very limited.
What I understood that DCT can help in image compression using Huffman, JPEG or MJPEG encoding.
My question, If we have an image as below:

which means that we have only three possibilities, 1, -1 or 0 in every pixel. they are distributed randomly, what I know that every column has 4 zeros and other values either 1 or -1. they are distributed randomly.
After performing DCT transformation, Does that compression for that image can be done efficiently,  Compared with the image which is full of 1 and -1 without any zeros.  How far can we compress that image perfectly? Finally, which encoding algorithm can be used perfectly in that case ?
Thank you so much in advance.

they are distributed randomly, what I know that every column has 4 zeros and other values either 1 or -1
compression for that image can be done efficiently, Compared with the image which is full of 1 and -1 without any zeros

Since your spatial representation is already very sparse, it's likely that the DCT will reduce the compressability.
Simple calculation:
naive storage
That's a matrix of $7\cdot 8 =56$ entries, each of which takes one of 3 states, so 2 bits per pixel, so 112 bits in total.
non-zero position + sign storage
We can calculate the column entropy rather easily:

For the non-zero positions, you have 8 choose 4 different options, $\binom84=70$; so the info in that is $-\log_2\left(\frac1{70}\right)$; You can store that in 7 bit without thinking much. Just write a table of all possible non-zero positions, and number that.
If we consider the sign of the nonzero entries as bits, that's 4 bits of information.

So, the entropy in each column is $\log_2(70)+4\approx 10.13$ bit.
So, per column you need 11 bits, so with 7 columns in total 77 bits, and would, without any compression, have represented your whole image.
If we look at the 7 columns as whole, that's $70^7< 2^{43}$ possibilities of non-zero positions, i.e. we can reduce the overall storage needs to $43+7\cdot4=71$ bits. Again, no entropy coding necessary.
DCT first, entropy coding later
If you do a DCT on each column, you'll need all 8 entries; every DCT bin can take at least 5 values, so it needs 3 bit each, so 24 bit per column, so 148 bits in total.
Since the $70\cdot16$ (see above) different possible input columns are equally likely, so are the $70\cdot16$ possible DCTs; the information in seeing any one of these thus is $-\log_2\left(\frac{1}{70\cdot16}\right)\approx10.13$ bit. Since they don't have any correlation, entropy coding can, at the very best case compress your seven columns to 70.90 bits, ie. 71 bits in storage. But that case can never be achieved – that would work with an infinitely large image, not one with but 7 columns.
conclusion
for your 8×7 image:

intuitive storage format that just saves the non-zero positions + signs: 71 bits
unachievable asymptotic compression = entropy of the source: 70.9 bits

So, in this case, there's an "obvious" storage format with negligible redundancy. Thus, entropy coding (Huffmann etc) makes no sense. And, since your original image is already very sparse, any transformation, but especially the DCT, only makes compression harder.
So don't do the DCT in this synthetic example.