What is the difference between FIR1 and FILTFILt in MATLAB? In MATLAB there are filters with 0 phase filtfilt and ordinary FIR1. When I filtered a group of simulated data, I found that filtFilt was superior to FIR1 in terms of phase delay. I know the difference. But would there be any essential difference between filtfilt-designed filter and FIR1-designed filter when written into FPGA? Or are they just different in Matlab when it is used to filter data?  Besides, could you please tell me how to set the order of filtfilt? I would like to compare the filtering performance difference between FILtFilt and FIR1 of the same order. Thanks!

filtfilt() applies forwards, then backwards filtering using any iir filter you feed it. As such, it only applies to finite length file processing, not realtime stream processing.
fir1 is a method to design fir filters that can subsequently be applied using filter, conv or (unusually) filtfilt.


Power contained in a random process $X(t)$ How do we calculate the AC and DC power of random process $X(t)$ , provided we have $R_x (\tau)$, and $S_x(f)$ ?

Assuming that the process is wide-sense-stationary, the total power in a random process with autocorrelation function $R(\tau)$ and power spectral density $S(f) = \mathcal F\{R(\tau)\}$ is given by
$$\text{Total Power} = \int_{-\infty}^\infty S(f) \,\mathrm df = R(0)$$ while the "DC power" (also equal to the squared mean of the process)  is given by
$$\text{DC Power} = \lim_{\tau \to \infty} R(\tau).$$ The DC power is also given by the magnitude or amplitude of the impulse (if any) at the origin in $S(f)$. Note that $S(0)> 0$ is not enough to claim that the  process has DC power -- for example, white noise has no DC power --; we need to have an impulse (a.k.a Dirac delta) at the origin in $S(f)$ for the process to have DC power.

What the AC power might be is left as an exercise for the reader to figure out.


How to Calculate Gaussian Kernel for a Small Support Size? How to calculate the values of Gaussian kernel? I think I understand the principle of it weighting the center pixel as the means, and those around it according to the $\sigma$ but what would each value be if we  should manually calculate a $3\times 3$ kernel?

Gaussian Kernel is made by using the Normal Distribution for weighing the surrounding pixel in the process of Convolution.

Since we're dealing with discrete signals and we are limited to finite length of the Gaussian Kernel usually it is created by discretization of the Normal Distribution and truncation.

I created a project in GitHub - Fast Gaussian Blur.
It uses many methods to approximate the Gaussian Blur Filter and evaluate their speed and quality.

The most classic method as I described above is the FIR Truncated Filter:

https://github.com/RoyiAvital/FastGuassianBlur/blob/master/ApplyGaussianBlur.m

Regarding small sizes, well a thumb rule is that the radius of the kernel will be at least 3 times the STD of Kernel.
If you chose $ 3 \times 3 $ kernel it means the radius is $ 1 $ which means it makes sense for STD of $ \frac{1}{3} $ and below.

Then just fill your Kernel By:

gaussianBlurRadius  = ceil(stdToRadiusFactor * gaussianKernelStd);
vGaussianKernel = exp(-([-gaussianBlurRadius:gaussianBlurRadius] .^ 2) / (2 * gaussianKernelStd * gaussianKernelStd));
vGaussianKernel = vGaussianKernel / sum(vGaussianKernel);


How to find out the transfer function of a FIR filter? $$h[n]=\begin{cases}a^n & \text{if } 0 \le n < N \\
0 & \text{otherwise}\end{cases}$$
And for which values of $a$ the filter is stable

I know that the transfer function will be

$$H(z)=\frac{z}{z-a}~,~|z|>a$$

how to find out the values of $a$ for it's stability ?


What is the steady state response of a system of two exponentials to sinusoidal excitation? What is the two exponential steady state response?

Here is an example solution for sinusoidal excitation of a system having a single exponential response to a impulse excitation:

Example:

Excitation: impulse: $L(t) = \delta(t)$

Impulse response of system: $I(t) = I_o e^{-t/\lambda}$

Excitation:  sinusoid: $L(t) = a + b \sin(\omega t)$

Sinusoid response of system:  $I(t) = I_o \lambda \left(a + \frac{b}{\sqrt{1 + (\omega \lambda)^2}} \cos(\omega t - θ)\right)$

$\tan(\theta) = \omega \lambda$

Question:  What is a similar solution for a system that has a double exponential impulse response?

Excitation: impulse: $L(t) = \delta(t)$

Impulse response of system: $I(t) = I_{o1} e^{-t/\lambda_1} + I_{o2} e^{-t/\lambda_2}$

Excitation:  sinusoid: $L(t) = a + b \sin(\omega t)$

Sinusoid response of system: ????

Thanks!!!!

The system with impulse response:
$$I(t)=I_1e^{-t/\lambda_1}+I_2e^{-t/\lambda_2}$$

with Laplace Transform:
$$I(s)=I_1\frac{1}{s+1/\lambda_1}+I_2\frac{1}{s+1/\lambda_2}$$

and the input:
$$L(t)=a+bsin(\omega)$$

with Laplace Transform:
$$L(s)=\frac{a}{s}+\frac{bw}{s^2+w^2}$$

will obtain the following output, in Laplace Domain:
$$IL(s)=I(s)L(s)=(I_1\frac{1}{s+1/\lambda_1}+I_2\frac{1}{s+1/\lambda_2}) (\frac{a}{s}+\frac{bw}{s^2+w^2})$$

with this very simple Inverse Laplace Transform, as the final solution:
$$IL(t)=
-\frac{I_1 \lambda_1 e^{-t/\lambda_1} (a \lambda_1^2 w^2 + a - b \lambda_1 w)}{\lambda_1^2 w^2 + 1}
- \frac{I_2 \lambda_2 e^{-t/\lambda_2} (a \lambda_2^2 w^2 + a - b \lambda_2 w)}{\lambda_2^2 w^2 + 1}
+ a I_1 \lambda_1 + a I_2 \lambda_2
+ \frac{\{b (sin(t w) (I_1 (\lambda_1 \lambda_2^2 w^2 + \lambda_1) + I_2 \lambda_2 (\lambda_1^2 w^2 + 1))
- w cos(t w) (I_1 \lambda_1^2 (\lambda_2^2 w^2 + 1) + I_2 \lambda_2^2 (\lambda_1^2 w^2 + 1)))\}}{(\lambda_1^2 w^2 + 1) (\lambda_2^2 w^2 + 1)}$$

As you can see, this kind of symbolisms can be readily solved in Wolfram Alpha through this:

inverselaplacetransform((c/(s+1/m)+d/(s+1/n))*(a/s+bw/(s^2+w^2)))


LTE system. How to compute precoder for resourse block guys!

I am studying a massive Mimo system and 3gpp LTE system and have questions about precoding in the LTE system. In LTE RB (resource block) is used to group consecutive subcarriers.  To compute precoder and reduce the complexity of computation it is possible to compute one precoder for RB.

In this post, the frequency selective channel was discussed and mentioned: "to average the second-order statistics."  Does it mean: covariance, covariance function? I have read that it is used for signal, but how to use for RB ... I dont know...

Does anyone give me more details?

Think of the RB as a subset of the OFDM channel set. That reduces the complexity of finding $N_\text{FFT}$ matrix decompositions of $N_\text{TX antennas}\cdot N_\text{RX antennas}$ matrices for each user to
$N_\text{carriers in the individual RB}$ decompositions for each user, which might well be smaller, because at user-level, often a significant amount of TX antennas can be ignored to begin with.


What Matlab Classifiers to use I'm currently trying to decide what classifier to use for my snore detection algorithm. Which uses some features like Zero cross rating, frequency band energy, first formants...

The classifiers which have been pointed to are:

Support Vector Machine
Naive Bayes Classifier
Knn Nearest Neighbour
Linear or Quadratic Classifier

However, could somebody describe what situation each classifier is used in,  or point me to a good paper, website for article which describes each classifier well.. I find the Mathworks website hard to understand the difference between each.

thanks

If you don't know what you're doing, naive bayes is hard to screw up. The naive part of naive bayes is that it assumes that all of your features are independent. They probably aren't independent but in many cases it doesn't matter so much that this assumption is violated. Naive bayes is also good at handling different feature types such as using both categorical variables and continuous variables for the same classification task.


Periodicity of sum of discrete signals In my lecture slides for school and from this website here

"The sum $z[n] = x[n] + y[n]$ of periodic signals $x[n]$ with fundamental period $N1$, and $y[n]$ with fundamental period $N2$ is periodic if the ratio of periods of the summands is rational, i.e.,
$\frac{N2}{N1}=\frac{p}{q}$
and $p$ and $q$ are integers not divisible by each other. If so, the fundamental period of $z[n]$ is $qN2 = pN1$."

I'm confused with a number of things. First if $y[n]$ and $x[n]$ were periodic, then wouldn't the ratio of their periods $\mathbf{always}$ be rational since each period is always integer? Also if each signal is not periodic then this would imply their sum could be $cos[n]+cos[n]=2cos[n]$ is not rational but the ratio of their "periods" is $1$ which is rational. So, I'm assuming that this requirement is just redundant if the original signals are periodic?

First if y[n] and x[n] were periodic, then wouldn't the ratio of their periods always be rational since each period is always integer?

No. I mean. Yes.

$x(t)$ might be e.g. $2\pi T$-periodic, and $y(t)$ could be $5T$-periodic ($T$ being the sample period). The choice of using capital $N$ for the period was a confusing one, but signals don't need to have integer sample-multiples as periods.

BUT, and that's where I'm with you, if you say $x[n]$ is periodic, it means the sampled signal $x[n]$ is periodic – which it can't be, if the original sample would have a rational-amount-of-samples period.

And thus, you're right: if the sampled signals are periodic, they need to have a rational ratio in periods, and thus, their sum will always also be periodic.


Result of conjugate symmetry property of DFT I know one of the properties of DFT for real-valued time series is conjugate symmetry.

But what does it imply? In the textbook it says that for a DFT of the length M, this makes M/2-1 spectral coefficients redundant. Could someone please explain this to me? How do we conclude this?

The textbook probably even has the proof; but it's not hard to do: Simply write down the sum that gives you the $k$th bin of the transforf $X[k]$ of the time series $x[n]$. Then, flip $k$'s sign and apply the conjugate to finish your proof.

To show:

\begin{align}
X[k] &\overset{\text{claim}}= X^*[-k]\\
X[k] &= \sum_{n=0}^{N-1}x[n]e^{-j2\pi k\frac nN}\\
X^*[-k] &= \left(\sum_{n=0}^{N-1}x[n]e^{-j2\pi (-k)\frac nN}\right)^*\\
&= \sum_{n=0}^{N-1}\left(x[n]e^{-j2\pi (-k)\frac nN}\right)^*&\text{Sum operator linear}\\
&= \sum_{n=0}^{N-1}x[n]\left(e^{-j2\pi (-k)\frac nN}\right)^*&x[n]=x^*[n]\text{ since real}\\
&= \sum_{n=0}^{N-1}x[n]e^{-(-j2\pi (-k)\frac nN)}\\
&= \sum_{n=0}^{N-1}x[n]e^{-j2\pi k\frac nN}\\
&=X[k] &\rule{1ex}{1ex}
\end{align}

