Does Hamming encoding preserve orthogonality? The question I have is this-- if I have two binary vectors, i.e. $x_1, x_2 \in \mathbb F_2^k$ that are orthogonal, i.e. $x_1^T x_2 = 0$, and let's say I Hamming encode each of them separately with an $(n,k)$ code to produce $y_1, y_2 \in \mathbb F_2^n$, respectively.

Can we make any general statements about the orthogonality of $y_1$ and $y_2$? Will they be orthogonal also?

I don't think you cannot state generally that $y_1$ and $y_2$ will be orthogonal. I'll try to sketch out my thinking:

Since the Hamming code is a linear code, each parity bit can be represented as a linear combination of the information bits; that is, each bit can be represented as:

$$
p_i = \mathbf{g_i^T}\mathbf{x}
$$

where $\mathbf{x}$ is the information bit vector and $\mathbf{g_i^T}$ is the generator for the $i$-th parity bit (it's one column from the code's generator matrix $\mathbf{G}$); it is a binary vector that represents the parity checks that go into $p_i$. When formulated this way, the parity bit can be thought of as the projection of the information bits onto the respective generator vector.

Assuming that you have a systematic code (all Hamming codes that I've seen before are), the encoded vectors $\mathbf{y_i}$ look like (assuming the code has $k$ information bits and $m$ parity bits):

$$
\mathbf{y_i} = \left[x_{i_1}\ x_{i_2}\ \ldots \ x_{i_k}\ |\ p_{i_1}\ p_{i_2}\ \ldots \ p_{i_m} \right]^T
$$
$$
\mathbf{y_i} = \left[x_{i_1}\ x_{i_2}\ \ldots \ x_{i_k}\ |\ \mathbf{g_1^T}\mathbf{x_i}\ \mathbf{g_2^T}\mathbf{x_i}\ \ldots \ \mathbf{g_m^T}\mathbf{x_i} \right]^T
$$

In order for the vectors $\mathbf{y_1}$ and $\mathbf{y_2}$ to be orthogonal, their dot product must be zero. It has the form:

$$
\mathbf{y_1} \cdot \mathbf{y_2} = \sum_{j=1}^k \left(x_{1_j} + x_{2_j}\right) + \sum_{j=1}^m \left(\mathbf{g_j^T}\mathbf{x_1} + \mathbf{g_j^T}\mathbf{x_2}\right)
$$

Since the information vectors $\mathbf{x_1}$ and $\mathbf{x_2}$ are orthogonal, the first sum above is zero, so we are left with:

$$
\mathbf{y_1} \cdot \mathbf{y_2} = \sum_{j=1}^m \left(\mathbf{g_j^T}\mathbf{x_1} + \mathbf{g_j^T}\mathbf{x_2}\right)
$$

$$
\mathbf{y_1} \cdot \mathbf{y_2} = \sum_{j=1}^m \mathbf{g_j^T}\left(\mathbf{x_1} + \mathbf{x_2}\right)
$$

There is no requirement that this condition be zero for an arbitrary linear block code. I haven't taken the specific structure of the Hamming code into account here, but I suspect that it doesn't happen to meet this criterion either.


How to make frequency axis for even and odd FFT length? Can I get help on how to make frequency axis going from negative frequency to positive frequency, (in Hertz), that will be the x-axis in an FFT result, but given either an even length FFT, or odd length FFT. I am having some troubles making it in MATLAB. (Assume you know the sampling frequency f_s).

You can make a positive frequency spectrum quite simply (where fs is the sampling rate and NFFT is the number of fft bins). In the Matlab implementation of the FFT algorithm, the first element is always the DC component, hence why the array starts from zero. THis is true for odd and even values of NFFT.

%//Calculate frequency axis
df = fs/NFFT;
fAxis = 0:df:(fs-df);

If you need to wrap the frequency spectrum, you need to take into account whether you have an odd numbered NFFT. There always needs to be a DC component, so . .

df = fs/NFFT;
fAxis = (0:df:(fs-df)) - (fs-mod(NFFT,2)*df)/2;

Notice how the calculation of the positive frequency axis is identical to above, but the FFT-shifted term changes to accomodate even or odd FFT lengths.

These code snippets were taken from a long answer posted on SO (that you might find interesting) found here:
https://stackoverflow.com/questions/9694297/matlab-fft-xaxis-limits-messing-up-and-fftshift/9699983#9699983


FIR to linear-phase FIR I am trying to convert an FIR impulse into linear-phase FIR impulse, I need to remove the phase response, so when it is applied to a sound it only changes the amplitude response.

Here is my code below:

clear all;
close all;
[x,fs]=audioread('impulse.wav')

s=fft(x)
ph=angle(s)

amp=abs(s)
s2=real(ifft(pol2cart(0,amp))) ; setting phase to zero

s2= circshift(s2,(round(length(s2)/2)))
plot(s2)

s2=s2/max(max(abs(s2)))

audiowrite('LinearPhaseImpulse.wav',s2,fs,'BitsPerSample',24)

Is this the correct way of converting an FIR to a linear-phase FIR impulse response?
How would one negate the phase so the FIR becomes linear-phase?

Converting an arbitrary FIR filter into a linear-phase FIR with the same magnitude response is generally impossible. As mentioned in msm's answer, this conversion must be done with an allpass filter, which must be an IIR filter, and only in those rare cases where you get appropriate pole-zero cancellations will the total filter be FIR. In general you'll end up with an IIR filter.

Having said that, there are possibilities to convert a given FIR filter into a linear-phase FIR filter with approximately the same magnitude response. By increasing the length of the linear-phase FIR filter, the approximation can be made arbitrarily good.

The method is similar to the one you implemented, but the crucial difference is zero-padding. You need to use an FFT length that is much greater than the length of the given FIR filter.

Below is an octave example explaining how it works. I chose a non-linear phase filter of length $30$, and approximated it by a length $39$ linear-phase FIR filter:

% come up with some non-linear phase FIR filter
[b,a] = butter(5,.3);
h = impz(b,a);
h = h(1:30);    % length 30
h = h(:)';

% do the FFT trick
N = 1024;       % FFT length (even!)
H = fft(h,N);
hlp = real(ifft(abs(H)));
hlp = [hlp(N/2+2:N),hlp(1:N/2)];

% windowing
Nw = 39;        % choose window length (odd!) = filter length
w = ones(1,Nw); % rect. window (choose another if you wish)
w = w(:)';
Nw2 = (Nw-1)/2;
hlp = hlp(N/2-Nw2:N/2+Nw2) .* w;

% plot
Hlp = fft(hlp,N);
n = (0:N/2)*2/N;
subplot(2,1,1), stem(h), hold on, stem(hlp,'r'), hold off
legend('NLP FIR', 'LP FIR')
subplot(2,1,2),
plot(n,20*log10(abs(H(1:N/2+1))),n,20*log10(abs(Hlp(1:N/2+1))),'r')
legend('NLP FIR', 'LP FIR') ,grid on

The plot below shows the impulse responses of the two FIR filters, and the magnitude responses. The magnitude approximation could be further improved by choosing a longer filter length for the linear-phase FIR filter.


Why is it that the input to an IFFT is first reordered in an OFDM Modulator I am back with my questions! I recently came across a paper titled "Low Latency IFFT Design for OFDM Systems Supporting Full-Duplex FDD" where they claim a new method of IFFT input mapping to reduce the output latency. However, they have mentioned how in a conventional OFDM, the null subcarriers are mapped in the middle and the data subcarriers to the side. I have attached a picture of this. But my question is why is it done in the first place? Is it something to do with the efficient use of butterfly operation?
Also, please mention a reference where it is cited that conventional OFDM systems use this method of reordering before performing an IFFT operation.
Thank you

Since Discrete Fourier Basis Vectors are $2\pi$-Periodic, hence, negative frequencies $-\frac{2\pi k}{N}$ are conventionally represented as $(2\pi - \frac{2\pi k}{N} = \frac{2\pi}{N} (N-k))$. It has nothing to do with Butterfly Structure of FFT Algorithm.

The above statement basically means that $-k^{th}$ tone is represented as $(N-k)^{th}$ tone in DFT/FFT. That is why you are seeing that DC and positive tones are placed as it is, but negative tone $-k$ is placed at $N-k$ before taking IFFT.
It is precisely the consequence of the convention that DFT coefficients are $X[k]$ where $k = 0,1,2,...,(N-1)$. $k$ goes from $0$ to $(N-1)$ and not from $(-\frac{N}{2})$ to $(\frac{N}{2} - 1)$.
So, DC + Positive $\frac{N_d}{2}$ tones are placed as first $\frac{N_d}{2}+1$ values, and then Negative $\frac{N_d}{2}$ tones are placed in reverse order from the end, i.e. $-1^{th}$ tone becomes $(N-1)^{th}$ value, $-2^{nd}$ tone becomes $(N-2)^{nd}$ value and so on to $-\frac{N_d}{2}^{th}$ tone becomes $(N-\frac{N_d}{2})^{th}$ value. The rest of the $(N - (N_d +1))$ tones are zero padded. This array of values becomes the IFFT input.


Convolution involving turning each pixel value to the most represented pixel value of the neighbourhood In order to correct gradual changes of intensities in the background of grey-scales images, I have been blurring them and then subtracting the original images from the convolved one.

In some cases, I have seen improvement using median blur instead of blur.

I was therefore wondering if I could go further and use a filter that convolves a grey-scale image so that each pixel will take the value of the most represented pixel value in the histogram of a rectangular neighbourhood.

Do you know the name of this convolution ?

Are you aware of any implementation in OpenCV ?

Thank you very much,

This is a non-linear filter, the operation being performed is not a convolution and cannot be represented by a filter kernel. It is sometimes called a mode filter, by analogy to the median filter, since the value taken by a pixel is the mode of the distribution of the neighboring values. OpenCV does median filtering (cvSmooth with CV_MEDIAN parameter), I don't think it does mode filtering though.


Calculation of actual analog input from bipolar ADC's output I have a $14$-$\textrm{bit}$ ADC which can take $-10\textrm{ V}$ to $10\textrm{ V}$ as input. So, $0$ to $16383$ (digital) represents the voltage from $-10$ to $10$. If I want to get the actual voltage from the digital values I am getting, what should I do? i.e. if $8192$ is coming then my program should identify it as $0\textrm{ V}$.

Please help me out. If my understanding is wrong, please correct me.

Your understanding is right. First subtract the bias, and then apply scaling to get analog actuality.


Where does convolution fit in DFT? I am trying to understand where does the convolution fit in Discrete Fourier Transform. I know that convolution is producing a third signal from two other signals. I also know that DFT transforms one signal from time to frequency domain.

So, what is the relationship between the two (if any)?

FFT and DFT are a faster way to implement convolution.
Since convolution in time domain corresponds to multiplication in the Fourier domain, you replace the sliding-window based implementation of convolution (in time domain) by a pointwise multiplication (in Fourier domain). This avoids you an extra inner loop over the input data and saves time.

With respect to boundary conditions, there is a minor difference with finite-length signals.
Using the convolution theorem to compute a convolution product in the Fourier domain implicitely assumes that the input data is periodic, i.e., when you reach the right end of the signal you re-enter by its leftmost part. This is also known as circular convolution.
If you need to have a convolution of explicitely non-periodic signals, you can add zeros at the beginning and end of your signal to avoid these border effects.


FIR filter delay and window length How long does it take for FIR filter to produce first output sample for given window / kernel / impulse (not sure what is best name) length?

Does it differ between minimum phase and linear phase variants? Lets say we have two FIR low pass filters,one minimum phase,other linear phase,both have 100 sample long impulse response/window length,how long will it take to obtain first output sample?

My guess is the minimum phase needs 100 sample long delay while the linear phase 50 samples.

For causal FIR linear phase filters the delay you ask for can simply be computed from the group delay and is given by (for an odd length symmetric FIR filter of length $L=2M+1$)
$$\tau = - \frac{d\phi(\omega)}{d\omega} = - \frac{d-M\omega}{d\omega} = M $$

So it will shift every input component by $M$ samples at the output.

NOTE: you can always produce your output sample $y[0]$ as soon as the first input sample $x[0]$ arrives. The delay is related with the information processing and not the sample processing;i.e., your initial $M-1$ output samples will be called as the filter transient response, before the steady output is computed.


Extending a waveform Apologies if this makes no sense or seems rather simple, I'm a Matlab newbie.

Basically, I've got a .wav audio waveform which I am resampling to various different frequencies and then trimming the length to a second (44100 samples). However, when resampling to certain tones, such as the note C3 or A5, I'm not left with enough samples to trim the note to an acceptable length.

Obviously this means that the sample needs extending, which I've worked out can be done taking a sample period s1 to s2 and adding that onto the end as many times as is necessary.

extendedSignal = x(1:s1-1);

for i = 1:4
    extendedSignal = vertcat(extendedSignal, x(s1:s2));
end

extendedSignal = vertcat(extendedSignal, x(s2:44100));

However, when it comes to finding a continuous period in which to sample, I cannot seem to find one for any note which does not sound like it is being repeated. Am I missing something really obvious?

Although it's not entirely clear to what you try to do (maybe there are better ways to accomplish your goal), I'll try to give some hints.

There exist "time stretching" techniques that can slow down a sample without changing its pitch, but they are not trivial to implement yourself. Some audio editors have plugins that support this.

As a poor man's implementation, if the sound has a single pitch, you can try something like the following (I've successfully used this in the past to extend a long, sung, single note):