Now, having showed that half of the coefficients are just the conjugate of the others: It's obviously sufficient to only consider one half; the other can be directly calculated from that.


How to determine sampling frequency for an x(t) signal avoiding aliasing? The antitransform of the function is given:

$ \hat{x}(f) = \frac{123 + i246\pi f}{246 - 24600 \pi^2 f^2 + i4920\pi f} $

I'm asked to determine which frequency can I sample using the function x(t) avoiding aliasing? (or which range of values)

Update the problem actually suggest five answers:
100, 10, 1, 0.4 or none of these. Turns out the answer is none of these.

The function you show has a frequency range of $-\infty$ to $\infty$ so it isn’t strictly band limited. To completely avoid aliasing, it needs infinite sampling.

Of course, you can use tables or Fourier identities to obtain the exact inverse.

Usually some acceptable criteria is specified for tolerable aliasing when specifying a sample rate.


Auto or cross correlation for synchronisation in MATLAB I have multiples records of  a motor sound, I would like to synchronize them first and then  calculate their spectrogram. My question is how can I use the autocorrelation or the crosscorrelation to synchronize them? I want to make sure that the only difference I see in the spectrogram is related to the frequency, not to the start position of the motor or the records.

I know to calculate the auto and cross correlation in MATLAB but I don't understand how to use them for synchronization purpose.

Thanks for any help!

Mostly it is a very simple task. Here is the example for shifted sinusoid:

t = linspace(0,6*pi,1000); % time vector for sinusoid
s = [zeros(1,100), sin(2*pi*t), zeros(1,100)];  % original signal
shift = 30; % we are shifting our sinusoid 30 bins to the right
ss = [zeros(1, shift) s(1:end-shift)]; % signal shifted
subplot(2,1,1)
plot(s); hold on;
plot(ss, 'r')

[c, lags] = xcorr(ss, s); % get the cross-correlation together with lags
subplot(2,1,2)
plot(lags, c)

pred_shift = lags(c==max(c)); % calculated the shift from maximum value of cross-correlation

The output is: 30 and you can use it to shift back the signal ss. In your case of engine signals it will be analogous.


Does the number of samples matter for FFT, and how to get a specific frequency visible? If I have a signal that has 740 samples, the sampling frequency is 1000 Hz, and I take an FFT of length 1024, are the frequencies I get now 0, fs/1024, 2*fs/1000, ..., 511*fs/1024? I mean, do I get this correctly: the frequency resolution does not depend on the number of samples I have, only FFT length and sampling frequency. If I had 1480 samples, the resolution would be still the same? This makes my head hurt.

Next, suppose that the signal I have has an interesting frequency at 26 Hz (strictly that one frequency). How do I get it appear properly in the frequency spectrum (my frequency bins are spaced 1000/1024 = 0.9766 Hz apart, so am I at dead end?).

The "information content" is related to the number of samples.  The length of the FFT merely interpolates the spectral frequency curve represented by the number of samples.  If an FFT result bin doesn't end up exactly centered on some desired frequency, you can interpolate it from the nearest bins, using a high-quality (Sinc kernel, et.al.) interpolator, since FFT bins have a greater-than-zero bandwidth (the shape of the transform of the default rectangular or other window).  Using a longer zero-padded FFT gives you more interpolation points, and thus a denser looking point plot, but no better spectral separation resolution than provided by the number of actual samples.


Adding white noise to complex signal (complex envelope) Good day everyone.
So, I run in to trouble - How to add white noise to signal samples?
For now I'm working with MSK signal - quadrature modulator and quadrature demodulator + Viterbi processor.  My signal samples at TX and RX side represents Q and I signals of modulator(real and complex part of the signal). Model works nice in ideal situation, without any noise and other interference.
So, I'm working with complex envelope of the MSK signal with amplitude equal to one. Currently, I would like to test my demodulator with some noise.
It looks like, I just need to add (+) samples from my noise generator to samples of transmitted signal. But I'm not sure if this is correct way,

Should I sum both Q and I channels to one sample of noise signal or not?
And I'm pretty confused how to determine signal to noise ratio after this. I would like to measure bit error ratio at some SNR points. As I read in some books changing $\sigma$ parameter of white noise generator will play a role here, but still how to determine SNR in numbers after this?

If you want a complex white noise sample with variance $N_0$, then you should generate two independent noise samples each with variance $\frac{N_0}{2}$. The two samples that you generate make up the real and imaginary parts of the desired complex noise sample. The resulting complex value will have the desired variance $N_0$. You should then sum this with your complex baseband signal to get the effect that you want.


Prove the following DFT for a periodic sinc sequence or ideal lowpass filter I need to prove the following DFT for a periodic  sequence or ideal lowpass filter:

$\sqrt{\frac{p}{m}}\frac{sinc(\frac{\pi np}{m})}{sinc(\frac{\pi p}{m})} \stackrel{DFT}{\longleftrightarrow} \begin{cases}
                    \sqrt{\frac{m}{p}},\: if \: |k-\frac{m}{2}|\ge \frac{p-1}{2}\\
                     0, \: otherwise
                 \end{cases}$

But actually, I do not know how to do it. Maybe someone gives me at least a hint?

Hints:

1) Definition of the unnormalized DFT:

$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-i\frac{2\pi}{N}kn}  $$

2) Geometric Summation Formula:

$$ \sum_{n=0}^{N-1} r^n = \frac{1-r^N}{1-r} $$

3) Exponential form of Sine:

$$ \sin( \theta ) = \frac{ e^{i\theta} - e^{-i\theta} }{ 2i } $$

4) Definition of unnormalized sinc:

$$ \operatorname{sinc}(x) = \frac{\sin(x)}{x} $$


Image Processing method for spiky selection I have an input as a 3D binary image and the preferred output below:

Input:

Preferred Output:

What image processing methods should I look for if I am to have only the spiky object(s) remain, just like the preferred output above?

There are more corners on the borders of your "spiky object", so one approach would be to tune a corner detector for this.

For example, I calculated the determinant of the structure tensor (Mathematica code below) of a distance-transformed image:

Binarizing with hysteresis yields this image, which should be a good starting point for the segmentation algorithm of your choice:

Mathematica code (src is the source image you posted)

At first, i calculate a distance transform of the input image. This creates contrasts over the whole object area (instead of just the border), so the whole object can be detected.

dist = ImageData[DistanceTransform[src]];

Next I prepare the components of the structure tensor. The filter size for the gaussian derivatives if 5, the window size is 20.

gx = GaussianFilter[dist, 5, {1, 0}];
gy = GaussianFilter[dist, 5, {0, 1}];
gx2 = GaussianFilter[gx^2, 20];
gxy = GaussianFilter[gx*gy, 20];
gy2 = GaussianFilter[gy^2, 20];

To calculate the corner filter at each pixel, I simply plug these into the symbolic determinant of the structure tensor:

corners = Det[{{dx2, dxy}, {dxy, dy2}}] /. {dx2 -> gx2, dxy -> gxy, dy2 -> gy2};

Which is basically the same as:

corners = gx2 * gy2 - gxy * gxy;

Converting this to an image and scaling it to 0..1 range yields the corner detector image above.

Finally, binarizing it with the right thresholds gives the final, binary image:

MorphologicalBinarize[Image[corners], {0.025, 0.1}]


How exactly do I plot a constellation Diagram? I am relatively new to signal processing and I would like to correctly draw a constellation diagram to check the results of a study. I need my signal processing to be as close as them to make sure I am doing it correctly.
I am trying to use constellation diagrams to fingerprint ADS-B signals which use a PPM (pulse position modulation) scheme. Considering I got some IQ data, it is sufficient to do a scatter plot of the IQ data to call it a constellation diagram ? Or do I need to collect only samples from IQ data that are in special positions(like the position of the pulses)? I am also asking this question because I do not get exactly the plot that the study shows for constellation diagram. I would be very thankful if someone can point why there is a difference or something I missed/misunderstood.
Here is the link of the study
You can find the plots on page 4, FIGURE 8. (What they call contour stellar image is just a constellation diagram but with colors based on density of points nearby from what I understood (FIGURE 7)).
And Using the following python code:
from scipy.stats import gaussian_kde
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)

xy = np.vstack([in_phase,quadrature])
z = gaussian_kde(xy)(xy) #This is just to get colors to see better

ax.set_xlabel("in_phase")
ax.set_ylabel("quadrature")

ax.scatter(in_phase, quadrature, c=z)

I get the following images (2 exemple):

Typically a constellation diagram is not just the raw RF values because that would be extremely disorganized. You would end up with a bunch of points in a gaussian distribution about the origin (thermal noise) and a "ring" around it that represents your signal at any particular phase during the sampling.
What you typically want to plot is the baseband signal, the RF after it is filtered, mixed against the carrier wave, and coherent (phase-locked). This will indicate which constellation values are being received with some noise (hence why in Figure 7 they are scattered about the data positions). I have no idea what is going on in Figure 8, but it's not a constellation diagram. It looks like something of the author's own invention to better analyze the given data.


Average spectrum representation - SonicVisualiser I try to have a signal (speech) represented in time and frequency domain  (waveform, spectrum and spectrogram). I use Sonic Visualiser and for spectrum representation it shows an analysis at a given point in time.
I wonder what representation should I use to represent a selected area of the waveform e.g. a word, in frequency domain (spectrums). Sonic Visualiser doesn't give averaging options.

Vasiliki

It would not be useful to average in time because vocalisations are represented directly by their spectrogram. (For example: http://www.u.arizona.edu/~ohalad/Phonetics/notes/Formants%20Spectrograms%20and%20Vowels.PDF or even http://www2.coloradocollege.edu/dept/ev/Research/Faculty/OVALItems/FordOsborneVocabulary/_SouthernVocabularyTable.html)

I would suggest that you export your spectrogram as an image and then isolate and layout the areas you are interested in, in some other software like GIMP. For more information, please see: http://www.sonicvisualiser.org/doc/reference/1.9/en/#imageexport and http://www.gimp.org/

Hope this helps.


What specifically does PCM refer to? What specifically is PCM? I've seen it refer to:

the sampling and quantisation of an analogue domain signal into digital samples

those sample quantisations alone in the digital domain i.e. in memory or storage i.e. .pcm files

the transmission of quantisations or indeed any bitstream as a digital signal with a digital line code through a DAC into the analogue domain

the process of 1+2

the process of 1+3

I like this paper a lot: The philosophy of PCM. It explains PCM very clearly, includes references to some of the earliest papers on it, and you can hardly argue with the qualifications of the authors :)
PCM was conceived as a three-step process. First, sampling (at frequency $f_s$ sufficiently above Nyquist). Second, quantizing (with enough resolution $\Delta$ to obtain the required quantization SNR). Third, coding, or mapping the quantized samples to $M=2^k$ pulse amplitudes, with $M$ chosen to obtain the required line pulse rate.
PCM, then, has parameters $f_s$, $\Delta$ and $M$. Each of these affect system performance in different ways.
Why use the term "modulation"? In communications, "modulation" is very general, and is used almost as a synonym to "mapping" or "converting elements from an input set to an output set". In fact, what you understand as "modulation" (multiplying by a carrier) is often called "upconversion" and "downconversion", to avoid confusion.