First I cut the sound in three parts:

Sound -> | Begin | Middle | End |

Always cut at zero crossings.

Then, if you make "Middle" short enough, like around 0.1 seconds or so, you can repeat it a number of times times to stretch your sound a bit:

| Begin | Middle | Middle | Middle | Middle | Middle | Middle | End |


Need some references on sound quality evaluation and syncing I have to do a project where I'm receiving the same sound signal from different channels (with different offsets), but they may not be consistent all the time. At any given time, I should decide which channel is providing the best transmission, and make the connection as persistent as possible.

I'm an electronic engineering student familiar related concepts. I appreciate it if you know of any reference that would help me.

From your comment I take it you are looking for speech quality measures. As far as I know, there exist a couple of objective measures, along with subjective procedures, since quality of speech is kind of a vague concept which may not be the same to everyone.

Since your application seems to require automatic evaluation, there are a number of measures you might want to check out:

ITU-T Rec. P 835 for assessment of listening speech quality (no available reference signal)
The PESQ for listening speech quality, mentioned by Florian already (when the reference signal is available)
The PEAQ for assessment of audio quality (when the reference signal is available)
The E-model für conversational speech quality (network parameters available)

All of these are listed in this book, along with a range of subjective measures.
It seems to me that for your application, the ITU-standard is most reasonable since no clean reference signal is required.

Edit: What also may be of interest of you is the frequency-weighted segmental SNR (fwsegSNR), which incorporates the articulation index (AI). This measure apparently correlates well with subjective measures.

An overview of objective speech quality measures is also given in this Loizou paper available online.


Image Scrambling versus Image Encryption What is the difference between image scrambling and image encryption or are they the same thing in the context of image processing? What about in communication systems in general?

Scrambling is used for making a signal more white (I mean flat spectrum). Encryption is used for protecting the data from unauthorized listener. If you just scramble the data it is quite easy for a professional to descramble it and obtain the source image, but if you encrypt it then it is difficult or impossible to reverse it for someone without the code.


Filtering Wind Noise From Audio Signal I have an audio signal where most of the signal of interest is above 400 Hz but it is corrupted by heavy wind noise. What type of filter would you recommend to try to remove the wind noise considering I have audio samples of just wind but I don't have uncorrupted samples of the signal of interest. A simple high pass seems to work ok but my boss wants me to try implementing some type of adaptive or spectral subtraction type filter.

Laminar air flow doesn't produce much in the way of sound.  Turbulence, which is usually excited by interaction that is, the disrupted flow by an interacting object.  Some of what is heard, is radiation from monopoles (dipoles, and quadrupoles) within the flow and some are resonances related to the object doing the disruption.  You can often look at a tree in high wind and visually correlate the intensity of sound with the intensity of motion of leaves and limbs.  If you stand near a steel tower in high wind, your will hear it "sing"

Part of what you record is the broadband excited vibration of the microphone's full structure.  This is why a windscreen on microphone works.  The screen dampens the wind incident on the microphones structure and is a poor radiator.  One model is a random stream of impulses exciting a linear filter.  The mic structure gives color to the wind sound.

This is a situation where a little care in data collection saves a lot of trouble down the line.

I read someplace once "Never invent gloves to fix a leaky valve"

You might also want to try to whiten the wind noise and see what happens.  It would be similar to inverse filtering on the structure resonance.

Another characteristic of the vehicle sounds are the doppler shifts at the time when the car is at closest distance to the microphone.


Why does this image look so bad when dithered with a certain palette? I'm trying to write software to create a mosaic of one big image made up of a bunch of smaller images. The result is okay, but I think it would be better if I used dithering.

My problem is that dithering is giving even worse results and I don't know why.

Here's the original:

Here's the same image quantized (no dithering) with the following palette (in (R,G,B) format):

(162, 143, 66)
(128, 148, 168)
(120, 99, 100)
(31, 39, 97)
(126, 116, 103)
(203, 35, 9)
(57, 81, 43)
(104, 101, 98)

Not great, but workable. About what I'd expect.

Now here's the same image dithered using Floyd-Steinberg:

Pretty much unrecognizable. What's going on here?

To make sure my dithering algorithm was implemented properly, here's the same image dithered using 2-color and 8-color palettes.

Both of these are a whole lot better than the other one. Why is that?

I think what's going on is something analogous to CCD sensor blooming.

The badly working palette does not cover the extreme color values of the original image. Floyd–Steinberg dithering works by error diffusion, where the error in color at a pixel being processed is spread to nearby pixels not yet processed, to the right and to some of the pixels on the next pixel row. With the kitten, you have near-black and near-white areas, and in those areas the algorithm – working like a bulldozer – will accumulate more and more error because the colors it can use will always fall short the same way, and when it finally comes out of the problem area, the accumulated error is so large that it takes many pixels to get rid of it, smudging the image.

In other words, if the image has out-of-gamut colors, you get problems. If a color is inside any polyhedron formed by any of the palette colors, then it is in gamut and can be reproduced by alternating colors from the palette. In other words, the color is in the convex hull of the palette colors. If it's not in the convex hull of the palette colors, it is out of gamut. Gamut mapping before dithering would resolve the problem.


How to find the output mean and autocorrelation of a non-linear system I need help with this question. I am sure this is the right StackExchange forum for this type of question.
Consider a nonlinear device such that the output is $Y(t) = aX^2(t)$, where the input X(t) consists of a signal plus a noise component, $X(t) = S(t) + N(t)$.
Determine the mean and autocorrelation function for $Y(t)$ when the signal $S(t)$ and the noise $N(t)$ are both Gaussian random processes and wide sense stationary (WSS) with zero mean, and $S(t)$ is independent of $N(t)$.
I know that for a WSS signal, the mean and autocorrelation depends on the time difference, $\tau$ only. Meaning that the mean of $S(t)$ and $N(t)$ is a constant and their individual autocorrelation is not dependent on $t$.
I've been working only with linear systems. I don't know how to solve this problem since $Y(t)$ a nonlinear system. I will be happy if you can help me with the solution
UPDATE 1
I followed the suggestion from @Hilmar and came up with the following:
1. The mean of $Y(t) = 0$ since $S(t)$ and $N(t)$ have 0 mean.
2. For the autocorrelation, I used this formula $R_{yy}(\tau)=R_{xx}(\tau)*h(\tau)*h(-\tau)$.
looking for $R_{xx}$, I used $R_{xx}=E[X(t)X(t+\tau)]$.
I ended up having $R_{xx} = E[S(t)S(t+\tau) +S(t)N(t+\tau)+N(t)S(t+\tau) + N(t)N(t+\tau)]$. Since $S(t)$ and $N(t)$ are independent, $R_{SN} = R_{NS} = 0$
$=> R_{XX}= E[S(t)S(t+\tau) + N(t)N(t+\tau)]$
$=> R_{XX} = R_{SS} + R_{NN}$

$R_{YY} = (R_{SS} + R_{NN})*h(\tau)*h(-\tau)$
Now I am stuck. From the definition of WSS signal, $R_{SS}$ and $R_{NN}$ are constants. i.e. they are dependent only on $\tau$. I can assume their sum to a some constant value, $C$I don't know how to move from here. I will need to find $h(\tau)$ I don't know how to go about it.

Put your second equations into your first equation, express $Y(t)$ as a function of $S(t)$ and $N(t)$
Apply the definition for mean and autocorrelation.
Simplify and solve


Calculate Horizontal Pixel Shift for Super Resolution I have a sample that is moving at a constant rate which is known.  Since I know the distance traveled I can equate that with the number of pixels for the pixels shift.  The only way I can think to do this is to do a camera calibration but I am working on a small scale and was wondering if there is another way to calculate the pixel shift for each image?

So you have an item in an image that is moving at a known constant speed in real units, but you want to know how to convert that speed into pixels per second (or whatever time unit)? If that's the case, then you can calculate it based on measurements - have the item move from one side of the image to the other, and count the number of pixels it has moved and in how long, and you have your rate.


bibliography about implementation of IIR filters I need to implement a 6th order IIR filter. The filter I need to implement is the A-weighting. I've simply used the bilinear function of Matlab with an fs of 20Khz. I need to implement this over the LPC1769 micro-controller. The application is mainly the filter and the communication with the ADC and the LCD to show the RMS value of the signal after being processed by the filter.

I've been reading the book "Real-Time Digital Signal Processing..." of Kuo, but I need more detail in the implemetation considerations.Could you recommend me some books or papers about this?

I would like to know more about how to implement filters, the implications of the ordering choise of the biquad stages, or the meaning of the parameters that are used in the fdatool of matlab at the moment of quantizing a filter.

Okay, according to the wikipedia article you cited (with an appropriate substitution in notation):

$$ \Big|H(j2\pi f)\Big|^2= 10^{\frac{2.00}{10}} \times \frac{12194^4\cdot f^8}{\big(f^2+20.6^2 \big)^2(f^2+107.7^2)(f^2+737.9^2)\big(f^2+12194^2\big)^2}$$