How to compute the impulse response from incomplete frequency response data? I need to compute the impulse response of a system, but I only have access to the frequency response data.

This data contains the output magnitude of the system (for a fixed input amplitude) at frequencies between 10Hz and 400Hz, sampled at each 5Hz (10, 15, 20, ..., 395, 400) with a 1000Hz sampling frequency. The magnitude past the area of interest tends towards zero, but was not measured directly.

I dont have phase information (or a way to obtain it), but I only care about magnitude information.

Here's a quick graph of the data I have:

(Amplitude units are undefined, but aren't meaningful since an arbitrary scaling factor is added during processing steps.)

How can I obtain an impulse response of the system with incomplete frequency response data?

I'm thinking of using the IFFT function in MATLAB, but I am unsure of how it will interact with the gaps and the missing information past 400Hz.

I assume I could pad the data with zeroes, but up to how many?

I can also obtain more data to "fill in" make the gaps smaller, but there will always be some since I could also sample at non-integer frequencies.

EDIT for extra context:
I am trying to implement a parallel filter design to equalize the magnitude response of my system, based on the works of Balasz Bank. The author provides MATLAB code to determine the parameters of the system, but said code requires the impulse response of the system as input, not the frequency response.

The code takes the impulse response as input, converts it to minimum-phase and then computes the filter parameters using a least squares method. (The actual code is in the roomcomp.m and parfiltid.m files)

The end goal of the equalizer is to get a flat frequency response for the amplitude (In this case, amplify 0-50Hz, attenuate 50-200Hz and amplify 200+Hz). The system takes in a 1kHz signal as input, equalizes it and outputs the result in real time at 1kHz.

I'd like to take phase into account eventually, but due to time constraints, for now any solution that results in decent equalization works.

Phase alone or magnitude alone reconstructions of the complete frequency response (and therefore the impulse response) is claimed in IEEE papers of Monson Hayes and Alan Oppenheim such as this https://www.researchgate.net/publication/3176816_Signal_Reconstruction_from_Phase_or_Magnitude)

They are specifically dealing with real systems whose Fourier transforms have Hermitian symmetry (redundancy), therefore half of the information (phase or magnitude) in principle should be enough to define the whole system. Yet the procedure is not very simple, iterative and requires some mature DSP level.

Check it out. It may help you.


what is silence in audio .wav file and how to compute decibels? I'm a beginner of audio signal processing, and trying to understand some basic things about audio signals.
I primary work with .wav files, so please allow me to use wav as my data to describe my questions.
Q1: Sample values in .wav file represent the waveform amplitude?
When we load data from .wav file, we will have the audio sample values (typically signed ints), and I searched on google, these values are the amplitude of the audio signal, right?
See example code
In [69]: from scipy.io import wavfile

In [70]: frame_rate, values = wavfile.read("data/test.wav")

In [71]: frame_rate
Out[71]: 48000

In [72]: values
Out[72]:
array([[  0,   0],
       [  0,   0],
       [  0,   0],
       ...,
       [385, 247],
       [412, 266],
       [359, 232]], dtype=int16)

Q2: What is silence?
If the answer to Q1 is YES, then given bit-16, I tried to create a numpy array with all values of 10000, expecting it should be loud sound. However, to my surprise it's just no sound at all. Why? I though value=0 is silence, but it's obviously not.
See example code
import numpy as np
from IPython.display import Audio

# 100k frames
x = np.arange(0, 100000)

# sine wave, max value=10000
sine = np.sin(x) * 10000
Audio(sine, rate=50000)  # i.e. 2sec audio, and I can hear the sound

# dummy, all values = 10000
dummy = np.ones(x.shape[0]) * 10000
Audio(dummy, rate=50000)  # same 2sec, but silence

Q3: How to compute decibels from the .wav amplitude values?
Given Q2, now I'm lost and don't know how decibels look like given amplitude values, I thought bigger amplitude values should mean higher decibels, but now it's apparently not the case.

Q1: Sample values in .wav file represent the waveform amplitude?

Yes.

Q2: What is silence?

That depends on how you define it. The signal that you generated actually has a lot of energy but it's all at 0 Hz. Your DAC cannot reproduce 0 Hz, your speaker cannot reproduce 0 Hz, air cannot transmit 0 Hz, and your ear can't hear 0 Hz. There is lots of energy but you can't get it to you ear or perceive it. If you want a loud signal, create a sine wave with a frequency 1000 Hz and an amplitude of 10000. That will be loud.

Q3: How to compute decibels from the .wav amplitude values?

That depends on what type of decibels and for what purpose. You can always calculate dBFS (dB re digital Full Scale) as
$$20 \cdot log(x/x_{max})$$
where $x$ is your sample and $x_{max}$ the clipping point of your AD or DA converter, i.e. 32768 for signed 16-bit integers.
Sound pressure levels, dBSPL, are much more complicated and require calibration and in many cases also a frequency weighting function.


Obtaining Correct (Least Squares Sense) Affine Transform Parameters Between Two Images I have two images that I want to compute the affine motion model parameters. The model that I use is $$x' = a_1x+a_2y+a_3$$ $$y' = a_4x+a_5y+a_6$$
To calculate those 6 parameters, I picked 6 points (overdetermined system) between two images, and calculated those parameters by using MATLAB (A = X\X_primes). However, when I apply the transformation, even the points that I picked (X) do not translate into the specified locations (X_prime). I inferred that my least square solution should be wrong, and I am looking for how to make it better, i.e. what kind of points should I choose on the image to obtain better parameters for affine motion model?

As a complement to the comments above, I am adding the slides for Image Registration by Leow Wee Kheng, where one can find:

  Given two images, how to register one with the other?

  Determine the corresponding points between the images. Manually mark corresponding points, or
  Detect and match features between views (see lecture on feature detection and matching).
  Determine the transformation between corresponding points. Assume that all pairs of corresponding points are related by the same
  transformation. Compute parameters of transformation given
  corresponding points.

Importance is laid on corresponding points. Indeed, even with an over-determined system, if the chosen pixels an weakly informative (on flat or noisy parts of the image), the matching might not be efficient. One can instead detect features (edges, corners, textural patches), and compute the registration based on them.

Sometimes, one does not have the same points, or the same number of features between images. I recently have tested a efficient algorithm (CDP, or Coherent Point Drift) for non-standard chromatographic images (made of peaks or blobs), that can be tuned from rigid to smooth deformations:

Point-Set Registration: Coherent Point Drift, 2010
Our application, for the record: BARCHAN: Blob Alignment for Robust CHromatographic ANalysis, 2017


Generating spectrum of a spectrum / identifying frequencies of ripples present in a spectrum I have a spectrum as shown. The x-axis is frequency in GHz. The spectrum shows ripples that we can visually quantify as ~50 MHz ripples. I am looking for a method to calculate the frequency of these ripples other than by visual inspection of thousands of spectra. Since the function is in frequency domain, taking FFT would get it back into time domain (with time reversal if I am correct). How can we get frequency of these ripples?

You are not doing anything "wrong" but sooner or later you will end up "re-inventing" the tool that is doing exactly that already. That tool is the Cepstrum.
You wish to estimate the frequency of ripples that show up in the frequency domain and you would like to do this via the Discrete Fourier Transform (DFT) and that is fine.
(Your assumption that if you take the DFT again it will send you back to the time domain is not entirely correct though)
Suppose your time domain signal is $x[n]$ and after the DFT it becomes $X[k]$.
Now, $X[k] \in \mathbb{C}$ but you don't care about the phase of the phase component (hope you understand why I wrote it this way). So, we take the magnitude of $X[k]$, let's call it $Q[k] = |X[k]|$. Notice here, $k$ still denotes frequency, capital letters still denote spectra.
$Q[k] \in \mathbb{R}$, so if you send it through the DFT, it does not care that this used to be complex or the result of another DFT or anything. It still sees a real "signal" and it decomposes it to a sum of sinusoids. If you $W = \mathcal{F}(Q)$, your $W$ will now be complex again, but this time it describes the fluctuations of $Q$.
Now, the thing is that because of the way the sums are put together within the DFT, the dynamic range of $Q$ will be all over the place. Specifically, the low frequencies (DC for example) can end up having HUGE sums and the high frequencies can end up having tiny little values. If you don't account for that, your $W$ will have a strong $\frac{1}{k}$ influence.
To account for this, instead of $W = \mathcal{F}(Q)$, you $W = \mathcal{F}(\log(Q))$. The effect of $\log()$ is to "compress" large values and boost small values. It is the same "trick" we do when we plot a spectrum in a logarithmic scale rather than linear scale.
But, if you do all this (which is a very reasonable thought progression), you have derived the cepstrum. Specifically, the Real Cepstrum.
Depending on your application's context, following up the cepstrum further might provide some more insights or uncover work that has already been undertaken in your field around what you are looking at.
Finally, I would also like to underline what Matt.L says, because you might realise that this "oscillation" is only described by a handful of harmonics and there might be a relationship that captures it that can be optimised by the data itself. This would be very useful, rather than trying to recover just one (dominant (?)) frequency.
Hope this helps.


Moving average filter for sound in MATLAB y[n] = b0x[n] +b1x[n-1]

I have this MA filter in matlab as following

 load handel;
 x = y;

 b0 = 0.5;
 b1 = 0.5;

 N = length(x); % Length of input signal
 y = zeros(size(x)); % Allocate space for outut

 y(1) = b0*x(1); % First sample, assuming x(0) = 0
  for n=2:N % Remaining samples
    y(n) = b0*x(n) + b1*x(n-1);
  end;

 soundsc(y,Fs);

I was asked to listen to the signal and describe it , well I have run the load handel without filter in matlab than i have run the signal with this filter , it sounds like that with filter the signal sounds more demped and less noise compared without filter and when i changed the b1 to b1 =-0.5 it hears only the descant the base of sound is gone . Now I want to know the reason behind.

What is the effect of the filter on the input signal which pass through the system?
what is the effect of the filter on input signal which pass through the system if we change b0= 0.5 and b1 = -0.5?

1.) The moving average filter works as a low pass filter. Consequently, it will attenuate the high-frequency part of the signal and passes the low-frequency part of the signal. Therefore, the output signal should sound flat.

2.) If we change b0= 0.5 and b1 = -0.5, the filter works basically as a differentiator, subtracting one sample from another sample. In this case, the filter works as a high pass filter, attenuating the low-frequency part of the signal. Therefore, when listening to the output signal, the high-frequency part of the input signal should appear amplified.


Color Artifacts in Fourier Transformed Image I am working with images and Fourier transforms.  I am trying to understand what might be causing some artifacts in my output image.