and the dB amplitude curve is this:

$$A(f)=10\log_{10}\left(\Big|H(j2\pi f)\Big|^2\right) $$

I hope you can see how I modified the language and that it is equivalent.  The purpose is to express this transfer function as an $H(s)$ and get where (and how many) the poles as zeros are.

Now you gotta factor some more shit out.  Oh, this is laborious and painful.  Sorry, unless someone wants to finish this answer, I will have to return to it.


How to measure a lung nodule in a CT Scan DICOM image? In this question, I would like to focus on intensity value of the CT Scan. First, take a look at the image below:

The upper image is the original image, while the lower image is the thresholded version. In order to measure volume of any shape, theoretically, it is possible to simply just count the number of voxels in the image. However, the outermost layer of the object (ex. nodule) shows darker intensity, while all voxels inside the object have very high intensity. If I simply count the voxels in the thresholded version, I will highly likely to get a result volume larger than the actual volume for the lung nodule.

I also see that there are variables such as window center (level) and window width, which can be used to adjust the DICOM image's intensity information. Different intensity can alter the result volume.

So here is the question: if I am to measure any given lung nodule, what should I do in order to achieve the best possible precision? When should we ignore the lower intensity voxels? Or must I do this in some other ways?

Aside from the whole discussion of pure signal processing:
What exactly do you define as "the nodule". This is usually a biological entity with difficult to define borders. The nature of the nodule is sometimes invasive growth and therefore even in histological sections poorly defined. The CT itself has therefore a threshold higher than histology and therefore the true border of a nodule cannot defined easily. On the other hand, healthy tissue around the nodule might be compressed and appear dense in CT. This might be depending on the respiration phase the scan was triggert on (best results in midinspiration, according to newer data). Or inflammation might obscure the true border of the nodule.

Another aspect is interpolation of the images. The technique is usually a spiral-CT, so you don´t miss any lesions. The sections are therefore calculated. This results less defined borders of objects. If the focus would be on border detection or near histology resolution, high-resolution scans are usually required. This is what you do for interstitial lung disease.
Unfortunately, they slice through the lung with a pretty large distance in between. Radiation dose would be extremely high, if you try to get a "full volume scan". But in this type of scan, you have to be aware, that You´s miss lesions between slices.

To come back to the initial question: I think you will have to validate your technique - whatever it is - with the gold standard. Which is histology sections. (Unfortunately, the lung isn´t easily sectioned...). Another option is an additional technique: like PET-CT (combining positron emission tomography with computer tomography), but the alignment algorithms are sometimes tricky.


Continuous time-domain electrical signal that only has a finite number of possible values I have recently started learning DSP and the book that I am reading challenges the readers to find a:

  continuous time-domain electrical signal that only has a finite number
  of amplitude values

I obviously could not think of an example and hence I am here.

Please suggest an example of continuous electrical signal that takes only finite number of discrete values ?

Well, a DC signal for sure.

Other than that, there is probably more context from the book required (especially about the assumptions).


Measure scalloping loss of window function Is there an easy/general way to calculate this? For instance, the rectangular window has a maximum scalloping loss of 0.3634, which can be derived in frequency domain from a normalized sinc(0.5).

The only way, given a precomputed window function, I can think of, is to apply the window function to a signal whose period is 2pi*(k + 0.5)/N and amplitude A, calculate a N FFT and measure the bin k. Bin k should have a magnitude of A * N * f, where f is a factor giving the worst-case scalloping loss, irrespective of N/sampling size, right?

Harris defines the scalloping loss as the ratio of coherent gain for a tone located half a bin from a DFT sample to the coherent gain at a sample (see the paper for the definitions)
$$
\mathrm{SL} = \frac{\left\lvert \sum_n w(nT)e^{-j\frac{\pi}{N}n} \right\rvert}{\sum_n w(nT)} = \frac{W\left(\omega_{1/2}\right)}{W(0)}
$$
Where $w$ is the window and $W$ is the window DFT

$$
\omega_{1/2} =\frac{\pi}{NT}
$$

Guess the value you are after is really $1-SL$

EDIT: I believe the left hand side of the equation is enough to answer @Shaggi's question, but to clarify a bit the notation of the rest, Harris defines DFTs as functions of $\omega_k$, so

\begin{align}
  W(\omega_k) & = W\left(\frac{2\pi k}{NT}\right) = \\
              & = \sum_n w(nT) e^{-j\omega_k nT}\\
              & = \sum_n w(nT) e^{-j 2\pi kn/N}\\
  W(\omega_{1/2}) & = \sum_n w(nT) e^{-j \pi n/N}\\
  W(0) &= \sum_n w(nT)
\end{align}

  F. J. Harris - On the Use of Windows for Harmonic Analysis with the
  Discrete Fourier
  Transform


Confusion with regards to STFT and MFCCs I'm wondering if someone can clear the following up for me as I'm having a difficult time understanding the concept.

I'm looking to implement and learn about the processes and uses using short-time Fourier transforms. I have had experience using DFT's and wanted to expand my knowledge.

From research that I have read online, as well as in some text books I understand that the STFT is calculated by the following:

Take the inputted signal of $x[n]$
Multiply $x[n]$ by a window function (e.g. Hamming, Hanning, Triangular, .., Barlett) by sliding the block size over the time axis
Compute the DFT for the resulting values of step 2

Is this correct? I know the following steps above are for usually a continuous-time STFT and in the case of a Discrete-time STFT the data can be broken up into chunks.

The other thing that is confusing me:

For a project last year, I had to implement a MFCC algorithm. Inside this algorithm I computed the "Mel Scale" in Triangular filter bank and then multiplied each of these values against the FFT (DFT) computed on the inputted signal in Discrete-Time. Is this therefore a use of the STFT, just using the Mel-Frame over Hamming etc.. Or, would I still need to multiply the result of the FFT (DFT) with the Triangular filter bank if I was to use a STFT in order to calculate the MFCC values?

Thank yo for reading.

I am going to break down this answer into two main parts, one for a step-by-step computation of the STFT, and the other in regards to how to compute the MFCC's, once you have the STFT.

STFT: Let us say you wanted to compute the STFT of your signal $x[n]$, of total length $N=100$, and block size $M = 10$, with $50$% overlap, and using a hamming window. I will give you a step by step, (and redundant) flow, to get you started:

Take the first $M=10$ samples out, that is, samples $x[0], x[1], x[2] ...  x[9]$.
Window them, by multiplying them, element by element, with your window coefficients.
Take the DFT of this result.
Take the next $M=10$ samples out, that is, samples $x[5], x[6], x[7], ... x[14]$.
Window them, by multiplying them, element by element, with your window coefficients.
Take the DFT of this result.

etc, etc. You get the idea. Let us say that you took a (same size) DFT in the above example. Then, you will end up with a $10$ x $19$ complex matrix. If you then simply take the absolute magnitude of the above matrix, you will end up with a real matrix, corresponding to the absolute magnitude of your complex matrix. This is called the Short-Time-Fourier-Transform, or STFT. We will use the STFT as the input into computation of the MFCCs.

Before that, let us get some formulas and terminology straight: Let us call our STFT matrix $S$. This matrix is composed of $B=19$ blocks, or columns, and $K=10$ frequency bands, or rows. To refer to a time-frequency point in $S$, we will say $S[k,b]$, where $b \in 1, 2, ... B$, temporal blocks, (columns) and $k \in 1, 2, ... K$ sub-bands, (rows).

MFCC:

Once you have the STFT computed, you can go ahead and use that as a stepping stone for computing the MFCC's. Regarding your question, you seem to be confusing the time domain windowing done on each block for the STFT, and the frequency domain windowing done in the MFCC. They are completely separate and unrelated.

Let us say that you computed your STFT using a hamming window, and you have a $K=500$ x $B=1000$ real matrix. Thus using our prior terminology, you have a real matrix $S_{K \text x B}$.  Thus, you have $500$ sub-bands, (frequencies), and $1000$ blocks, (time).

Let $Y_{b}[k] = |S[k,b]|^2$. That is, we are letting $Y_{b}[k]$ represent the power-spectral density of a particular column/block of the STFT matrix $S$. Obviously, each $Y_{b}[k]$ is a $500$ x $1$ column vector.

Clearly, we have $1000$ $Y$'s. That is, we have $Y_{1}[k], Y_{2}[k], Y_{3}[k], ... Y_{1000}[k]$.

Now, you want to compute the MFCC of each column. Take the first column. This is $Y_{1}[k]$. Forget about all the other columns, and forget about the hamming window. The first column is now used in your MFCC algorithm. The first step in MFCC was getting a PSD column, which you have, because we have $Y_{1}[k]$.

The second step is the computation of, say, 10 power estimates. This is computed by windowing this column, by 10 overlapping windows, as per the mel-scale, like so:

Note how each (triangular) window is mostly zeros. So now you literally do a dot product of each window, with $Y_{1}[k]$. You will get $10$ numbers, for this column. Now take the log of those $10$ numbers. Now take the DCT of this result. Congratulations! You just computed the MFCC for this column. You just computed the MFCCs for $Y_{1}[k]$.

Now, go back, and repeat for all the other $Y_{2}[k]$. Then repeat for $Y_{3}[k]$, etc.


Is there a common name for the first order IIR averaging filter? Is there a term to describe this IIR averaging/smoothing filter?

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

This is often called a leaky integrator, a special case of a first-order IIR lowpass filter. They are discussed in more detail in several previous questions:

Is a leaky integrator the same thing as a low pass filter?
Is there a technical term for this simple method of smoothing out a signal?
Designing a practical integrator
"Low-pass filter" in non-EE, software API contexts


What are the statistics of the discrete Fourier transform of white Gaussian noise? Consider a white Gaussian noise signal $ x \left( t \right) $.
If we sample this signal and compute the discrete Fourier transform, what are the statistics of the resulting Fourier amplitudes?

Math tools

We can do the calculation using some basic elements of probability theory and Fourier analysis.
There are three elements (we denote the probability density of a random variable $X$ at value $x$ as $P_X(x)$):

Given a random variable $X$ with distribution $P_X(x)$, the distribution of the scaled variable $Y = aX$ is $P_Y(y) = (1/a)P_X(y/a)$.
The probability distribution of a sum of two random variables is equal to the convolution of the probability distributions of the summands. In other words, if $Z = X + Y$ then $P_Z(z) = (P_X \otimes P_Y)(z)$ where $\otimes$ indicates convolution.
The Fourier transform of the convolution of two functions is equal to the product of the Fourier transforms of those two functions. In other words:

$$\int dx \, (f \otimes g)(x) e^{-i k x} =
\left( \int dx \, f(x) e^{-ikx} \right)
\left( \int dx \, g(x) e^{-ikx} \right) \, .
$$

Calculation

Denote the random process as $x(t)$.
Discrete sampling produces a sequence of values $x_n$ which we assume to be statistically uncorrelated.
We also assume that for each $n$ $x_n$ is Gaussian distributed with standard deviation $\sigma$.
We denote the Gaussian function with standard deviation $\sigma$ by the symbol $G_\sigma$ so we would say that $P_{x_n}(x) = G_{\sigma}(x)$.

The discrete Fourier transform amplitudes are defined as
$$X_k \equiv \sum_{n=0}^{N-1} x_n e^{-i 2 \pi n k /N} \, .$$
Focusing for now on just the real part we have
$$\Re X_k = \sum_{n=0}^{N-1} x_n \cos(2 \pi n k /N) \, .$$
This is just a sum, so by rule #2 the probability distribution of $\Re X_k$ is equal to the multiple convolution of the probability distributions of the terms being summed.
We rewrite the sum as
$$\Re X_k = \sum_{n=0}^{N-1} y_n$$
where
$$y_n \equiv x_n \cos(2\pi n k / N) \, .$$
The cosine factor is a deterministic scale factor.
We know that the distribution of $x_n$ is $G_\sigma$ so we can use rule #1 from above to write the distribution of $y_n$:
$$P_{y_n}(y) = \frac{1}{\cos(2 \pi n k / N)}G_\sigma \left( \frac{y}{\cos(2 \pi n k / N)} \right) = G_{\sigma c_{n,k}}(y)$$
where for brevity of notation we've defined $c_{n,k} \equiv \cos(2\pi n k / N)$.

Therefore, the distribution of $\Re X_k$ is the multiple convolution over the functions $G_{\sigma c_{n,k}}$:
$$P_{\Re X_k}(x) = \left( G_{\sigma c_{0,k}} \otimes G_{\sigma c_{1,k}} \otimes \cdots \right)(x) \, .$$

It's not obvious how to do the multiple convolution, but using rule #3 it's easy.
Denoting the Fourier transform of a function by $\mathcal{F}$ we have
$$\mathcal{F}(P_{\Re X_k}) = \prod_{n=0}^{N-1} \mathcal{F}(G_{\sigma c_{n,k}}) \, .$$

The Fourier transform of a Gaussian with width $\sigma$ is another Gaussian with width $1/\sigma$, so we get
\begin{align}
\mathcal{F}(P_{\Re X_k})(\nu)
&= \prod_{n=0}^{N-1} G_{1/\sigma c_{n,k}}\\
&= \prod_{n=0}^{N-1} \sqrt{\frac{\sigma^2 c_{n,k}^2}{2\pi}} \exp \left[ \frac{-\nu^2}{2 (1 / \sigma^2 c_{n,k}^2)}\right] \\
&= \left( \frac{\sigma^2}{2\pi} \right)^{N/2} \left(\prod_{n=0}^{N-1}c_{n,k} \right) \exp \left[ -\frac{\nu^2}{2} \sigma^2 \sum_{n=0}^{N-1} \cos(2\pi nk/N)^2 \right] \, .
\end{align}
All of the stuff preceding the exponential are independent of $\nu$ and are therefore normalization factors, so we ignore them.
The sum is just $N/2$ so we get
\begin{align}
\mathcal{F}(P_{\Re X_k}) &\propto
\exp \left[ - \frac{\nu^2}{2} \sigma^2 \frac{N}{2}\right]\\
&= G_{\sqrt{2 / \sigma^2 N}}
\end{align}
and therefore
$$P_{\Re X_k} = G_{\sigma \sqrt{N/2}} \, .$$

We have therefore computed the probability distribution of the real part of the Fourier coefficient $X_k$.
It is Gaussian distributed with standard deviation $\sigma \sqrt{N/2}$.
Note that the distribution is independent of the frequency index $k$, which makes sense for uncorrelated noise.
By symmetry the imaginary part should be distributed exactly the same.

Intuitively we expect adding more integration should reduce the width of the resulting noise distribution.
However, we found that the standard deviation of the distribution of $X_k$ grows as $\sqrt{N}$.
This is just due to our choice of normalization of the discrete Fourier transform.
If we had instead normalized it like this
$$X_k = \frac{1}{N} \sum_{n=0}^{N-1} x_n e^{-i 2 \pi n k /N}$$
then we would have found
$$P_{\Re X_k} = G_{\sigma / \sqrt{2N}}$$
which agrees with the intuition that the noise distribution gets smaller as we add more data.
With this normalization, a coherent signal would demodulate to a fixed amplitude phasor, so we recover the usual relation that the ratio of the signal to noise amplitudes scales as $\sqrt{N}$.


how I can use the FFT data to compute the convolution? First,I have a vector A of length N,through process

A1 = fft(A,N);
A2 = Filter(A1);
A3 = ifft(A2,N);

Now,I need compute the convolution of vector A2 with vector B,vector B is another vector with N elements，I noticed that seemed that I could use

ifft(fft(A3,2N).*fft(B,2N))

to calculate convolution. fft(A3,2N) means padding A3 with N zeros.

In fact, A2 = fft(A3,N)

So, the problem is what can I do to avoid calculate A3 and fft(A3) by using A2?

Or maybe the problem is how can I use A2 to calculate fft(A3,2N)?I wan to know if I can reduce some calculations in the whole process.thanks.
Here is a matlab test

clear;
N = 4;
A = 1:4;
B = 3:6;
A1 = fft(A);
A2 = (A1).*2020+1234;%filter(no sense)
A3 = ifft(A2);
A4 = fft(A3,2*N);
B2 = fft(B,2*N);
x = xcorr(A3,B)%ans: 0.1952    0.4051    0.6958    1.0470    0.7676    0.5050    0.2424
x2 = ifft(A4.*conj(B2),2*N)%ans: 1.0470    0.7676    0.5050    0.2424         0    0.1952    0.4051    0.6958
x3 = ifft([A2(1:N/2) A2(N/2+1)/2 zeros(1, 2*N-length(A2)-1) A2(N/2+1)/2 A2(N/2+2:end)].*conj(B2),2*N)
%ans: 3.5624    4.2973    5.4391    6.3191    6.2232    5.2077    4.0659    3.4665

You can pad zeros in frequency domain. But you need to take care about Nyquist bin when N is even.

A3 = [A2(1:N/2) A2(N/2+1)/2 zeros(1, 2*N-length(A2)-1) A2(N/2+1)/2 A2(N/2+2:end)];

Here Nyquist bin is A2(N/2+1)/2 and zeros are padded in the middle of the spectrum

For N-odd

A3 = [A2(1:(N-1)/2+1) zeros(1, 2*N-length(A2)) A2((N-1)/2+2:end)];

This is only for the case if you need to interpolate in time-domain. Almost irrelevant to the problem you have.

Posted the solution to you problem below.


Why/How a Gaussian Channel is both band-limited and Power limited? All physical channels are band-limited, such that,

  X(f) = 0 for |f| > W

A Gaussian Channel is limited in both power and bandwidth. For this reason we encounter phrases, such as, "A band-limited Gaussian Channel" (for e.g., in the opening statement of Channel Capacity theorem).

I have two questions:

Why/How a Gaussian Channel is band-limited and Power limited?
What are physical implications of using such a channel?

Please help/guide me to understand this.

What is usually meant by a Gaussian channel is a channel with additive (white) Gaussian noise (AWGN), which simply means that the received signal will equal the transmitted signal plus some unknown noise with a certain power and a Gaussian probability density function. A pure AWGN channel is neither band-limited nor power-limited. However, these restrictions are present in real-world systems, so it is usually assumed that the channel is a band-limited AWGN channel. Furthermore, it is normally assumed that the transmitter can transmit with some given maximum power. This latter fact will result in a bound on the achievable signal-to-noise ratio (SNR).

So, to answer your question, the finite band-width and the power limitation are realistic constraints, but they have nothing to do with the AWGN channel, which only specifies the type of noise (Gaussian) and how it influences the signal (by addition).


Equation Derivation: Help with intermediate steps As per the attached image above, I don't know how the equation in the first line transformed to the form in second line.

I tried solving it by multiplication but the result was different.

This is not directly related to signal processing, but similar calculations are ubiquitous in DSP. So, here we go:

$$\begin{align}e^{-n/D}\frac{1-e^{(n+1)/D}}{1-e^{1/D}}&=\frac{e^{-n/D}-e^{1/D}}{1-e^{1/D}}\\&=\frac{e^{1/D}}{e^{1/D}}\cdot\frac{e^{-(n+1)/D}-1}{e^{-1/D}-1}\\&=\frac{1-e^{-(n+1)/D}}{1-e^{-1/D}}\end{align}$$


Why does DFT have only $N$ components? Why does the DFT have only $N$ components in it? I can see that after N components the frequency component is periodic and repeats with the same values but that does not seem to explain why we can have only $N$ components. Is this related to the sampling theorem? How does this all tie together?