I am starting with a 512x512, RGB image of Lenna.  I FFT the image, and run a low-pass filter on it.  Unchanged for everything inside the circle, 0 for everything outside the circle.  Here is what the mask looks like:

So far so good.  When I inverse FFT back to an image, I get bright green artifacts:

These look out of place to me, what might be causing them?  Is it the sharp cutoff in my filter?  I would have guessed that would creating ringing artifacts, but nothing like a sharp color blob.  Also why green only?  It's not like there's a lot of green in the original.

I'm using numpy on python 2.7.6 in 64-bit Linux, if that matters.

Some things you should keep in mind:

Your image is real, keep your DFT Symmetric.
Numerical issues might cause that after the Inverse DFT imaginary parts are lefts.
Use either abs or real to get rid of them.
One should understand what it means to "Zero" out by mask DFT values.
It is not the "Ideal LPF". It causes Ringing / Halos due to Gibbs Phenomenon. One should use better designed LPF.

I hope it helps.


Considering the FFT of Real & Complex Signals I've been implementing a website to perform the FFT of various signals, real & complex.

Examining the first example, a real signal $x[n] = 10 cos(2\pi\times4n)$, I got the following FFT:

Which was exactly what I expected - two nice peaks of half amplitude at $\pm 4$

So I then extended to the FFT of the complex signal $x[n] = 10e^{j2\pi\times4n}$ (or equivalently $x[n] = 10cos(2\pi\times4n) + j10sin(2\pi\times4n)$). This is shown in the time domain as;

However, this time I got the following as the FFT, however, it now only has a single peak at $+4$, rather than the $\pm4$ mirrored peak I was expecting, and received in the real signal.

After reading through a wide number of articles relating to the conjugate symmetric property of the FFT with regard to real, even signals that states $x[n] => F*(\omega) = F(-\omega)$.

However, this doesn't help me understand what the FFT of a complex signal should look like - i.e. if my figure is correct?  How should an FFT for a real/complex/imaginary signal appear in terms of mirroring and symmetry?

NB: My own way of justifying this so far is that for a real signal, $cos(2\pi ft) = \frac{1}{2}(e^{2\pi ft}+e^{-2\pi ft})$ which creates negative frequencies - hence the mirroring, whilst the FFT of $e^{2\pi ft}$ directly does not...however I may be completely wrong!

There's nothing wrong here -  complex sinusoids like your signal really have only one peak in frequency domain! This is the fundamental idea of why we use the Fourier transform for periodic (even complex) signals.

You can think of it this way:

the cosine has two peaks, one at +f, the other at -f. That's because Euler's formula actually says $\cos x = \frac12\left(e^{ix}+e^{-ix}\right)$. You can see the two complex sinusoids that lead to your two peaks.

Now, doing the same for the sine, we see that $i \cdot\sin x =i\frac1{2i}\left(e^{ix}-e^{-ix}\right)$. That means, the sine kills one of the two peaks of the cosine, if you will so.

Everything's alright :)


Potential issues arising from too stable discretization When numerically simulating a system, usually some kind of discretization is necessary, obtained by some kind of z-transform, such as, for instance, the bilinear transform $s\mapsto \frac{2}{\triangle t}\frac{z-1}{z+1}$, which is kind of the same as using the Trapezoidal rule to approximate integration.
It also has the very nice property that it preserves BIBO-(in)stability from the s into the z-domain. However, when it is implicitly used as a differentiator, there is a very real possibility of unwanted numerical oscillations.

Now take, for example, the Euler backwards method, corresponding to the substitution $s\mapsto \frac{z-1}{\triangle t}$. This method reduces numerical oscillations quite effectively, which would (ignoring for now the slower convergence) seem to be quite a big advantage. However, if we do look at the BIBO-stability areas before and after the z-transformation, this method leads to an effectively increased stability region, i.e. there are some transfer functions that are unstable in the s-, but stable in the z-domain.

1) That the latter method is able to dampen numerical oscillations seems intuitively connected to a greater "numerical stability". However, as I understand, not all systems where numerical oscillations occur after using the bilinear transform are inherently unstable. A simple RL-circuit with transfer function $\frac{U}{I}=\frac{1}{R+Ls}$ will a current source attached will, for instance, show numerical oscillations when the source is suddenly switched on or off. So are these phenomena - increased stability domain after the transformation, and damping of numerical oscillations - connected?

2) In literature, I do find a lot about choosing a suitable z-transformation such that the resulting discrete system is stable. However, there was no regard given to possible cases where the original system might have been unstable. To me, this seems like something of an oversight - when simulating an unstable system, I would expect the discrete system to also be unstable. Because, well, the output signals for stable and unstable systems would possibly differ a lot. Have I misunderstood something there? Or is there a reason why one can safely assume the systems in question (namely, electrical power systems) to be stable?

I do know that many algorithms (say, eg, EMTP and all related) usea mix of the above, but while this decreases the difference in stability regions, it is still there.

Edit:
By numerical oscillations using the Trapezoidal rule I mean phenomena such as those outlined in http://www.ece.uidaho.edu/ee/power/ECE524/spring14/Lectures/L39/numerosc.pdf , or similar (found using google, there seem to be different texts about the same).

The Simscape Power systems Toolbox from Simulink uses the Tustin Backward-Euler method to fix this problem.

I hate to simply post links, but here it is https://www.mathworks.com/help/physmod/sps/powersys/ug/simulating-discretized-electrical-systems.html

"To perform a discrete simulation, open the powergui block and set Simulation type to Discrete, and specify the sample time. The electrical system is discretized using the Tustin/Backward Euler (TBE) method. This method combines the Tustin method and the Backward Euler method. It allows you to simulate snubberless diode and thyristor converters. It eliminates numerical oscillations seen with the Tustin method and provides better accuracy than the Backward Euler method.

The TBE discretization method combines the accuracy of the Tustin method and the numerical oscillation damping property of the Backward-Euler method. It allows you to simulate power electronic circuits with virtually no snubber, using purely resistive snubbers with a very large resistance value resulting in negligible leakage currents."

EDIT : From what I understand, they generally use the tustin method but switch to backward Euler when there is a quick variation such as a switch being open in series with an inductor for example.


What is DC level in white gaussian noise? Am studying unbiased estimators and keep seeing this term "DC Level". What is the expansion of DC and what is a DC Level? Even the Wikipedia page on WGN says nothing about it.

DC level in signal processing refers to the average or the mean value of a signal. So a zero-mean signal will have an average value of zero over its domain of definition.

When this concept is extended to the stochastic processes, we shall distinguish between time average and ensemble averages. A stochastic process can be viewed as an ensemble of time functions and a sequence of random variables. When a stochastic process is assumed to be ergodic, then its time averages over a sample function equal its ensemble averages over the random variables. In such a case the zero-mean random process creates time functions which have zero averages, hence zero DC levels.


Autocorrelation - Stochastic vs deterministic processes On one hand, there is the statistical definition:

$R_x(\tau) = \frac{\operatorname{E}[(X_t - \mu)(X_{t+\tau} - \mu)]}{\sigma^2}$

and on the other hand, there is this other definition (e.g. see Dilip's comment here):

$R_x (\tau) = \int_{-\infty}^{\infty} x(t) x(t+\tau) {\rm d}t$

Why are these two definitions different? Is the first one defined for stochastic processes, while the second one defined for deterministic signals?

If so, does the Wiener-Khinichin theorem hold for the deterministic variant? (the Wikipedia article only discusses it in the context of a stochastic process)

You are correct that these two definitions are for different processes.

$R_x\left(\tau\right) = \frac{1}{\sigma^2} \mathbb{E}\left[\left(X_t - \mu\right)\left(X_{t+\tau} - \mu\right)\right]$

is defined for a stochastic process $X$ distributed according to a law $\nu$.

Whilst for a deterministic, finite-energy signal, $x\left(t\right)$ the autocorrelation is:

$R_x\left(\tau\right) = \int_{-\infty}^\infty x\left(\tau+t\right)x\left(\tau\right)dt$

The Wiener-Khinchin theorem will hold for any square-integrable deterministic $x\left(t\right)$ that is:

$\int_{-\infty}^\infty |x\left(\omega\right)|^2 d\omega < \infty$

so that the Fouier integral exists and converges in the $L^2$ sense (i.e. signals with finite energy).

The two definitions are different because the integrals used in the definitions are different. The deterministic case implicitly assumes a Riemann integral - that the process is nice in some way - it is continuous, and if differentiable the derivatives are bounded etc. However for random processes, $dt$ is difficult to define - the process could jump from $-\infty$ to $\infty$ over any inteval I could possibly define, so the Riemann limit isn't appropriate.

Instead we use the Lebesgue integral, where we can define the expectation of a random process: all that's required is that the range is smooth (for example $\mathbb{R}$), and we can still take the expectation with respect to some measure, modulo some restrictions on the random process (so that the integral converges in a monotone or dominated sense).

Integration of Random Process is a tricky business, and the definitions are written differently to keep people mindful of what they are working with.

One other note - the deterministic case isn't so interesting, as we can take the FT of a deterministic signal: the whole point of the WK theorem is for random processes - where the Fourier integral may not exist. I find it a particularly satisfying theorem, and the work to understand it definitely worthwhile.


What is meant by optimal estimator and how to determine optimality? Considering an estimation problem of estimating a scalar deterministic parameter $a$ from the observations $y$ which are corrupted by randomvariable $w$. The observations are $y[n] = a + w[n]$

Least Squares estimator can be used to estimate $a$ when $w$ is a White Gaussian Random Variable. This estimation method is known to be optimal. Why?

What if $w$ is from Poisson Distribtuion or some other non-gaussian, then would the estimator for $a$ be better or worse than the one found using $w$ as Gaussian r.v?

So Least squares estimator is as it literally - the estimator which brings the mean square error to minimum. In the case of Gaussian white noise it has a simple and analytic solution. I recommend you develop it yourself, if your'e comfortable with matrices calculus it is not that hard.

You can generate tons of estimators by defining different cost functions. Other popular cost functions are l1-norm ($\sum |x_i -  y_i|$), likelihood ($P(y|x)$).

It's hard to define "better" if the case isn't the same, For example, does the case of very high variance poisson process is the same as the same variance gaussian process? I believe it is not, and in addition one is non negative integer while the other is real (and you can use it for you advantage).


Image processing and compression I have two questions related to this image;

Is it possible to  identify any spatial domain features in the frequency domain?

My understanding: what I read around they are mostly the same as spatial domain but shows $\sin$ or $\cos$ ripples for every position in the frequency domain but not sure if it correct something like this:

Second question: when I have information about the images, the frequency spectrum and histogram, how do I select the compression method? Do I look for specific pattern or do I need something else? Or is it something entirely else?

For the first question, first. In short, yes you can, somehow. It depends a lot on how you process images, and the nature or morphology  of the images you analyze.

First, it is well etabished that amplitude and phase spectra of natural images are often difficult to analyze (see classical textbooks). I would emphasize three classical issues: nature of the processing, of the images, of the features.

First, the processing itself. Here, you can spot a bright cross on you FFT magnitude spectrum. It is probably due to the periodicity assumed by the FFT. This bright cross, which has a lot of energy, might obfuscate more subtle details. Therefore, some processing steps can be useful  before interpreting features: filtering, windowing, cropping, etc.

Second, the morphology of the data. Fourier is great at analyzing stationnary data, either deterministic or stochastic. When the data has  non-stationnary aspects, discontinuities, the power of approximation of Fourier methods decreases a lot. This does not entail that Fourier is bad, but its interpretability becomes way more complicated. However, in the Fourier plane, you can spot diagonal lines, or cones, or bumps, that concentrate spatial information. Here, with proper preprocessing, you could expect more concentrated information from the mixed sines. For instance, here, the bright line from the top-left quadrant to the bottom-right one is likely to be related to the left-to-right upward edge of the hat, and the parallel stripe patterns above.

Third, the features. Here, you are looking at visual features from the amplitude spectrum only. A lot of information also resides in the highly non-linear part of Fourier: the phase. Moreover, crafted features can be more involved, like look at at specific scales  (eg wavelets or related tools) or layers in deep learning, etc.

There is a lot of ongoing works around your question. My take is that complex-valued, scale-related features can be very efficient. But the margin, you know...

Related posts:

Neural Networks and Complex Valued Inputs
Why are wavelet transforms (Multi Resolution Analysis) used more in practice for compression rather than Fourier series?

For the second question. Early image compression (1970's, 1980's) have used by the amplitude spectrum and the phase. And the story was: for natural images, phase requires more bits than amplitude. Since then, non-complex (eg real) transforms have been more efficient (DCT-II, biorthogonal wavelets). But I do believe complex transforms will have their word again.


Linear Predictive Coding - Meaning of spectrally Flat I am a student with very little background on signal processing and am studying the Linear Predictive Coding algorithm. I came across the term spectrally flat in this article,
https://ccrma.stanford.edu/~hskim08/lpc/ and am not sure what it means. I know that spectral is talking about frequency.

  There are two possible signals for the source: an impulse train or
  random white noise. These signals model pitched sounds and
  plosive/fricatives respectively. The common characteristic for both
  impulse train and white noise is that they are spectrally flat; all
  spectral information is modeled in the filter.

So my question is:

1) what does spectrally flat mean ?