The DFT is discrete in both time and frequency domains with the same number of samples in each domain, this means the input to the transform and the output of the transform are both discrete and both have $N$ samples. This is defined in the formula for the DFT where we have $N$ samples in time indexed as $n = 0$ to $N-1$ and $N$ samples in frequency indexed as $k = 0$ to $N-1$. The sampling frequency is located at $k = N$.
$$  X[k] = \sum_{n=0}^{N-1}x[n]e^{-j k \frac{2\pi}{N} n}$$
This is by definition, as for example the DTFT is continuous in frequency given by
$$  X(e^{j\omega}) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}$$
As we see in the DFT formula, each sample $k$ in $X[k]$ is a correlation (sum of products) of the time domain waveform $x[n]$ with the frequency tone given by $e^{-j k \frac{2\pi}{N} n}$ -- observe that the form $e^{j\phi}$ is a phasor with magnitude $1$ and angle $\phi$, so we see that for each value of $k$ that we pick, $e^{-j k \frac{2\pi}{N} n}$ is just a phasor spinning at a different rate defined by $k$. This causes the correlation to maximize when $k$ is equal to or close to the frequency of $x[n]$.
The reason for $N$ samples is also quite intuitive since as the OP saw the transform mathematically behaves the same way as if the time domain waveform was periodic: This is similar to the Fourier Series Expansion (FSE) which decomposes any single valued analytic function into an infinite series of discrete frequencies, with each frequency an integer multiple of the fundamental, with the fundamental given by $f=1/T$ Hz where $T$ is the duration in seconds of the time domain signal: That makes sense as the FSE also behaves as if that time domain waveform continued to infinity as a periodic waveform; if we repeat an arbitrary waveform over $T$ seconds, it will have a fundamental frequency at $1/T$ Hz. It also must be discrete, since each of the individual components given by the higher harmonics will also repeat only over duration $T$, and we recreate the waveform by summing each of these component (so any component that didn't repeat over $T$ would violate our ability to do that).
Thus knowing we have $N$ samples in time, the fundamental frequency would then be at $1/N$.  We also know that the sampling rate is at bin $k = N$ in frequency, thus the fundamental frequency is at $k=1$, and there must then be $N$ total frequency samples (some may be equal to 0), each one the integer multiple of the fundamental, just like the FSE!


DCT vs DFT why do we need/want phase? Why do we use complex exponentials in the Fourier transform, why do we want the phase part? As opposed to in DCT where there is no phase and only magnitude?

Moreover, what does this phase concretely mean for the signal reconstruction/decomposition? (I know it is related to a shift of a signal, i.e. i know that sin and cos are phase shifted by pi/2 rad u.e. 90 degrees but i can t really make sense of it in Fourier?

A DCT is identical to a DFT, twice as long, of the input data concatenated with its mirror image.   Data concatenated with its mirror results in symmetric data.

Since a symmetric vector is strictly even, there are no odd (sine or imaginary) components in the DFT result. Just cosine (or real or even) components.  Thus you can reconstruct the original input without the all-zero imaginary component.  (In practice, there might be some numerical rounding noise).

Not so with a DFT of just the original data (not doubled and mirrored), which likely is not exactly symmetric (even).  Thus you need the full complex result of a DFT to reconstruct its input (unless it just happens to be perfectly circularly symmetric around the 1st vector element, x(0)).


Can channel capcity be explained in simple terms? Can channel capacity be explained in simple terms - without mutual information and such probabilistic concepts? In its most general form, what are all the parameters that it depends on?  What is its dependence on physical parameters such as the electrical parameters if we talk of information conveying through electrical signals?  Can we measure/calculate the channel capacity of any channel given its physical parameters?

There is a limit where simplifying a scientific concept starts becoming inaccurate.

However, the capacity of a channel is a really simple idea that is completely illustrated by the following equation:

$C = B\log_2(1+\frac{S}{N})$

The equation says that the maximum amount of information (in bits per second) that can be transmitted in a channel error-free is proportional to the Bandwidth $B$ in Hertz, and is also proportional to the logarithm of base-2 of 1+ the signal to noise ratio.

So let us break this down to two points, the bandwidth and signal to noise ratio.

The bandwidth is the range of frequencies that you use in your channel.  Let us take an example where your channel has a really small bandwidth ($B=0$) where you only have one frequency that you can operate on. This translates to a single sinusoid of a certain frequency. Note that you cannot turn it on or off, or even change its phase because that will result in a range of other frequencies. Therefore communication is impossible on this channel.

The signal to noise ratio, is a measure of how strong a signal is to how strong the noise is and can take values between 0 and $\infty$. A signal power equal to 0 or a noise power equal to $\inf$ will make $\log_2(1+S/N)$ go to zero where communication is impossible.

On the other hand, 0 noise power, or an infinite signal power (given the bandwidth is not 0) will mean that you can send as much data as you like through the channel without any errors. Which is also the case for an infinite bandwidth and a non-zero value of $\log_2(1+S/N)$.


getting frequency response from input and output signal How does one get the frequency response of a filter, given an input signal and the signal output by the filter?

What you are looking for is probably an empirical transfer-function estimate (ETFE). Several methods can be used, but Welch's averaged periodogram [1] is arguably the most used [2].

The plant output data is usually generated using Gaussian white-noise excitation, although more informative input signals can be generated by experiment design, if prior information about the plant is known [3]. The ETFE of the plant $\widehat{G}(k)$ is found as the quotient of the cross power spectral density estimate of the input and the measured output $P_{yu}(k)$, and the power spectral density estimate of the input $P_{uu}(k)$, i.e.,
\begin{equation*}
 \widehat{G}(k) = \frac{P_{yu}(k)}{P_{uu}(k)} .
\end{equation*}
In Welch's method, the time-series data is divided into windowed segments, with an option to use overlapping segments. Then, a modified periodogram of each segment is computed and the results are averaged. Welch's method for generating an ETFE corresponds to the function tfestimate in MATLAB. One of the advantages of Welch's method is the flexibility in terms of the number of frequency samples and excitation signal used.


Removing unwanted peaks from signal I have a signal which looks something like:

The 12 main round peaks is what I'd like to keep. However, the peak detection algorithm I apply is picking up the in-between noisy peaks. Any strategies I could use to get rid of this noise?

I tried smoothing the signal, but the noisy peaks are still there.

Thanks

Maybe you could tell us more about the peak picking algorithm you are using!

Some ideas:

Use median filtering to remove noise (rather than a linear filter).
If you have prior knowledge of the shape of the peaks you want to detect, use a correlator (matched filter) for this shape.
Discriminate peaks according their amplitude.


When can we write Heisenberg uncertainty Principle as a equality? We know that Heisenberg uncertainty Principle states that
$$\Delta f \Delta t \geq \frac{1}{4 \pi}.$$

But (in many case for Morlet wavelet) I have seen that they changed the inequality to an equality. Now my question is when are we allowed to change the inequality to an equality:
$$\Delta f \Delta t = \frac{1}{4 \pi} $$
why =

It is important to define the time and frequency widths $\Delta_t$ and $\Delta_{\omega}$ of a signal before discussing any special forms of the uncertainty principle. There is no unique definition of these quantities. With appropriate definitions it can be shown that only the Gaussian signal satisfies the uncertainty principle with equality.

Consider a signal $f(t)$ with Fourier transform $F(\omega)$ satisfying

$$\int_{-\infty}^{\infty}f^2(t)dt=1\quad\textrm{(unit energy)}\\
\int_{-\infty}^{\infty}t|f(t)|^2dt=0\quad\textrm{(centered around }t=0)\\
\int_{-\infty}^{\infty}\omega|F(\omega)|^2d\omega=0\quad\textrm{(centered around }\omega=0)$$

None of these conditions is actually a restriction. They can all be satisfied (for signals with finite energy) by appropriate scaling, translation and modulation.

If we now define time and frequency widths as follows

$$\Delta_t^2=\int_{-\infty}^{\infty}t^2|f(t)|^2dt\\
\Delta_{\omega}^2=\int_{-\infty}^{\infty}\omega^2|F(\omega)|^2d\omega$$

then the uncertainty principle states that

$$\Delta_t^2\Delta_{\omega}^2\ge\frac{\pi}{2}\tag{2.6.2}$$

(if $f(t)$ vanishes faster than $1/\sqrt{t}$ for $t\rightarrow\pm\infty$)

where the inequality is satisfied with equality for the Gaussian signal

$$f(t)=\sqrt{\frac{\alpha}{\pi}}e^{-\alpha t^2}\tag{2.6.3}$$

The equation numbers above correspond to the proof below which is from Wavelets and Subband Coding by Vetterli and Kovacevic (p.80):


What is the unit of SNR? SNR is defined as:

$$SNR=10 \log_{10}\left(\frac{\sigma_s^2}{\sigma_n^2}\right)$$

For instance, when we say the SNR of 10 dB is 10 times higher than that of 20 dB, it is higher in what unit? Amplitude, power or what?

Also, could you relate as to why it's 10 times higher for the above example.

SNR stands for Signal to Noise Ratio. It is a ratio and as such does not have any units, it describes the proportion of signal to undesired noise. There is no single correct measure of SNR, it differs depending on the application.

In the equation you have given, the SNR is broken down in the following way:

1) Calculate the power ratio

$$\frac{\sigma_s^2}{\sigma_n^2}$$

Your equation requires us to know the variance of the signal and the noise, dividing one by the other gives the power ratio. To answer one of your questions, this is the ratio of power, not amplitude, we could achieve the same result by taking a ratio of the root mean square of the signal and noise - another measure of power.

2) Express the ratio in decibels

$$10\log_{10}$$

This is simply a convenience, the SNR may be a very large number so expressing the value on a logarithmic scale can make the information more manageable. Converting a unit to decibels simply involves taking the log10 of the value and multiplying by 10.


Stability of a system in time-domain Consider a system for which the input $x(t)$ and output $y(t)$ are related by the differential equation $$\frac{d^2y(t)}{dt^2} + \frac{3}{2}\frac{dy(t)}{dt} - y(t) = x(t) \tag{1}$$Determine $h(t)$ for each of the following cases:

The system is stable.
The system is causal.

My attempt: I know that this question can be solved easily using the Laplace transform but I'm trying to solve it in time-domain. Since it's assumed that system is linear, causality is equivalent to:
For any time $t_0$ and any input $x(t)$ such that $x(t) = 0$ for $t \lt t_0$, the corresponding output $y(t)$ must also be zero for $t < t_0$. This is the initial rest condition and we can easily solve $(1)$ by methods like this.
Stability for an LTI system is equivalent to $$\int_{-\infty}^{+\infty}|h(\tau)|d\tau<\infty$$
How this condition can be used to solve $(1)$ in time-domain? Can we get initial conditions using this condition?

Solving the characteristic equation
$$s^2+\frac32 s-1=0\tag{1}$$
gives the following homogeneous solution:
$$y(t)=c_1e^{-2t}+c_2e^{t/2}\tag{2}$$
The requirements of causality and stability are taken into account by choosing at least one of the constants $c_1$ and $c_2$ to be zero for either $t>0$ or $t<0$. If we're looking for a causal solution we know that $y(t)=0$ for $t<0$, i.e.,
$$y_c(t)=\big[c_1e^{-2t}+c_2e^{t/2}\big]u(t)\tag{3}$$
which can be interpreted as choosing $c_1=c_2=0$ in $(2)$ for $t<0$.
For a stable solution, we require $c_2=0$ for $t>0$ because stability means that there can't be exponential growth. For the same reason we require $c_1=0$ for $t<0$:
$$y_s(t)=c_1e^{-2t}u(t)+c_2e^{t/2}u(-t)\tag{4}$$
Eqs $(3)$ and $(4)$ are the general forms of the solutions given the requirements of causality or stability, respectively.
In both cases, the constants $c_1$ and $c_2$ are determined by requiring that $y''+\frac32 y'-y$ equals a Dirac delta impulse with weight $1$.
Note that when calculating the derivatives of $(3)$ and $(4)$ you need to use the product rule:
$$\big[f(t)u(t)\big]'=f'(t)u(t)+f(t)u'(t)=f'(t)u(t)+f(t)\delta(t)\tag{5}$$
Solving such problems using the Laplace transform is indeed much less tedious.

As an example I'll show how to arrive at the values of $c_1$ and $c_2$ for the causal solution $(3)$. Define $f(t)=c_1e^{-2t}+c_2e^{t/2}$. With $y_c(t)=f(t)u(t)$, the derivatives of $y_c(t)$ are
$$\begin{align}y_c'(t)&=f'(t)u(t)+f(t)\delta(t)=f'(t)u(t)+f(0)\delta(t)\\y_c''(t)&=f''(t)u(t)+f'(t)\delta(t)+f(0)\delta'(t)=f''(t)u(t)+f'(0)\delta(t)+f(0)\delta'(t)\end{align}$$
We need to satisfy the differential equation $y_c''(t)+\frac32 y_c'(t)-y_c(t)=\delta(t)$. Consequently, the term $f(0)\delta'(t)$ must vanish:
$$f(0)=c_1+c_2=0$$
Furthermore, the coefficients associated with the Dirac delta impulses must add up to $1$:
$$f'(0)+\frac32 f(0)=1$$
This results in the requirement $$-2c_1+\frac12 c_2+\frac32 (c_1+c_2)=1$$
which leads to $c_1=-\frac25$ and $c_2=\frac25$.
In exactly the same way you can derive the coefficients of the stable solution $(4)$, which results in $c_1=c_2=-\frac25$.


What is the meaning of Mutual Information beyond the numerical calculation? Beyond the raw equation for calculating mutual information, what does it mean in physical terms?
For example:
From Information Theory, we know that entropy is the smallest loss-less compression scheme that we can use on a alphabet with a specific probability distribution.

What would that mean in terms of Mutual Information?

Background:
I'm trying to calculate the mutual information of uni-gram words and determine which of two books did they come from.

essential $$ I(book; word)$$

Two additional intuitive takes on mutual information:

When two random variables are independent the joint distribution $p(x, y)$ and the product of the marginal distributions $p(x)$ and $p(y)$ are identical. One could thus assess the degree of independent between two random variables by computing a probabilistic distance between $p(x) \times p(y)$ and $p(x, y)$ - this distance being 0 when the two variables are independent. A common probabilistic distance between variable is the Kullback-Leibler divergence. If you take the Kullback-Leibler divergence between the joint distribution and the product of the marginals of two random variables, you end up with... mutual information.
From a compression / coding perspective, imagine you are given a sequence of $N$ pairs of observations $(x, y)$. You want to compress them into a file. Two strategies: storing all the (x) in one compressed file then independently all the (y) in another compressed file ; vs compressing the pairs. Using an optimal coder, the file size in the first case is $N \times H(X) + N \times H(Y)$, while in the second case the file size is $N \times H(X, Y)$. The second approach is more efficient if there's a relationship between the two observed variables! How many bits have we saved per observation? $\frac{N \times H(X) + N \times H(Y) - N \times H(X, Y)}{N} = I(X,Y)$! So mutual information tells us how many bit per observation do we save by coding two data streams jointly rather than independently.

I am not sure about your example, though... Mutual information is computed between two random variables (distributions). I can see how "book" can represent the distribution of words in a book ; but I am not sure what "word" means here. Mutual information also requires "paired" observations to be computed.


How does Adobe After Effects generate its "audio spectrum" effect? I'm trying to replicate the "audio spectrum" effect from Adobe After Effects. An example can be seen in this video:

Obviously, it has to be some variant of a fourier transform, but I've tried taking the fft of an audio sample in realtime and my implementation looks more like a big block of bars:

I imagine some kind of intermediate function is being applied to the fft, but I can't really figure out what it might be.

I've had a try at reproducing the effect, and I think these are some of the key elements:

You need a high-resolution FFT (large size; windows overlapping in time so that this doesn't come at the expense of frame rate) and to discard all but the lowest-frequency bins. You can tell this from the video because there are only rarely any harmonics visible and then only a second harmonic. The music you're hearing contains many more harmonics, so the display must be showing only the lower frequencies, not all the way up to $f_s/2$.
You need to dynamically adjust the vertical offset so that only the significant peaks are shown, hiding the broad-spectrum sounds whether they are continuous 'noise' or sharp drumbeats. What I did was take the median of the current frame's data and subtract it (as well as a manually set additional term), but there may be better strategies.
If still needed once you've discarded the high-frequency components, apply a frequency-varying offset so that the levels are even across the spectrum (rather than decreasing with frequency as your picture does). Whether this is actually needed depends on the frequency content of your input signal — thus the frequency response of the microphone it was recorded with, and so on.
The nice rounded rather than spiky peaks should be apparent once you've zoomed in on the low frequencies, but they are also dependent on your choice of window function. Use a window function with a large main lobe.

Pseudocode:

fft_length = 32768
display_length = 1024
freq_varying_offset = fft_length * 0.000005
for frame from 1 to ...:
    let start = frame / video_frame_rate * audio_sample_rate
    let samples = audio_data.subarray(start, start + fft_length)

    let buf = log_fft(window(samples), fft_length).subarray(0, display_length)
    let offset = fixed_offset - median(buf)
    for i from 0 to buf.length:
        buf[i] = buf[i] - offset + i * freq_varying_offset

    // now draw the bars on the video frame using the contents of buf


Visualising OFDM subcarriers in the frequency domain Im struggling to understand why in some images the frequency domain shows clear bars with amplitude and phase and with 0 frequency variation (it is clear that the frequency of the waves does not change in the image on the left).

Yet in others it shows it sloping or curving upwards:

Why isn't this a bar when the frequency of the signal clearly doesn't vary?

Same thing in OFDM. The constellation points should map to a single 312.5KHz aligned frequency and modify their amplitude and phase, yet when you look at the frequency domain on the OFDM carriers, they're sinc shaped waves instead of being bars. How is this possible?

Note that this illustration is just that, an illustration: The spectrum on the right can't be explained by the time-domain signal excerpt on the left alone. What's shown on the right is the spectrum that's typical for a tone with some phase noise, no matter where that comes from. But since frequency is the time derivative of phase: no, if you get something like the right spectrum, you don't only have a tone of a single, constant, frequency.

  Why isn't this a bar when the frequency of the signal clearly doesn't vary?

Because only periodic signals have line spectra. Your OFDM signal isn't periodic (for longer than an OFDM symbol), so it can't have a line spectrum.

Since any information transmission requires change of some parameters of a wave to actually contain any information, no information-carrying signal can have (with sufficiently long observation) have a line spectrum. End of story!

You'll notice that when reading up on pulse shaping, for modulations¹ where there's no DC offset², the pulse shaping defines the spectral shape.

In OFDM, the pulse shape is a rectangle the length of the OFDM signal. Hence, the sinc-shaped spectrum. (and the trick about OFDM is that the zeros of the individual sincs happen to fall exactly on the maxima of the other sincs, so that they don't interfere, hence the "O" in OFDM.)

¹ no matter whether we're considering a single-carrier system or just one subcarrier of an OFDM system

² i.e. all practical digital modulations aside from a few power detection-based specialities.


Is there anything called uncorrelated filter? I understand to an extent various filter like low pass filter, high pass filer, Wiener  filter Kalman filter etc. I also understand some of this filter will decorrelate/uncorrelate the signal. The question is: is there any filter which does nothing but uncorrelated a signal? What I mean by decorrelate is the following. Consider $x$ is the signal and we apply a filter $w$ to $x$, after the operation we get say $y$ now, I want $y$ to have zero auto-correlation for say $n$ number of samples, preserving the edges as surely those would form correlated section. An obvious answer (as pointed out in comment) of 0 as $w$ is not considered interesting.

I do not know if you are referring to a Whitening filter (check this). A Whitening filter can be viewed as a tool converting the correlated sequence into a white sequence. The application of this can be found in Wiener filter and speech compression. For example, in speech compression, the speech signals are highly correlated, which means they are highly redundant. Whitening them by converting them into white sequences (No redundancy) can reduce the complexity to compress the speech.


What is the significance of the linear phase in the phase response of an M-point moving average filter? I have plotted the magnitude and the phase responses of an M-point moving average filter, the following are the plots when M = 10:

We can observe that corresponding to every lobe in the magnitude spectrum we have a line in the phase spectrum - a linear phase response is observed corresponding to every lobe. (See the below picture)

I tried to observe the magnitude and phase responses for different values of M and observed a similar pattern.
Is there any reason behind this pattern? Is there any significance for this?

The significance is a constant time delay for all frequencies. Time delay is the derivative of phase with respect to frequency, so given a linear phase as shown, the time delay is constant. Notice that the abrupt steps in phase actually only occur when the magnitude goes through zero so does not effect the nature of it being a "linear phase" filter.
An M-point moving average filter will have a time-delay (for all frequencies) of $M/2$, so if you increase the number of points, the phase slope versus frequency will increase.
Why do we like linear phase filters so much? Please see this post on that:
Why is a linear phase important?
Any finite impulse response (FIR) filter with symmetric or anti-symmetric coefficients (The OP's case with coefficients all 1 is symmetric) will be a linear phase filter.


wavelet analysis normalized to have unit energy I am analyzing some time series of sea surface water temperature and am using wavelet analysis to do so. I am computing the continuous wavelet transform and then removing some specific frequencies and then converting the series back from the frequency to the time domain. I am following the methods of torrence and compo 1998 (Practical guide to wavelet analysis) to perform this operation but am a bit confused as to why my results are not as expected. As an example, I will fist decompose the time series into the frequency domain and then reconstruct the signal without altering any of the coefficients:

In the attached figure we have the original series in red and the reconstructed series in blue, straight away we can see that the the reconstructed signal is much lower than the original series. The main reason I think this occurs is that the wavelet transform normalizes the series to have unit energy (in order to compare the coefficients across scales). Therefore, I was wondering if it was possible to convert the reconstructed signal back to have the same magnitude (?) as the original series i.e. reverse the normalizing to have unit energy stage?

Take the mean of the original signal and add it to the reconstructed signal. IIRC wavelets have no DC response, and so the reconstructed signal will have 0 mean.


Compressive Sensing: What Class of Signals Are Exactly Model Sparse? In traditional compressed sensing, one would require $ m = O ( K.logN)$ measurements to recover a signal of size $N$, and sparsity $K$.

In this paper  Richard Baraniuk and his co-authors show that it is actually possible to get rid of the $log N$ factor in the number of sampled measurements, provided the $K$-sparse signals follow some additional characteristics. In traditional CS, in general, a $K$-sparse signal ( or a $K$-compressible one) lives on (or near) one of the canonical $K$-dimensional subspace of the basis of measurement. In model sparse signals, the additional requirement is that only a subset of these subspaces are allowed. These subspaces are such that they obey some well-formed characteristic. For instance, a Tree-sparse signal is one for which the top wavelet coefficients are to form a connected and rooted subtree in a wavelet tree decomposition.

I have to admit that I am still trying to understand their paper; so it may be a bit naive to ask it, but here is the question.

What sort of signals can one represent with these tree-sparse forms?  The above paper gives a few sample signals to that effect. But is there a general characterization of these signals?
For example, consider a signal comprising of, say $O(N/d)$ gaussians, with peaks randomly located among the N points ( $d < N$ is a constant in this case).  Can this be represented/approximated by one of these Tree-sparse models?

Any general pointer to literature that studied these objects?

In the case of general Compressed Sensing, we assume that we have no prior information about the location of the sparse support of a given signal. If we knew the support exactly, that is, if an oracle told us where the non-zero coefficients are located, we would only need $K$ measurements to recover our signal accurately. Since we do not generally have such a genie, we require more measurements to both determine both the support and the signal values on that support, hence $M$ on the order of $O(K \operatorname{log}N)$.

The authors in the Model-Based Compressed Sensing technique make the observation that if an a priori assumption can be made about the location of signal support, fewer measurements are required (since we know more about the support, statistically, we need to deduce less information from our measurements).

When a signal obeys a tree-structured, or hierarchical, model for its support, we have an easy way of predicting the support. This knowledge allows us to lower the number of required measurements for accurate reconstruction. This kind of model is very common since many classes of interesting natural signals can be represented compactly using Wavelets. The multiresolution nature of wavelet analysis groups significant (i.e. non-zero) coefficients into a tree-structure: significant coefficients at one resolution scale are highly correlated with spatially co-located coefficients in neighboring resolution scales.

(from: http://www.dominikdeak.com/coding/wavelet-compression/figure_17.png)

In the example you give about a mixture of Gaussians, it is not immediately clear that Model-Based CS can be used, since, as you say, the parameters (the peaks) of these Gaussians are randomly distributed within the signal. This kind of signal can be sparsely-represented (pretty much ideally) with some redundant dictionary of Gaussians of different parameters.
The sparse support when represented within this dictionary would not necessarily obey
any kind of model, the measurements alone would be needed to determine the support (Gaussian parameters).

On the topic of Compressed Sensing recovery of Gaussian Mixture Models, I ran a quick google search and came up with the following on arXiv:

Yu and Sapiro, "Statistical Compressed Sensing of Gaussian Mixture Models."


calculating the stop-band edge frequency of a low-pass FIR I have three information about a low-pass FIR filter.

pass-band edge frequency
pass-band ripple
stop-band ripple

I want to get the stop-band edge frequency I've been googling for the last hour to find a formula to get, but without success .
can any tell me how can I calculate it if have those 3 values ?
thanks in advance

The question is what you really want to do. Because the answer to your question is that there is no way to exactly compute the stopband edge from the three values you have. Looking at your specs I guess you're talking about a Parks-McClellan equiripple FIR filter, for which there are heuristic formulas relating the width of the transition band and the maximum ripple to the filter order (see this answer).

Usually in a design problem the passband and stopband edges would be given, and the maximally allowable error (ripple) would be given also, and you would need to figure out (by trial and error and/or by some heuristic formulas) what the required filter order would need to be. The problem you asked about doesn't occur in practice.

As a last comment, if you were given the actual filter coefficients then, by definition, the stopband edge is of course the lowest frequency for which the filter's magnitude response equals the given stopband ripple value (for a low pass filter, that is).


Why this two FIR filters have the same amplitude but different phase response Why a filter of the form

$$
1+a_0 z^{-1} +a_1 z^{-2}+ a_2 z^{-3}
$$
has the same amplitude response when in the reversed form
$$
a_2+a_1 z^{-1} +a_0 z^{-2}+z^{-3}
$$
but the phase response is different. I don't get it, and can't understand what the differences in terms of the relative group delays mean.

If you have a causal length $N$ FIR filter with impulse response $h[n]$, and with frequency response $H(e^{j\omega})$, and you invert it on the time axis, you get a new filter

$$g[n]=h[-n]\tag{1}$$

This filter is non-causal, but since it is an FIR filter, it can be made causal by shifting it to the left by $N-1$ samples:

$$f[n]=g[n-(N-1)]=h[N-1-n]\tag{2}$$

This is what you're doing. Now since

$$H(e^{j\omega})=\sum_{n=-\infty}^{\infty}h[n]e^{-jn\omega}\tag{3}$$

you have

$$G(e^{j\omega})=\sum_{n=-\infty}^{\infty}h[-n]e^{-jn\omega}=\sum_{n=-\infty}^{\infty}h[n]e^{jn\omega}=H(e^{-j\omega})$$

and

$$F(e^{j\omega})=G(e^{j\omega})e^{-j(N-1)\omega}=H(e^{-j\omega})e^{-j(N-1)\omega}$$

From (3) you can see that for real-valued $h[n]$ you have $H(e^{-j\omega})=H^*(e^{j\omega})$ (where $*$ denotes complex conjugation). So you finally get

$$F(e^{j\omega})=H^*(e^{j\omega})e^{-j(N-1)\omega}\tag{4}$$

and since $|e^{-j(N-1)\omega}|=1$, and $|H^*(e^{j\omega})|=|H(e^{j\omega})|$ you have

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

So the magnitude responses of the original filter and the time-reversed filter are identical. The phase responses are different: from (4) you have

$$|F(e^{j\omega})|e^{j\phi_F(\omega)}=|H(e^{j\omega})|e^{-j\phi_H(\omega)}e^{-j(N-1)\omega}\tag{6}$$

from which

$$\phi_F(\omega)=-\phi_H(\omega)-(N-1)\omega\tag{7}$$

follows. Since the group delay is the negative derivative of the phase, the group delays are related by

$$\tau_F(\omega)=-\tau_H(\omega)+N-1\tag{8}$$


Phase change rate of human speech I have a .wav file of human speech, sampled at 48kHz, with useful signal up to 5kHz, and almost no noise. If I do Hilbert transform in MatLab, I get the Analytic signal:

$$ s_\mathrm{a}[n] = s_\mathrm{real}[n] + j \cdot s_\mathrm{hi}[n], $$

where $s_\mathrm{hi}[n]$ is a Hilbert transformation of $s_\mathrm{real}[n]$ from the .wav file. Now for each sample I can calculate the phase:

$$ \phi[n] = \operatorname{Arctan}\left( \frac{s_\mathrm{hi}[n]}{s_\mathrm{real}[n]} \right) $$

and take time derivative from it. That will be the speed of phase changing. I don't want to call it "instantaneous frequency", that's too controversial.

Is any smarter way of doing this in MatLab? Probably this:

Fs/(2*pi)*diff(unwrap(angle(hilbert(s_real))));

I did this in Adobe Audition plugin, in Delphi. Way too many spikes, every time when the envelope $\sqrt{s_\mathrm{real}^2[n] + s_\mathrm{hi}^2[n]} $ was around zero.

I will not address the question whether or not it is meaningful to compute the instantaneous frequency of a speech signal. Instead I will show you a better method for computing the instantaneous frequency from a given analytic signal. This method avoids the phase unwrapping problem by directly computing the instantaneous frequency from the real and imaginary parts of the (discrete-time) analytic signal

$$s[n]=s_R[n]+js_I[n]\tag{1}$$

where

$$ s_I[n] = \mathscr{H} \big\{ s_R[n] \big\} $$

with $\mathscr{H}\{\cdot\}$ denoting the (discrete) Hilbert transform.

First, note that the formula for computing the phase given in your question only works if the real part of the analytic signal is greater than zero. In general you have to use the atan2 function implemented in most programming languages. In Matlab you can of course use angle, as shown in your example code.

Second, the phase becomes irrelevant and undefined if the magnitude of the respective complex number is (close to) zero. So it's pointless trying to compute the phase angle of a complex number with a negligible magnitude. This is true in general, so it doesn't matter how smart your method for computing the instantaneous frequency is.

Let's define the instantaneous frequency of a discrete-time analytic signal as the forward difference

$$f[n]=\frac{\phi[n+1]-\phi[n]}{2\pi T}\tag{2}$$

where $\phi[n]$ is the phase angle of the analytic signal $s[n]$, and $T$ is the sampling period. Note that other definitions such as the backward difference or central difference are also possible and useful.

We will not use $(2)$ directly for computing $f[n]$ because it requires phase unwrapping, which is especially problematic in the presence of noise. With

$$s[n]= \Big|s[n]\Big|e^{j\phi[n]}\tag{3}$$

we get

$$s[n+1]s^*[n]=\Big|s[n+1]\Big|e^{j\phi[n+1]}\cdot \Big|s[n]\Big|e^{-j\phi[n]}=\Big|s[n+1]\Big| \Big|s[n]\Big|e^{j(\phi[n+1]-\phi[n])}\tag{4}$$

(where $^*$ denotes complex conjugation), and, consequently,

$$f[n]=\frac{1}{2\pi T}\arg\Big\{s[n+1]s^*[n]\Big\}\tag{5}$$

The argument should be computed using the atan2 function. The result of the argument computation is only defined for $s[n]\neq 0$ and $s[n+1]\neq 0$.

Note that using backward or central differences in the definition of the instantaneous frequency would have resulted in the following expressions:

$$\begin{align}f_b[n]&=\frac{1}{2\pi T}\arg\Big\{s[n]s^*[n-1]\Big\}\\f_c[n]&=\frac{1}{4\pi T}\arg\Big\{s[n+1]s^*[n-1]\Big\}=\frac12\big(f[n]+f_b[n]\big)\end{align}$$

Expressing Eq. $(5)$ in terms of the real and imaginary parts of the analytic signal gives:

$\begin{align}
f[n] &=\frac{1}{2\pi T}\arg\Big\{\big(s_R[n+1]+js_I[n+1]\big)\big(s_R[n]-js_I[n]\big) \Big\}\\
\\
&=\frac{1}{2\pi T}\arg\Big\{s_R[n]s_R[n+1]+s_I[n]s_I[n+1]+j\big(s_R[n]s_I[n+1]-s_R[n+1]s_I[n]\big)\Big\}\\
\end{align}\tag{6}$

More details about instantaneous frequency (and instantaneous bandwidth) can be found in the paper "The calculation of instantaneous frequency and instantaneous bandwidth" by A.E. Barnes.

Also take a look at this related answer.


Different representations of frequency space of 2D image FFT I'm learning images processing using FFT.  In my test example provided below the input pixel values are clamped 0-1 (0-255), but I do eventually want to process floating point heightfield pixel values.

The software I'm using (Houdini v18) provides forward and inverse FFT functions and as a base line I can successfully convert an image to frequency space and back.

However when I look at the 2D representation of the frequency space it looks different to any examples I've found online.

Test image:

This is the result of Houdini's FFT with zero frequencies at center (height offset is pixel intensity):

The function returns 2x "images" representing real an imaginary values.

This however is the frequency space representation I find everywhere online:

From what I understand the radial type FFT images represent frequency between 0-2PI in u and v and the pixel value is the magnitude.  And that most images are offset so 0 is centered.

Do I need to convert the real and imaginary components to frequency and magnitude and plot those?  If so how?

EDIT:  My end goal is to apply radial pass filters to the FFT so I just need to apply any spacial transforms to get the FFT to that state.

The second image you find online as you mentioned is magnitude and the first image must be either real part or the imaginary part (but not both). If you calculate complex value magnitude for each point, i.e. $Mag_{i,j} = R_{i,j}^2 + I_{i,j}^2$ you get something similar to the following image.

I = imread('cameraman.tif');
F = fft2(I);
F_Centered = fftshift(F);
Mag = abs(F);
imshow(Mag);

but taking log10 harnesses the large magnititde of some frequency bins so for visulasation the normalization does not round them to zero.

imshow(log10(Mag),[])

Thus applying the filter on the above image which is only processed for visualization is not correct and inverse FFT of the values would not have much meaning.

To apply your filter you get the fft of the image, F(image) , then you need to make a same size mask image (all 1 and the frequency section you need to be removed must be zero, also remember symmetry feature of Fourier) and then apply do a element wise multiplication of fft values with this mask.

Also, I would suggest to watch this short lecture by Hoff on this subject:

https://youtu.be/02c6ohQV2TA


How does SC-FDE modifies $E_b$ (in AWGN function) In lot of matlab scripts of SC-FDE simulations I see $E_b = N$ (being $N$ the length of the block to be transmitted through the channel and Eb energy per bit) to generate the AWGN. I would like to understand why is that, and more specifically:
Assuming, instead of a block length $N$, I have 2 sequences length $N/2$. I FFT() each one, concatenate the frequency-domain sequences and send the result (of length $N$) through the channel. What would my $E_b$ in AWGN function be?

Here $E_b=N$ means $N$ pulses (i.e. PHY waveforms) of unit energy represent one information bit.
I suppose it is due to the S/P (serial to parallel) block that translates each bit of information to a vector of size $N$. So it is implicitly assumed that the energy of each component of the vector is one.

Regarding your second question, you should consider the number of information bits. If you have one bit of information to generate the two vectors of size $N/2$, then again $E_b=\frac{N/2+N/2}{1}=N$.


Optimising an input signal for a system to minimise distortion of output signal I have a discrete time-series signal which I am able to process before passing it through a system which distorts said signal. This system cannot be altered and is non-linear, but a signal can be passed through as many times as is needed (for iterative optimisation algorithms, for example). I would like to optimise the signal so that the output of the system is as similar to the original input signal as possible.

I'm currently unclear on exactly what problem I'm trying to solve is and would appreciate advice on potential approaches? I understand the concept of system identification, which gives me insight into what the system is doing, but am unsure how to apply this to the signal to try and counter these effects.

System is non linear so all bets are off.

Jokes aside though, wouldn't you just characterize the system, then generate an inverse system that the signal passes through first such that the effects of the original system are cancelled out?

Edit: another potential solution is to 'linearize' the system by finding the characteristics of it, then scaling your input signal such that the system behaves in a linear way. Then follow the inverse idea that I gave in my original comment.


How to calculate the power of a finite length signal? I am confused with these concepts. If the signal is expressed as r(t), I know the power of the signal is given by:

But if the the length of signal T is finite and cannot approach infinity, how can I calculate the power? Can I say the power is approximately equal to the following equation?

If the signal $r(t)$ is of finite length (and if integral of its square during its definition interval is also finite) then its average power given by
$$ \bar{P_x} = \lim_{T \to \infty} \frac{1}{T} \int_{-T/2}^{T/2} |r(t)|^2 dt = \lim_{T \to \infty} \frac{1}{T} \int_{0}^{L} |r(t)|^2 dt = \lim_{T \to \infty} \frac{[\text{constant}]}{T} = 0$$ will be zero...

Btw, you may compute an instantaneous power of the signal, simply as the square of it $|r(t)|^2$...


Object detection using several images how do people combine different information (e.g., RGB and depth image or RGB and thermal image) to make their object detection algorithms better?

If the single modality images A and B are perfectly registered, I can think of:

Write one object detector working on A, write one working on B and then write one regressor that takes the outputs from A and B and learns a final output.
Instead of working with feature vector f to describe a certain region (e.g. a bounding box) in the image A one could learn a prediction model based on a feature vector F = [f;g] where g is the corresponding feature for image B.

These approaches do not make as much sense if the images are not nicely registered, but there is a certain amount of misalignment.

So again: how do people combine different modalities for object detection and in particular deal with misregistration?

Usually, people tend to go with the second option, since its a joint approach. In either case, the methods fusing these two data modalities (RGB+Depth etc) together are termed multi-modal methods.

Also, another 3rd option is not to fuse on the feature level but to jointly minimize an energy function or jointly voting for the correct hypothesis.

One successful algorithm which works on RGB+Depth is LineMOD. The publications are here and here. (Modified Line). LineMOD takes the second path, and fuses the depth and RGB features. It also has an OpenCV implementation (a basic one).


Discrete Fourier transform When you perform a discrete Fourier transform (DFT), and compute the DFT coefficients, $c_k$ and $X(K)$, what do these values tell us about the signal it came from?

Wikipedia would be a better place to get a full answer, but in short: Discrete Fourier Transform is defined over a finite set of $N$ samples $x_n$ by the formula
$$
X_k=\sum\limits_{n=0}^{N-1} x_n e^\frac{-i2\pi kn}{N}
$$
The resulting DFT coefficients contains information about the frequency content of the original signal.


Where to get transtion matrix for Kalman filter? I am trying to understand the Kalman filter, and I am struggling to find out how to choose the transition matrix (denoted as $\mathbf F(k)$ or $\mathbf A(k)$. This matrix is used to update the weight error covariance matrix as follows

$$\mathbf P(k+1) = \mathbf F(k+1) \mathbf P(k) \mathbf F^\mathrm T (k) + \mathbf Q$$

In the paper On the Intrinsic Relationship Between the Least Mean Square and Kalman Filters it is description of "The Kalman filter for deterministic states.", what does not use the part with prediction, so it is without the transtion matrix.

But there is also the general Kalman filter, with the prediction according to transtion matrix. I have not understand how this matrix should be obtained.

I quess, that this transition matrix is from state space formulation of the model. But what should I do if I do not have a model? Is this mean that the general Kalman filter is usable only If I have exact model of the system?

Simple: if you do not have a model, you cannot apply the Kalman filter. Or you could and make up a model, but you cannot expect any of the optimality properties of the filter to hold. Based on that paper, it seems they assume that $\mathbf{x}_k$ is known or at least directly measureable...


Is this waterfall characteristic of QAM or OFDM (or something other) modulation? I recorded this signal at 2.45GHz, it is coming from a camera.

Is this QAM or OFDM ?

Based on those side-lobes I would guess OFDM but I am not sure.

Also there seems to be an interesting preamble at the top.

The sharp rectangular shape of the spectrum, with the immensely flat top definitely says "OFDM, with suitable whitening/PAPR reduction"; what constellation is used on the individual carriers is a bit hard to tell.

Can you, with a no-samples-left-out spectrogram (e.g. gr-fosphor sink, or by playing back your recorded file slower and using a suitably high update rate and suitable low number of FFT points in your QT waterfall sink, determine how long the OFDM "burst" is?

Take a bunch of samples from the end, and a bunch of samples from the start (like, one fourth of the burst duration), and calculate a cross-correlation. The maximum correlation should give you the time when what is used as the cyclic prefix at the end of the burst starts; that way, you'll be able to make an educated guess at the symbol number that makes up the CP, and by that the OFDM symbol length/IFFT size. Take frame, and try to manually resample it, so that both CP starts are exactly that (guessed) IFFT size apart. FFT the hell out of the frame; because your resampling had the purpose of time-stretching your signal so that the length of it actually (hopefully) matches the transmitter's sampling, you'd have a timing-synced demodulator.

Look at the spectrum and the constellation plot of the resulting signal; maybe you can already guess the constellation used to modulate the OFDM subcarriers from that.


IQ modulation principle I am new in the signal processing. I hope I can get help. I want to ask a question regarding IQ modulation work .

Why can’t I simply multiply a complex base-band signal with a cosine to modulate it?

If you have a complex baseband signal and you modulate it with a cosine then you shift the signal in frequency but your passband signal is still complex-valued, so it cannot be transmitted over a single channel (you would need two channels, one for the real part and one for the imaginary part).

You can, however, transmit the same information over a single real-valued channel by modulating with a complex carrier and taking the real part of the result.

Let $x(t)$ be the complex baseband signal, and let $\omega_c$ be the desired carrier frequency (in radians). Your suggestion results in a passband signal

$$s_1(t)=x(t)\cos(\omega_ct)=\frac12 x(t)e^{j\omega_ct}+\frac12 x(t)e^{-j\omega_ct}\tag{1}$$

If $x(t)$ is complex-valued, then $s_1(t)$ is also complex-valued, and if $x(t)$ has bandwidth $B$, then for transmitting $s_1(t)$ you need two real-valued channels (for the real and imaginary part of $s_1(t)$) with bandwidth $2B$.

However, if you define a passband signal

$$s_2(t)=\text{Re}\left\{x(t)e^{j\omega_ct}\right\}=x_R(t)\cos(\omega_ct)-x_I(t)\sin(\omega_ct)\tag{2}$$

where $x_R(t)$ and $x_I(t)$ are the real and imaginary parts of $x(t)$, respectively, then you can transmit $s_2(t)$ over a single real-valued channel with bandwidth $2B$.


Is it possible that adding noise to a signal can increase detection possibility? I am simulating a feature detector in noisy environment which consist of AWGN Noise and Impulsive Noise. But I am getting a strange result as Signal with Impulsive+AWGN noise has better detection probability then with Signal with only AWGN. I know its wrong at some point. How Is it possible that signal with more noise like Impulsive noise has better detection possibility. Please share your experience.

Our ears have a hearing threshold- i.e. a minimum amount of sound power that is required for them to detect.  Thus, a pure tone that is less than the threshold will not be detected.  If white noise is added to the tone such that the total power of the tone and noise exceeds the threshold, then the ear can hear the tone.

Perhaps you have some sort of similar non-linear effect going on.


USB microphone & interface selection As a part of my thesis work -in which I am using both acoustic and visual information- I am looking forward to implement a real time sound source localization sub-system. From the geometric principles, I know that at least 4 non-coplanar microphones are needed for 3D localization.

However, as I understand, using only a laptop to interface 4 microphones is impossible. I read that one sound card can only support 2. I also search through 4 microphone usb interfaces and found out quiet expensive solutions. I have also encountered to PS Eye, I'm not sure its sensory information is accessible through MATLAB or CLAM/C++ though.

As I only need an approximate information, I then decided to use only the azimuth angle information (planar confidence) of the source. For that, I know 2 microphones are enough. But I read about another theoretical problems for this case, like the need of calibration, harder sound seperation etc.

I have to ask is there a cheap interfacing solution, using 3 or 4 USB microphones? If not, would it be feasible to build up something from scratch and use amplified electret mics? I am really confused about hardware selection. Neither my advisor, nor the academic papers have too much practical information.

Thanks in advance.

Just get a multi-channel recording interface and some regular analog microphones. For example http://www.sweetwater.com/store/detail/Scarlett18i20/ has plenty of channels and microphone pre-amps.


Minimizing profile curvature of 2D image I have 2D grayscale image representing terrain heights. I am in search
for (fast) algorithm to minimize profile curvature of this image on selected scales (in scale-space sense).

Profile curvature is function of partial derivatives:

$$P_c = - \frac{(dx^2 \times dxx + 2 \times dx \times dy \times dxy + dy^2 \times dyy)}{((dx^2 + dy^2) \times (dx^2 + dy^2 + 1)^{1.5})}$$

and its geometrical meaning is nicely explained by this image:

I came up with completely ad-hoc iterative approach to minimize curvature
at specified scale, but I am little bit lost when trying to apply this accross
all scales (I am trying to utilize Gaussian / Laplacian image pyramids). It works by adjusting heights by adding height offset according to profile curvature at given point (height offset is negative in convex areas and positive in concave areas).

Results for single scale are not completely hopeless though:

Original data:

Processed data:

Result of this procedure should be formation of sharp V-shaped ridges
and valleys.

Is there any mathematically sound approach how to solve this kind of problems ?
Any tips are appreciated !

I am in search for (fast) algorithm to minimize profile curvature of this image on selected scales (in scale-space sense).

  ...
  ... to minimize curvature at specified scale, [...], trying to apply this accross all scales (I am trying to utilize Gaussian / Laplacian image pyramids). It works by adjusting heights by adding height offset according to profile curvature at given point (height offset is negative in convex areas and positive in concave areas).

Modifications of curvature are basically modifications of a "signal's" dynamics.

In effect, what you are trying to say is "If the signal is rising (or falling) faster (or slower) than some threshold then...do something different". This "do something different" can be anything from "don't respond" to "respond slower" or "respond faster" and so on.

In one dimension, the Douglas-Peucker algorithm does more or less what you are after.

The algorithm takes only one parameter which you can relate to the scale of observation but essentially modifies the local slope so that it remains within certain limits.

Variants that do surface simplification, or mesh simplification exist too (or this, in the case of contours) which use some expression of "cost" or "stress" to decide which element to remove during the simplification.

If you are "tied" to a height map (i.e. implicit 3D data), then you are looking at removing data (i.e. decimation) and interpolating / filtering. For example, you can use your $P_c$ to "mark" the areas of high or low curving slope in one pass and compute the "fill-in" values of each patch on a second pass. For examples of these approaches, you might want to take a look at this or this or more generally this.

Hope this helps.