2) why does an impulse train have a flat spectral ?

1) what does spectrally flat mean ?

Spectrally flat means that the signal has energy at most / all frequencies, and this energy is roughly the same level across the frequency range.

  2) why does an impulse train have a flat spectral ?

As you're probably aware, the Fourier transform of an impulse train is also an impulse train. That means that the claim by the link that the impulse train is spectrally flat is not correct.

The reason an impulse train is used in modeling speech is because it gives a tonal quality for non-fricative sounds (mostly vowel sounds). That tonal quality is hard to generate by filtering if the source being filtered is white noise.


Why sampling rate is not an option with normalized frequency in Matlab Filter Designer My Matlab version is 2017b. When I use Filter designer with normalized frequency specifications, I cannot set the sampling rate. May you tell me why the sampling frequency cannot be specified? As I know, the normalized frequency mapped $[0,\pi]$ to $[0,1]$, so the sampling rate is not dropped out.

...the normalized frequency mapped $[0,\pi]$ to $[0,1]$...

The normalized frequency $F$ is defined as $f/f_s$ where $f$ is the "real" frequency (i.e. not the normalized one) and $f_s$ is the sampling frequency. The map is thus from $[0, f_s]$ to $[0, 1]$ (or to $[0, 2\pi]$ in normalized pulsation $\Omega = 2\pi F$). When using normalized frequency $F$, you thus specify frequencies as a percentage of the sampling frequency. Hence, the sampling frequency itself is not directly needed.


On computing the number of independent values from a signal's bandwidth? Can anyone help me understand where the value of $2$ comes from in the following statement?

A time function t whose frequency response is white out to a bandwidth B has
about 2 · B · T independent values when measured for the time interval T.

This is in the context of a serial stream of digital (e.g. binary) data. The function t represents this data stream.

It comes from the Nyquist rate.  The Nyquist sample rate for a signal with bandwidth $B$ is $2B$.  In other words, that is the lowest sample rate that will contain all of the information in the signal without distorting it (i.e. without losing information).

The question implicitly assumes that the sample rate is $2B$, and thus that in time $T$ you will get $2B*T$ independent samples.  If the sample rate was greater than $2B$, then the signal would not be "white" anymore and the samples would have some correlation.


Converting floating point number to Q-notation fixed point number in C/Java? I am trying to convert from a floating point number to a Q-notation signed 16-bit fixed point number (e.g. Q15). How can this be done in C or Java? I searched online but could not find any available code.


Measuring similarity of music, using lossy compression I have two sound clips $C_1,C_2$ which are pretty similar.  I'd like to measure how perceptually similar they are, i.e., how similar a human would perceive the two to be.  Is there a way to use a lossy compression algorithm (e.g., a MP3, AAC, or Ogg Vorbis encoder) to compare the two clips?

It occurs to me that that audio compressors already contain a good deal of knowledge about psychoacoustics and human perception of sound, built into them.  Is there a good way to use them to measure how similar the two clips are?

Maybe something like $L(C_1 || C_2) / (L(C_1) + L(C_2))$, where $L(x)$ is the length of the compression of sound clip $x$, and $C_1 || C_2$ is the result of concatenating the two clips?  Or maybe find the highest bitrate such that
$F(C_1)$ is close to $F(C_2))$ via a simple metric (e.g., L2 norm applied to the FFT spectrum), where $F(C)$ is the result of compressing $C$ at that bitrate followed by decompressing it?  Or something else along these lines?  Has anyone studied this?

If it is relevant, the two clips are fairly similar: one was obtained by a transformation of the other.  They are aligned in time and have the same length.  Each is relatively short (at most a few seconds).  I've done some searching but haven't found any reference or research paper that discusses this sort of approach -- maybe I just haven't found it yet, though.

Is there a good way to use them to measure how similar the two clips
  are?

This is a very sane intuition but something I have rarely seen in the literature - outside the peripheral idea of doing feature extraction from compressed streams (and the motivation here is just the reduced computational burden).

I think the major reason is that audio codecs do not store a perceptual representation of the audio data. Instead, they store a very coding-centric representation of audio - entropy-coded, quantized, coefficients of a transform - chosen so that the distortion between the original and the encoded signal falls below the threshold of perception. In a sense, looking at a compressed audio stream does not tell you anything about how a human would perceive the audio signal. It just tells you that anything that would fall below the perception thresholds has been taken away. The codec "carves" around the meaningful audio signal, but never touches it.

  where L(x) is the length of the compression of sound clip x, and C1||C2 is the result of concatenating the two clips?

This makes sense on strings - I think this idea is found in some papers by Cilibrasi or Vitanyi - but this is not how audio codecs work. The audio codecs you mention are designed to process streaming audio with low latency, so they only "see" a short window of the signal at a time and do not attempt to eliminate long term redundancy. Just like JPEG for images. Copying 10 instances of the same note one after the other will result in a file 10x as large. Audio codecs have a fixed bit-rate, or an adaptive bit-rate that is adapted according to the content of the narrow window of signal that the codec sees.

The redundancy eliminated by an audio codec is very short-term.

Your idea would require a much more sophisticated "objet oriented" codec that attempts to decompose the audio into a hierarchy of objects, such as notes. This idea is still in its infancy in the academic world.

Another way to look at it is to check if the decoding process bears similarity with a generative model of the signal. This is the case for some string compression algorithms, where the compression process recovers something akin to a generative grammar of the string (Sequitur, factor oracles...). This is the case for LPC-based speech codecs in which the LP coefficients embed knowledge about the articulation, and the residual knowledge about prosody and voicing/unvoicing - the compressed stream can be seen as data for a speech synthesizer. But this is not the case for general purpose audio codecs like the ones you mention - the decoder bears no similarities with the music production process.

  Or maybe find the highest bitrate such that F(C1) is close to F(C2)) via a simple metric (e.g., L2 norm applied to the FFT spectrum)

You might do the opposite - look at the bitrate under which the clips become indistinguishable, but ultimately it wouldn't be that much easier than comparing sinusoidal representations pruned by a perceptual function.


DFT of pure sinusoidal wave I'm writing a program in which you can synthesize waves by adding to a sound's Fourier transform, and then inverse the transform to get the modified sound. In order to do this, I need to know what to add to the DFT to synthesize a pure wave. I've tried to learn about Fourier synthesis from many sources, but they all talk about the Fourier series instead of the Fourier transform, and they all say that for a pure wave all you need is a value in the coefficient of that wave's frequency, and 0 everywhere else.
But the Fourier transform being a continuous function. And you can see here what the Fourier transform of a pure 12KHz sine wave looks like (on a logarithmic scale). As you can see, it's not just an instantaneous peak and 0 everywhere else.
So I tried to do the math myself. To put it formally, say we have a sampled signal $S_t=A cos(\frac{-2 \pi k}{N} t + \phi)$, where , $N$ is the number of samples in our signal, $0 \le t < N$, $A$ is an amplitude  $0 \le k < \frac{N}{2}$ determines the wave's frequency, and $\phi$ is the phase of the wave. The DFT of this signal would be a sequence $F_0,...,F_{N-1}$ where:
$$ F_r = \sum_{t=0}^{N-1}{{S_t}e^{\frac{-2 \pi i t r}{N}}} $$
My hope was that since our signal is a pure wave, there will be a direct formula for calculating $F_r$ that can be computed much quicker than by directly doing the math above, and that that formula will be the shape of that spike in the image I linked, because that's what I'm really interested in. For the sake of simplicity, I assumed that $\phi=0$, since once I figure out how to solve that case I'll probably be able to generalize it. I won't bore you with all the math because it's not important and you probably already know the answer. What I got in the end is that $F_k = F_{N-k} = \frac{NA}{2}$, and for all other $0 \le r < N, F_r = 0$.
Shocker. So in the end I found that it does equal to 0 everywhere except the frequency of the wave. But then I don't understand where that spike shape from the graph I linked is coming from. I'm still after it though, because just adding to the one sample which corresponds to the frequency I'm interested in doesn't work for me (I've tried).
So I guess my question is: how come plotting the Fourier transform of a pure wave doesn't actually produce zero in all samples except the one corresponding to that wave's frequency, and what's the formula for that spike that you do see.

First of all, welcome to DSP SE.
What you see in the image you have linked is termed (spectral) leakage. When you are dealing with the Fourier series you deal with a periodic continuous function which is "decomposed" into a (possibly) infinite sum (series). Then, when you go to the Fourier transform, you have a non-periodic function (which you could possibly presume to be periodic in some interval if, for example, it starts and ends on the x-axis) which is decomposed into an infinite sum of spectral components (you have energy in every point on the spectrum).
Before moving further on, please note that this is a very brief explanation about the Fourier series and Fourier transform with a lot of details omitted and the mathematical interpretation presented here in more of a convenient way than the absolutely correct. For more information on the topic, I would strongly suggest the textbook Mathematical Methods for Engineers and Scientists 3 - Fourier Analysis, Partial Differential Equations and Variational Methods by Dr Kwong-Tin Tang (the first part of the book is relevant here).
Now, if you go from the continuous domain to the discrete domain (we will deal with the discretisation of the free variable here, which is time and not the amplitude as is the reality when you deal with digital signals) you move from the Fourier transform to the Discrete Fourier Transform (DFT). Once again, omitting the "technical" details, in DFT when the signal contains frequency (spectral) components that do not make an integer number of complete periods in the duration of the signal to be analysed, you end up with the phenomenon called leakage (see link above). This is due to the "nature" of the transform, in which the signal is assumed to be periodic. Thus, when you have spectral components that do not complete an integer number of periods in the duration of the signal, you will end up with an amplitude (for this component) other than zero at either the end or the beginning of the signal. If you try to "copy and paste" the signal before and after your original signal (in order to make it periodic for the purpose of the transform) you will realise that you end up with some discontinuity at the point where the original and the "copy-and-pasted" signal starts. In order to reconstruct this discontinuity, you have to introduce an infinite amount of spectral components, which show up as energy in frequencies around the main spectral component (in the case you only have a "pure" (co)sine wave). You can see an example of a single frequency in the picture below.

Now, another concept you will most probably encounter quite often (if not always!) when dealing with DFT is the term window. In the simple case where you apply no windowing function (bear with me a wee bit more for an explanation) to the signal is like applying a "rectangular" window where you multiply all the values of your signal with one, thus, effectively doing nothing to them.
If you apply a windowing function to your signal you effectively suppress the signal at its edges in order to avoid the appearance of discontinuities such as those visible in the second plot on the right-hand side of the image above. Thus you somehow suppress the leakage effect. Please keep in mind that the simplest way of applying a windowing function is to element-wise multiply the samples of your signal with the windowing function in the time domain.
The signal on the right-hand side of the above image after the application of a Hann window would look like in the image below

And for a "clearer" representation of the spectrum before and after the windowing you can see this image

You can see that the peak is not so well localised but the leakage has been constrained to a smaller "neighbourhood" than before. For more information on windowing and window functions, you can have a look at the Wikipedia page (link above at the word "window") and any elementary Digital Signal Processing textbook such as the Digital Signal Processing - Principles Algorithms and Applications by Proakis & Manolakis or Introduction to Signal Processing by Orfanidis (which is freely distributed).


Usefulness of Matrix Notation for Linear Periodically Time Variant Transformations I am quit puzzled about a notation I found in Gardner's 1986 book titled Introduction to Random Processes. It is in the chapter on Cyclostationary Processes at section 12.4 which pertains to Linear Periodically Time-Variant Transformations:

  (...) many modulation systems can be modelled as the scalar response
  of a multi-input LPTV transformation with stationary excitation. This
  includes amplitude modulation, phase and frequency modulation,
  quadrature amplitude modulation, pulse-amplitude modulation, all
  synchronous digital modulations such as phase-shift keying,
  frequency-shift keying, and so on. Consequently, the study of
  cyclostationarity is facilitated by general formulas that describes
  cyclic spectra in terms of the parameters of LPTV transformations.
  This includes cyclic spectra that are generated by LPTV
  transformations of processes that exhibit no cyclostationarity, as
  well as cyclic spectra that are transformed by LPTV transformations of
  processes that do exhibit cyclostationarity. Let us consider the LPTV
  transformation $$ Y(t) = \int_{-\infty}^\infty
\textbf{h}(t,u)\textbf{X}(u)du $$ for which $\textbf{X}(t)$ is a
  (column) vector excitation, $Y(t)$ is a scalar response, and
  $\textbf{h}(t, u)=\textbf{h}(t+T, u+T)$ is the (row) vector of
  impulse-response functions that specify the transformation.

The subject of LPTV Transformations is new for me and I don’t understand the usefulness of using non-scalars for $\textbf{h}(t, u)$ and $\textbf{X}(t)$.

Is this a commonly used notation?

Why should a matrix product be used?

I have been trying to make sense of this for some days now and I have really made no progress whatsoever. I hope someone here has an idea...

With thanks!

I don't have the book, so my answer unfortunately cannot be very detailed. However, just based on the excerpt, this could be useful if you are modulating a number of different data streams $\mathbf{X}(u)$ onto a single transmit signal. For example, this would be the case if you have a number of frequency bins that you will use for different data, as would be done for OFDM (orthogonal frequency devision multiplexed) transmissions or DMT (discrete multitone) as is used in DSL modems. I'm personally curious why $\mathbf{h}(t,u)$ is not a matrix and $Y(t)$ a vector, since it is perfectly plausible to have a system with multiple transmitting antennas.


Regarding Bode plots; $H(s)$ and $H(j\omega)$ In circuit analysis, I understand the use of Laplace Transforms to obtain the impedance of a linear RLC circuit, ie transforming from the time domain to the frequency domain. In most texts I have seen and classes I have taken, the Bode plot looks at the rational transfer function of $H(j\omega)$, which is effectively the Fourier Transform, and not $H(s)$.

Since $\omega$ is considered to be the variable that represents frequency, what is $\sigma$ if $s = \sigma + j\omega$? Does $\sigma$ have a physical interpretation?

Note that each complex pole $s_{\infty}=\sigma + j\omega$ of the transfer function $H(s)$ contributes to the system's impulse response a complex exponential of the form

$$e^{s_{\infty}t}=e^{\sigma t}e^{j\omega t}$$

The term $e^{\sigma t}$ is real-valued and represents exponential damping (assuming that the system is causal and stable, i.e. $\sigma<0$), and the complex term $e^{j\omega t}$ represents an oscillation at frequency $\omega$. So $\sigma$ is a damping constant.


Design of the low pass filter when decimating I am filtering an audio signal into various sub-bands in C and would like to decimate the signal by 2 in the lower frequency bands to reduce processing overhead. I get that this involves low-pass filtering the signal and then throwing away every other sample, but how do I go about designing the low-pass filter optimally?

At the moment I have a signal sampled at 44.1 kHz. Say that for signals in bands limited below 8 kHz, I low-pass filter with a 2nd order butterworth filter with a cutoff of 9 kHz and then decimate. This works fine as I filter well below the new Nyquist limit of 44100/4, but it feels very unscientific.

What is the best filter type for this purpose. Looking at the matlab documentation, the decimate function uses an 8th order Cheby Type 1 filter by default, but I have no idea why. Is there some rule-of-thumb that I'm missing?

You can determine the "best" filter by analyzing what type of errors is created and how much of each error your specific application can live with.

If your original content is already perfectly band limited, than you don't need a low pass filter at all and can directly throw away the extra samples.

The decimation low pass filter typically determines how the band limiting of the content is done. Here are the factors to consider

Aliasing: since the low pass filter is never perfect, the decimation will
create some amount of aliasing. The amount depends on the steepness of the filter and the energy content of the signal at the frequencies that need removing.
Pass band: The low pass filter isn't perfectly flat in the pass band so there is some amount of ripple you need to live with. The pass band must be below the new Nyquist frequency. The higher the filter order, the closer you can get and the wider the pass band will be.
Phase response: Any IIR filter will change the phase response of your signal. The steeper the filter, the higher the maximum group delay. If need to maintain the phase, you need to use a linear phase FIR filter (ideally in a polyphase configuration).

The choice of filter is a trade off between pass band width, pass band amplitude distortion, pass band phase distortion, aliasing noise and computational complexity.If you have a spec for all of these, you can design the "best" filter. An example would be "pass band +-0.5dB up to 90% Nyquist, aliasing at -80dB or below, linear phase".


PSD of a windowed signal Suppose PSD of a signal $x(n)$ is $ PSD_{xx}(\omega)=1$. Now if a window function $w(n)$ of length $ 0\le n \le N-1$ is applied to $x(n)$, what will the new PSD be?

Windowing is multiplication in the time domain which corresponds to convolution in the frequency domain, so the spectrum  of the windowed signal is the spectrum of the original signal convolved with the spectrum of the window. The PSD is simply the magnitude squared of the spectrum of the windowed signal, i.e.
$$P(\omega)=\left \| X(\omega)\ast W(\omega)  \right \|^{^{2}}$$
where $\ast$ is the convolution operator.
Hence the desire to have the spectrum of the window look like an impulse as that has the least impact on the PSD.


MFCCs and chord recognition This is my first time working with Mel-frequency cepstral coefficients and I have two questions....

Are some of the 13 coefficients strictly related to some sound properties like frequencies, amplitude, pitch or something ? For example if I'm interested only in sound pitch, are there some coefficients that carry that information and I can exclude the rest of the coefficients that carry other information to reduce the size of the vector ? Or those coefficients means nothing individually and they must be used all together as a whole ?
I'm trying to learn Chord recognition methods. I was thinking for calculating MFCC from a short Wave file where chord is played. (It is just a chord played on piano, there is no ambient noise or other instruments). And after I get the 13 MFCCs, cluster them to a centroid which describes the appropriate chord. (24 centroids, 12 for minor chords, 12 for major). I'm just wondering, are the MFCCs good descriptor for the present pitches in the sound ?

Thanks in advance

MFCC's are related strictly to spectrum shape and what's more - these are relaying on the mel scale. Additionally each coefficient is a DCT value from fitted cosinusoids to log-energies. Thus interpretation of pitch would be troublesome for you.

Chord recognition is usually based on chromagram calculation or so called chroma features. Check out work and publications of guys from QMUL, i.e: Real-time Chord Recognition for live performance - in here you will find another references.

You might also find their software useful: Sonic Visualizer together with lot's of so called Vamp plugins. You can check out the code and even use command line tool to extract features you want.


Why is the frequency response different based on the amplitude of the impulse I feed into the filter? I am working on a 2-pole IIR high-pass filter.  The filter is fixed point and it is causing me some issues.  I believe I'm just misunderstanding something.  I am implementing everything in Octave for now, but eventually it will be implemented in a 32-bit processor.  The corner frequency is very low, $2.7\textrm{ Hz}$.

When I test the filter by passing in an impulse and then doing an fft on the output, I get different results based on the magnitude of the impulse.

First, I made the filter using signed 64-bit integers just for testing.  So the impulse I feed into the filter looks like this:

Fs = 51200; %sample frequency
impulse = int64( [2^31, zeros(1,Fs*10-1)] ); %multiply by 10 to get 0.1 Hz resolution

Here is the fixed-point impulse response compared against the floating-point impulse response.

If I rerun the code but with a larger impulse I get the desired frequency response.

Fs = 51200; %sample frequency
impulse = int64( [2^63, zeros(1,Fs*10-1)] ); %multiply by 10 to get 0.1 Hz resolution

They match so closely you can barely see the red, floating-point response.

Why is there this difference just from changing the amplitude of the input?
I need to implement this in 32 bits, but when I do, the amplitude of the input is restricted to $2^{31}$ so I end up with the distorted impulse response.  How can I fix this or work around this?

The lower amplitude the input signal has compared to the resolution of the filter coefficients the weaker is its ability to shape the frequency of the input signal. It also depends on the span of the filter coefficients.

You get the same result if you feed the filter with 32-bit coefficients with a signal having an amplitude covering only the lower 16-bits. If you want to work around this effect you need to amplify the input signal.


Gradient of Total Variation (TV) Norm in Total Variation Denoising In this link, it says that the gradient is as follow

  The Gradient of the TV norm is
  $$
\mathrm{Grad}J(f)=\mathrm{div}\left(\frac{\nabla f}{\lVert\nabla f\rVert}\right).
$$

From this other link, it also mentions that the derivative is as follow

  Formally, we can write $\partial J(u)=-\mathrm{div}\left(\frac{\nabla u}{\lvert \nabla u\rvert}\right)$

I know how to calculate divergence, but I don't understand how the gradient of total variation is related to divergence.

I am by no means an expert on total variation, however I think you should check out this Wikipedia page. It doesn't directly answer your question, but I believe the lemma below illustrates the relationship between total variation and divergence.

There, it gives a lemma that follows from the Gauss-Ostrogradsky theorem and provides a proof for it,

$\int_{\Omega}f\,div(\phi) = -\int_{\Omega} \nabla f \cdot\phi$.

I think for some intution about divergence, it might be helpful to read Wikipedia's explanation.


Why does increasing FFT length (narrowing bandwidth) not decrease Noise Power per bin? I was under the impression that zero-mean white guassian noise had a constant Power Spectral Density which means that smaller bandwidths should decrease the power in that band. I'm getting surprising results for average power in each FFT bin, it's staying the same irrespective of how long I make the FFT. Am I not measuring power correctly?

Here is a reduced for of my MATLAB script (should be Octave-compatible). Please excuse my odd habit of time vectors going up-down and frequency and misc ones going left-right.

%pkg load communications
% ^UNCOMMENT^ if using OCTAVE
close all

Q=8;        % Lowest power of 2 to iterate though in FFT loop
R=12;       % Highest power of 2 to iterate though in FFT loop
N=1;        % Linear PSD of noise
c=0;        % Loop counter

K=2.^(Q:R)  % Loop vector

for k=K
    n=sqrt(N).*(randn(k,1)+j.*randn(k,1))./sqrt(2); % complex noise
    P=(abs(fft(n)')).^2./k; % Noise Power in frequency domain
    c=c+1;
    e(c)=mean(P); % Mean power measurement across bins
    %v(c)=var(P);
    figure
    plot(10.*log10(P))
end
e
%v

I feel like I'm missing something fundamental about this. Does FFT binning not truly band-limit signal power because of the way the algorithm works? I'm basically trying to get an output analogous to a spectrum analyzer in that is shows something proportional to true power.

A spectrum analyzer typically has more than one mode. One mode normalizes the frequency output such that a pure tone of amplitude “one” with frequency centered on a bin center has an amplitude of “one” in the spectrum display.

For a pure tone, nonrandom, centered on a bin frequency of amplitude “one” , the corresponding uniformly weighted  N sample DFT bin magnitude will be “N”, so all the bins are reduced by dividing by N.  (before magnitude squaring)

An actual Spectrum Analyzer will include a fudge factor correction for nonuniform windows, and average multiple DFTs .

I haven’t seen too many spectrum analyzers that don’t have a DFT under the hood.  The last analog swept frequency analyzer I used was a very long time ago. I have seen a few modern versions at frequencies where converters aren’t good enough but DFT versions dominate most of the market.


How to Pick a Journal for Publishing a Research Work? I got my MS in EE way back in 1994 and have not presented at a conference, or published in a journal before, but I had and still have an intense interest in the applied math that I learned. I have been working my own research project for a few years, and I want to publish a paper on my work. What I did is develop an algorithm similar to a discrete Kalman fiter, and it gives better estimates under certain conditions. I am looking for advice on what journal would a good choice. It would take to long for me to learn LaTeX (is that what it's called), and I don't have the education needed to develop the sophisticated math that I see in many journals.

Getting work published in a high impact factor journal is quite a slog. If you're starting from not having presented at a conference etc. before, that's going to be an even bigger ask.

However, it's not impossible.

For IEEE journals, short notes about signal processing topics can be sent to the IEEE Signal Processing Letters journal.

However, starting there might be too much.

What I'd suggest is to look at starting with an online journal such as arxiv or PLOSone.

Even before submitting there, I'd suggest posting a question here that might be relevant to the situations where the algorithm might be better than others, and also posting an answer with your algorithm.

The aim with any of these submissions is to get feedback about your work and improve it or its presentation.


Frequency spectra of a sampled process I have a process $X$ with frequency spectra

I sample this process with sampling frequency $f_s = 2$.

What will the frequency spectra $R_z(f)$ of the sampled process $Z$ look like? I realize that aliasing occurs since the original process contains frequencies above the Nyquist Frequency (in fact only contains frequencies above the Nyquist Frequency). The process has no frequencies within the interval where no aliasing occurs, i.e, $-1 \leq f \leq 1$. Therefore I reasoned that the frequency spectra of the sampled process would be $R_z(f) = 0, \quad 1\leq f\leq 1$ but the answer is $R_z(f) = 2, \quad 1\leq f\leq 1.$ Why?

Sampling in the time domain is equivalent to periodic repetition in the frequency domain with $1/T_s$ as the period. So formally you'd have something like
$$R_{sampled}(f) = \sum_{n=-\infty}^{\infty}R_z(f-2\cdot n)$$

If you evaluate this, you actually get $R(f) = 2$. It's easier to see if you draw it: take the initial spectrum, shift it by -6,-4,-2,0,2,4,6... and sum it all graphically


What is the bandwidth of product of two signals? For signal $x(t)$ with bandwidth $B_x$ and signal $y(t)$ with bandwidth $B_y$, what will the bandwidth of the signal $z(t)=x(t)y(t)$ be?

The bandwidths simply add.

You can break down both signal into their sinusoidal components and pairwise multiply them. For every pair of sines you get the sum and difference frequencies. This is a simple consequence of the multiplication theorems for sine waves. The highest frequency of the product will be the sum of the highest frequencies in the individual signals. This is equivalent to a convolution in the frequency domain.


Shannon formula for channel capacity Shannon formula for channel capacity states that
$$C=\log\left(1+\frac{S}{N}\right)$$
If this formula applied for baseband transmission only? Is it applied with passband transmission?
Are there any changes in the formula with different types of modulation, for example, OFDM??

This formula is for baseband with real symbols in an AWGN channel.
However, it is trivially easy to adjust it for quadrature modulation or complex symbols: just double the capacity.
OFDM and other wideband modulations operate over frequency-selective channels, which have different capacity. There are capacity formulas for many channels in the literature.


How to compute total interference power at receiver antenna I would like to know what is the proper way to sum interferences.
I have, for example, let say  0.4 mW, 0.035 mW and 0.2 mW are those interferences.
Do I have to sum them resulting 0.635 mW or make a conversion to log scale, for example, dBm and add these dBm log values?
The results are different; if I convert the three interferences 0.4 mW, 0.035 mW, and 0.2 mW to dBm and then sum them, it results in -25.52 dBm, which back to linear scale is 2.8$\mu$ W

The proper way to sum the interference is to sum them in power, unless they are coherent, in which case they would sum in magnitude. Coherent would mean the interference are all at the same frequency and locked in phase. This is unlikely so if not state explicitly you would assume the are non-coherent and a power sum should be used.
To sum them in power you would simply add the powers and then convert the sum to dBm. If starting from dBm, you would convert them to power quantity first, add and then convert back to dBm as demonstrated below using the OP's quantities:
0.4 mW is -3.98 dBm
0.035 mW is -14.56 dBm
0.2 mW is -6.99 dBm
The total power would be 0.4 + 0.035 + 0.2 mW = 0.635 mW which is $10\log_{10}(0.635)= -1.97$ dBm.
If we started with the dBm quantities we could compute this same result as follows:
$$10\log(10^{-3.98/10}+10^{-14.56/1}+10^{-6.99/10})$$


Upsampling algorithm too slow I have to come up with an upsampling algorithm for a project and cannot use Matlab's built in upsample function. I came up with the following function but it just runs forever when I use a 13 second, 16-bit song sampled at 44.1kHz. Please let me know how I can improve the efficiency of this algorithm:

    x = desired_Fs/original_Fs; %desired and original sampling frequencies
    len = x * length(y); %y is the input signal
    q = x - 1;

    for i = 1:x:len
       for j = length(y):-1:i+1
            y(j + q) = y(j);
       end
       for k = i+1:1:i+q
            y(k) = 0;
       end
    end `
    upsampled_audio = y;

First, you should use informative variable names. For example, use "upsample_ratio" instead of x, "inSignal" instead of y etc. This will make your code much easier to read and understand.

Second, you should avoid using for loops in matlab if possible, and especially don't use nested for loops, since they take very long time in matlab.

Third, you overuse the variable y. Use intermediate variables to keep temporary and result values.

Forth, your algorithm concept is wrong in case the upsampling ratio is not a whole number.

Specifically for your code, start by replacing the inner loops with array operations:

vec = length(y):-1:i+1;
y(vec+q) = y(vec);
y(i+1:1:i+1) = 0;


Does each pixel of the original image has a corresponding coefficient for every scale and every direction calculated by the DWT? I want to take a pixel and find its coefficients calculated by the discrete wavelet transform for every scale and every direction. But as the DWT subsamples images, I don't see how to do that. Do I have to interpolate the coefficient matrices produced by the DWT?
In the non-subsampled contourlet transform for example this is not an issue because there is no subsampling. For each pixel of the original image I take the coefficients that have the same position in the generated images, but that's not possible with the wavelet transform because images get smaller and smaller in size.

No, by essence, as you guessed correctly. There is a direct counting argument: given a $2N\times 2N$ image, the first level of the  DWT generates  three $N\times N$ wavelet subbands. And through the pigeonhole principle, you cannot have a bijection between a $4N^2$ and a $3N^2$ coefficient set.

As for the contourlet, you can resort to undecimated wavelet versions, for which several implementations are possible. Some allow more directions than those obtained through traditional tensor product.

However, I consider that it is misleading to believe that one isolated pixel could be easily followed across scales or directions and mapped to one coefficient: this requires smoothing or directional filters taking into account weighted combinations of the surrounding pixels. In other words, one should not expect arbitrary precision on location AND arbitrary scale or orientation determination.


Terminologies: Adaptive, recursive and iterative Considering non-linear filtering technique like Extended Kalman filtering with Expectation Maximization (EM). EM is an iterative technique but what is Kalman filtering? Is Kalman Filtering called iterative approach?
Adaptive signal processing algorithms like Least Mean Square and Recursive Least square (RLS) estimate weights at each time instant and then update them. Since, RLS has a recursive nature and it falls under adaptive technique, then does this mean that adaptive is another name for recursive?
What is the difference between adaptive and iterative terminologies? Are Kalman filtering and its nonlinear versions adaptive?

Q1: In short, is adaptive=recursive=iterative?

Q2: Where is EM and Kalman Filtering approaches categorized under and why?

Q1: In short, is adaptive=recursive=iterative?

Filtering is applying a filter $f$ to an input signal $x$ to get an output signal $y$:

$$
y(t) = f(x(t), t)
$$

The filter $f$ is called iterative if its next value can be calculated from other time-dependent values of either $x$ or $y$ or both.

The filter $f$ is called recursive if it depends on $y$ for some time $s \le t$ (where $t$ is the current time).

If the way $f$ uses past values of $y$ to filter $x$ changes, then the filter $f$ is also adaptive.

For example, the filter:
$$
y(t) = x(t)
$$
is iterative, but not recursive nor adaptive.

For example, the filter:
$$
y(t) = x(t) + \alpha y(t-1)
$$
is both iterative and recursive, but not adaptive.

For example, the filter:
$$
y(t) = x(t) + \alpha(y,x) y(t-1)
$$
is both recursive and adaptive (the filter coefficient $\alpha$ depends on some function of past inputs $x$ and past outputs $y$).

  Q2: Where is EM and Kalman Filtering approaches categorized under and why?

Kalman filtering is definitely iterative and recursive as the new updates will depend on past values of the input and output (well state estimates).

The standard implementation of the Kalman filter is not adaptive.  However, there are simple extensions that make it adaptive.


Upper Bound for the DFT (FFT) Coefficients of a Bounded Signal I have a discrete signal $x[n]$ of length say $1024$, such that :
$$-1 \leq x[n] \leq 1, \qquad \forall n$$
and let $y[k]$ be its DFT.
Is there an upper bound for $|y[k]|$ ?  (i.e. for the DFT of a bounded signal)

You can obtain a bound on the magnitude of the DFT of $x[n]$ ($|x[n]|\le 1$) as follows:

$$\big|X[k]\big|=\left|\sum_{n=0}^{N-1}x[n]e^{-jkn2\pi/N}\right|\le\sum_{n=0}^{N-1}|x[n]|\le N\max_n|x[n]|=N\tag{1}$$

Note that for signals $x[n]=e^{j2\pi nl/N}$, $l\in\mathbb{Z}$, the bound is tight. However, for most other signals with $|x[n]|\le 1$ you will find that the maximum of $|X[k]|$ is substantially smaller than the bound given by $(1)$.


How does subcarrier correlation affect the BER of OFDM? I am studying the OFDM system and would like to examine it under different subcarrier correlation scenarios. I have calculated BER for the uncorrelated scenario and now I have rewritten it for the correlated channel.

If there is no noise, BER shoud be 0. It works for the uncorrelated scenario.  Why does it not work for the correlated scenario?

OFDM subcarriers are packed relatively tightly together. If you look at the original OFDM signal in the frequency domain, you may wonder why adjacent subcarriers are not interfering with each other. The answer is that subcarriers are orthogonal to each other. Even adjacent subcarriers, have 0 influence on each other, and are independent in that sense.

It can be proven that OFDM is the most efficient way to pack these subcarriers together in the sense that it is the smallest inter-subcarrier spacing whereby the subcarriers can all be orthogonal to each other.

So then, in the real world, during transmission, things may happen besides getting white noise added, that cause the subcarriers to lose orthogonality to each other, and become correlated at the receiver. Probably your textbook would explain what may cause the subcarriers to get correlated?

In that case, there would be inter-subcarrier interference (also known as inter-carrier interference), and that would give you some non-zero BER.


Why does fft energy leakage was solved if we take a sinus with Fd=100KHz with fs=1000KHz sampling rate We do an FFT with 256 members we get energy leakage.
Why choosing Fd=113.28Khz solves this problem?
Thanks

See my answer here on how to work with the units in question:
DFT exercise in the book Understanding digital signal processing 3 Ed
Then apply the fact that leakage does not occur when the DFT frame covers a whole number of cycles of your tone signal.

$$ Fs = 1000000 \frac{Samples}{Second}  $$
$$ Fd = 113280  \frac{Cycles}{Second} $$
$$ N = 256 \frac{Samples}{Frame} $$
$$ \frac{Fd}{Fs} \cdot N = \frac{\frac{Cycles}{Second}}{\frac{Samples}{Second}} \cdot \frac{Samples}{Frame} = \frac{Cycles}{Frame}  $$
$$ \frac{113280}{1000000} \cdot 256 =  28.99968 $$
Which is pretty darn close to a whole number.


Is it common to impose the sparsity on the Fourier coefficient itself? In compressive sensing, I see many works to impose the sparsity on the wavelet coefficients (e.g., by minimizing the L1 norm of such coefficients.) Another example in MRI is to impose sparsity on the gradient of the unknown image (total variation).
However, I'm interested in imposing the sparsity on the Fourier coefficients itself. I couldn't find for the moment. If there is any reference, could anyone let me know? Or, I'm not sure whether it doesn't make sense.

This question is typically the subject of a paper like Robust Uncertainty Principles: Exact Signal Reconstruction From Highly Incomplete Frequency Information, Emmanuel J. Candès, Justin Romberg, 2006:

This paper considers the model problem of reconstructing an object
from incomplete frequency samples. Consider a discrete-time signal
$f\in \mathbb{C}^N$ and a randomly chosen set of frequencies $\Omega$.
Is it possible to reconstruct from the partial knowledge of its
Fourier coefficients on the set $\Omega$?

From citations within, and the quantity of citing papers, there is a possibility to find many useful references. Typically, a potential issue with non-local Fourier décompositions is that discontinuities in signals/images induce slow decay in frequency coefficients. Wavelets are usually better at dealing with sharp variations.
However, some modalities of imaging or spectroscopy are indeed acquired in a Fourier domain, and in practice, the compressibility of data (decaying coefficients) can suffice (over strict sparsity).


finding inverse Z transform with usage of tables (LTI and causal sequences) problem is as follows

try to find the inverse Z

$Z^{-1} (\frac{3}{z+2}) = ???$

with the usage of z-transform tables

Ok, so in order to find something from the table, I thought that we expand with z and try to see if we  get closer to any kind of expression in the tables

$Z^{-1} (\frac{3z}{z*(z+2)}) = Z^{-1} (3*\frac{1}{z}\frac{z}{z+2}) $

then we get constant out

$=3*Z^{-1}(z^{-1}*\frac{z}{z+2})$

then we see that right side fraction can be found in the table, but the combined sum expression can also be found with further manipulation (as can be seen later)

So, I think we declare temp variable here to clean it up (this portion is found in the table green portion)

$F(z) = \frac{z}{z-(-2)} \iff f[k]=(-2)^k $

then we can make it so

$3*Z^{-1}(z^{-1}*F(z))$

then we find the direct expression in the table (we have delay of one unit because of $z^{-1}$ term, it is the orange portion)

$3*Z^{-1}(z^{-1}*F(z)) = 3* u[n-1]*f[n-1]$

so it looks like answer should be

$3* u[n-1]*f[n-1] = 3*u[n-1]*(-2)^{n-1}$

my teacher hasnt provided the solutions guide yet, despite that he promised to do so, and I reminded him about it, so I was just checking if this is the right solution and thinking...

I did the same problem, but I used polynomial long-division style to get the following series representation

$3*z^{-1}-6*z^{-2}+12*z^{-3}-24*z^{-4}...$

so based on that I guess that the $x[n]=3*(-1)^n*2^n $

where $ n>=0$

so, I was just a little bit confused does there actually exist the delay of one unit in the original input sequence (result of inverse Z), sometimes I get confused about the causal sequence when you multiply with the unit step sequence, but then, it seems that sometimes also it can just be stated that the original sequence must have been causal sequence, so you just define n>=0 or something like that???

Your solution is correct. This is the way I see it (assuming that the sequence is causal/right-sided):

$$\mathcal{Z}^{-1}\left\{\frac{3}{z+2}\right\}[n]=3\mathcal{Z}^{-1}\left\{z^{-1}\frac{z}{z+2}\right\}[n]=\\=3\mathcal{Z}^{-1}\left\{\frac{z}{z+2}\right\}[n-1]=3(-2)^{n-1}u[n-1]$$


Questions about autocorrelation, correlation, and randomness I'm learning about autocorrelation functions and trying to calculate some from a time series. Conceptually, what is the difference between the following two autocorrelation functions?

$$g_{1}(\tau)=\frac{\big\langle (A(t)-\left \langle A \right\rangle)(A(t+\tau)-\left \langle A \right\rangle ) \big\rangle}{\big\langle (A(t) - \left \langle A \right \rangle )^{2} \big \rangle}$$

$$g_{2}(\tau)=\frac{\big\langle A(t)A(t+\tau)\big\rangle}{\big\langle A(t)^{2} \big\rangle}$$

My understanding is that the autocorrelation for the second equation will asymptotically go to zero if there is no correlation between the events in the time series. Is that the case? What about for the first equation? And why?

For the first equation, is the average value $\left \langle A \right \rangle$ the average over all of $\tau$, or do you need to take a running average and calculate it for each time $t$ when doing the calculation?

What happens if the time series is correlated? Will it never approach zero (in the case of the second equation)?

Is there any a priori way of determining how far to take the lag time $\tau$?

Finally, how does all of this relate to randomness?

1) g1 assumes that there is a non-zero mean and that always needs to be subtracted off for the correlation calc to make sense.  g2 assumes that the mean is zero. g2 is incorrect if the mean is not zero.

2) the mean is taken over the whole sample because the autocorrelation calc has the underlying assumption that the process is stationary so the mean doesn't change during the process.

3) As the lag gets larger, then, if the process is stationary, then the autocorrelations should doe down approach zero. Note that there are "long memory" processes that don't have this characteristic but I don't think that they are
considered stationary processes. Someone feel free to correct me if I'm wrong about long memory processes not being considered stationary.


SIFT - why s+3 scales per octave? I have a problem with SIFT that I do not understand. Lowe [1] proposed in his work the s=3 levels of scale are enough for one octave. Afterwards, he mentioned that you need to compute s+3 levels.

Why there are 3 and not 2 additional levels required. I understand that you require one additional level above and one additional level below the scales since you search for extrema in neighbored scales.

For what is the third additional level of scale?

Thank you very much in advance!

[1] Distinctive Image Features from Scale-Invariant Keypoints
    D. G. Lowe
    Int. Journal of Computer Vision 60(2) (2004), pp. 91--110

We must produce s + 3 images in the stack of blurred images for each
  octave, so that ﬁnal extrema detection covers a complete octave.

For $s=3$ this means you will have $s + 3 = 6$ blurred images (the Gaussian images shown in the paper in Figure 1 on the left). Having $6$ Gaussian images will result in $5$ DoG images (shown in Figure 1 on the right). This will allow you to do the extrema detection on $s=3$ scales (using the method shown in Figure 2).