Convolution of one signal with an evenly spaced signal I do not specialize in signal processing so I wonder if there is any references to the following procedure.

Let $[n]=\{0,1,\ldots,n\}$.

Consider the function $f:[n] \to \mathbb{R}$ and $g:[kr] \to \mathbb{R}$, such that $\operatorname{supp} g = \{0,r,2r,3r,\ldots,kr\}$. Assume for simplicity, $n/r$ is a integer.

We are interested in computing the linear convolution $f*g$.

If $r$ is large, one way to do the computation is partition $f$ into $r$ functions $f_0,\ldots,f_{r-1}$, such that $f_i(t) = f(tr+i)$, and $f_i:[n/r]\to \mathbb{R}$.

We also define $\bar{g}:[k]\to \mathbb{R}$ as $\bar{g}(t) = g(rt)$.

Find $f_i*\bar{g}$ to each $i$, then we can read off $(f*g)(tr+i) = (f_i*\bar{g})(t)$.

This make sure we apply $r$ convolutions on sequence of length $n/r$ and $k$. Instead of $1$ convolution on sequence of length $n$ and $kr$.

This is known as polyphase decomposition. It is often used as en efficient implementation of filtering combined with decimation or interpolation.


Crash course on turning engineering formulae into computer code? I'm trying to teach myself about the basics of feature extraction from digital signals. Specifically, I want to try monophonic then polyphonic pitch recognition. I don't have any university level math or computing science experience (I'm self-taught), so I find the formulas in many texts and websites to be a road-block because I don't understand what most of the Greek symbols mean. E.g. I had to search and find that $ \Sigma $ is called Sigma and that it means "summation"...

As a self-taught recording engineer I have a decent understanding of what needs to happen from a high level:

Separate digital signal into component frequencies
Analyse the relative strength of each frequency bin
Evaluate the strongest frequencies and compare against known values to determine most likely chord

When I start trying to do any reading about these theories I can't translate it into code because I have no idea what the mathematical explanations are trying to tell me.

Maybe this is a question for a math forum, but I figured that a group of DSP enthusiasts and experts could help me fill in the gaps between my coding knowledge and my basic understanding of sound. I get that a Fast Fourier Transform takes a signal and transforms it from the time domain into the frequency domain, but I don't understand the process when I see it laid out in mathematical terms.

How does someone in my position move forward? What kind of resources exist to help someone with high-school algebra (but a "can-do" attitude!) learn enough to get comfortable with basic signal processing? Is there a less "math-intensive" way to gain an understanding of what needs to happen that could lead to a better understanding of the symbols and formulas commonly used?

I don't usually recommend this series, but Linear Algebra For Dummies might be right for you. Chapter 19, for example, covers all the Greek symbols.

You should also complement this with an introductory book on DSP, like Applied digital signal processing : theory and practice and DSP for audio An Introduction to Audio Content Analysis: Applications in Signal Processing and Music Informatics


Using recursive least square filter I am fairly new to the signal processing world and that being said I have little to no experience. The problem that I am having is that I am not quite sure how to use dsp.RLSFilter. So far I have only used highpass filter and it was straight forward - I had to just decide on the cutting frequency, type of highpass and sampling frequency, whereas for the RLS filter I have a ton of parameters to choose from. For example, how to decide on the method to calculate the filter coefficients? Furthermore, in the documentation it is stated:

  Call step to filter each channel of the input according to the properties of dsp.RLSFilter. The behavior of step is specific to each object in the toolbox.

and in the example step is not used:

rls1 = dsp.RLSFilter(11, 'ForgettingFactor', 0.98);
filt = dsp.FIRFilter('Numerator',fir1(10, .25)); % Unknown System
x = randn(1000,1);                       % input signal
d = filt(x) + 0.01*randn(1000,1); % desired signal
[y,e] = rls1(x, d);
w = rls1.Coefficients;
subplot(2,1,1), plot(1:1000, [d,y,e]);
title('System Identification of an FIR filter');
legend('Desired', 'Output', 'Error');
xlabel('time index'); ylabel('signal value');
subplot(2,1,2); stem([filt.Numerator; w].');
legend('Actual','Estimated');
xlabel('coefficient #'); ylabel('coefficient value');

So what is the difference between using step and using the method in the example? When I try using the same the same way as in the example I get the following error message: Array formation and parentheses-style indexing with objects of class 'dsp.RLSFilter' is not allowed.  Use objects of class 'dsp.RLSFilter' only as scalars or use a cell array. I tried using num2cell on x and d, however, I had 0 success.

Unlike a standard high pass filter where you set a cut-off frequency and other design parameters for a fixed filter result with a pass band ripple, stop band rejection, phase response etc.. the "recursive least squares filter" is an adaptive filter commonly used for channel equalization.

As a high level working example of a LMS equalizer, please see my post here where I derived and used the Wiener-Hopf equations to solve for the LMS filter coefficients. This would be an approach to equalization when the coefficients are solved for a block of data and then applied; most applicable to fixed filtering solutions such as compensating for distortions introduced in the RF portion of a receiver (compensating for analog filters for example).

Compensating Loudspeaker frequency response in an audio signal

To adaptively solve for the same coefficients, such as compensating for changing channel effects, two common algorithms are referred to as "LMS: Least-Mean Squared" and "RLS: Recursive Least Squares". I compare the two in my slide copied below, and for more details please refer to the following references

Rappaport, T.S. Wireless Communications Principles and Practice, Second Edition, Prentice Hall 2002

Proakis, J.  “Adaptive Equalization for TDMA Digital Mobile Radio,“ IEEE Transactions on Vehicular Technology, Vol. 40, No. 2 pp 333-341, May 1991

Sayed, Ali, Fundamentals of Adaptive Filtering, Wiley, 2003

Bingham, J.A.C., The Theory and Practice of Modem Design, 1988

The comparison of the implementation for the LMS and RLS adaptive equalizer is shown in the figures below, where vectors are indicated by single bars and matrices with double bars:

LMS: Easy!

RLS: Fast!

Finally to note that a least-squares equalizing filter is NOT the best choice for a channel with severe frequency selective fading as noise enhancement will result where deep fading nulls occur. For this condition, a decision feedback equalizer is often a better choice if time-domain equalization is to be used.


Do DCT each output coefficients have to be orthogonal? I have a question: Do DCT each output coefficients have to be orthogonal?
Here DCT is DCT-II type. It is used by the jpeg compression algorithm.

It doesn't have much meaning to ask if the coefficients are orthogonal. If you do a 1-d DCT transform then you get a 1-d vector as output. You need at least two vectors to talk about orthogonality.

The DCT transform however, is an orthogonal transform. This means that there exists an orthogonal matrix that transforms the input data into the dct coefficients.


Signal processing for audio and speech I have started learning about signals and I am interested in sound signals.

There are some questions that I need to resolve. For example in music, difference of notes like 'sol' and 'la' are about difference in frequency. But which features of signals reveal differences between piano and guitar. What sound features distinguish between two different syllables?

Can you give me a good learning source about audio processing?

In line with a previous similar question here are my suggestions:

There are so many nice books but I believe you should first have a look at the science of sound from Rossing for getting the most broad view on the subject.

Then you can look at the following books , each dealing with a particular dimension of the problem:

1- Acoustics_BERANEK

2- Elements of Acoustical Engineering_OLSON

3- Music, Physics and Engineering_OLSON

4- Art of Digital Audio_WATKINSON

5- Audio Engineering_SINCLAIR

6- Digital Audio Effects_ZÖLZER

7- Digital Audio Signal Processing_ZÖLZER

8- Fundamentals of Acoustic Signal Processing_TOHYAMA

9- Audio Signal Procesing and Coding_SPANIAS

10- Audio Programming Book_BOULANGER

11- PsychoAcoustics_ZWICKER

12- Acoustics and PsychoAcoustics_HOWARD

EDIT:
I'm sorry that I've forgotten to mention about books on speech processing, so here you can find a list:

1- Speech and Audio Signal Processing_GOLD

2- Speech and Hearing for Communication_FLETCHER

3- Discrete Time Processing of Speech Signals_PROAKIS

4- Fundamentals of Speech Recognition_RABINER

5- Introduction to Digital Speech Processing_RABINER


Equalisation of FFT spectra The image shows the frequency spectra of a signal recorded using a SDR. The SDR is of DirectDownSampling type and has some compensation in the signal chain to compensate for the DC spike which would appear at 370MHz the center frequency. But unfortunately this has caused a dip at the centre frequency. Now since the Tx channel is varying the signal power fluctuates which is shown in image 2 and image 3. The null at 370MHz is at about -22dB but this null should be at about -17dB with the other nulls. I want to be able to dynamically compensate for the 370MHz null during varying channel conditions bringing it up to -17db and also compensate for the slope on the flat tops adjacent to the 370MHz null. I have been trying of fit a line to the top of the flat tops with little success.

has some compensation in the signal chain to compensate for the DC spike which would appear at 370MHz the center frequency.

Yes, a typical problem of IQ demodulators is the LO leakage/DC spike. You really can't do much about that. You don't mention the SDR you're using, but many can be used with offset tuning, completely moving the LO out of your band of interest and thus completely solving that problem

  But unfortunately this has caused a dip at the centre frequency.

Why is that a problem? Looking at your spectra, I'd clearly say you've got four subcarriers working there, and the centre frequency happens to end in the gap between the center two. There's really not much you'd need to fix.

  Now since the Tx channel is varying the signal power fluctuates which is shown in image 2 and image 3. The null at 370MHz is at about -22dB but this null should be at about -17dB with the other nulls.

Why? It's a null. There's no data here. Simply leave it alone!

  I want to be able to dynamically compensate for the 370MHz null during varying channel conditions bringing it up to -17db

That's a bad idea. You'll break more than you fix when trying to equalize based on a cancelled component. You get overproportional noise amplification that way!

Instead, you'd want to base your equalization on something that comes through well.

This is exactly the discussion you probably read when you first were introduced to frequency-domain equalization: Zero-Forcing looks easy, but will amplify the noise on channel-suppressed frequencies. You'd be doing even worse: you'd be selectively taking the noise on a suppressed frequency and use that to equalize. That means you'll multiply a lot of noise onto your whole spectrum. Bad, bad idea.

  also compensate for the slope on the flat tops adjacent to the 370MHz null. I have been trying of fit a line to the top of the flat tops with little success.

Why should there even be a flat top here, or a straight line? This is clearly a multipath channel, and a curve fit isn't what you'd use to equalize that.

From the look of your signal, I'd say this is probably one larger or four low-subcarrier count OFDM transmissions; if that is the case:

The whole point of OFDM equalization is done by the OFDM receiver on each subcarrier individually, and is not something that you should even do before.

All in all, the frequency selectivity of this channel looks pretty benign.


FWGN model implementation vs explaination (textbook pages attached) I came across the FWGN model for multipath channel simulation in the book: "MIMO-OFDM Wireless Comuunications with MATLAB"

Here is an explanation about Clarke/Gans model:

Some code to generate a Chanel using this model is also in the book:
here is some of it:

function [h,Nfft,Nifft,doppler_coeff]=FWGN_model(fm,fs,N)
% FWGN (Clarke/Gan) Model
% Input: fm= Maximum Doppler frquency
% fs= Sampling frequency, N = Number of samples
% Output: h = Complex fading channel
Nfft = 2^max(3,nextpow2(2*fm/fs*N)); % Nfft=2^n
Nifft = ceil(Nfft*fs/(2*fm));
% Generate the independent complex Gaussian random process
GI = randn(1,Nfft); GQ = randn(1,Nfft);
% Take FFT of real signal in order to make hermitian symmetric
CGI = fft(GI); CGQ = fft(GQ);
% Nfft sample Doppler spectrum generation
doppler_coeff = Doppler_spectrum(fm,Nfft);
% Do the filtering of the Gaussian random variables here
f_CGI = CGI.*sqrt(doppler_coeff); f_CGQ = CGQ.*sqrt(doppler_coeff);  % <--- Why sqrt ?
% Adjust sample size to take IFFT by (Nifft-Nfft) sample zero-padding
Filtered_CGI=[f_CGI(1:Nfft/2) zeros(1,Nifft-Nfft) f_CGI(Nfft/2+1:Nfft)];
Filtered_CGQ=[f_CGQ(1:Nfft/2) zeros(1,Nifft-Nfft) f_CGQ(Nfft/2+1:Nfft)];
hI = ifft(Filtered_CGI); hQ= ifft(Filtered_CGQ);
% Take the magnitude squared of the I and Q components and add them
rayEnvelope = sqrt(abs(hI).^2 + abs(hQ).^2);
% Compute the root mean squared value and normalize the envelope
rayRMS = sqrt(mean(rayEnvelope(1:N).*rayEnvelope(1:N)));
h = complex(real(hI(1:N)),-real(hQ(1:N)))/rayRMS;

I have noticed that this implementation multiplies the noise with a square root of the filter response in the frequency domain, why did the writer use the square root of the doppler spectrum ?

A similar representation is also present here.

EDIT:
Another thing, the text doesn't mention the $-\pi/2$ phase shift nor does it implement it in the code. Why is that ?
Is that just multiplication by -j ? If so, why negative ?

I assume that the function Doppler_spectrum computes the desired power spectrum of the noise. Since the power spectrum of a filtered noise process gets multiplied with the squared magnitude of the filter's frequency response, the filter's frequency response needs to be square root of the desired power spectrum (assuming white noise at the filter's input).

In formulas, if a random process $X(t)$ with power spectrum $S_X(\omega)$ is filtered by an LTI system with frequency response $H(\omega)$, the power spectrum $S_Y(\omega)$ of the output signal $Y(t)$ is given by

$$S_Y(\omega)=S_X(\omega)|H(\omega)|^2$$


Difference between $\tt fft$ and $\tt dftmtx$ in MATLAB I have the following MATLAB code:
audioFile= audioread('myFile');
%% Reshape according to the given distortion period
distortionPeriod = 512;
audioFile= length(audioFile) / distortionPeriod;
audioFile= reshape(audioFile, [newLength, distortionPeriod]);

I noticed the following difference:
 audioFileFft = fft(audioFile); % incorrect transform
 audioFileDftMat = audioFile * dftmtx(distortionPeriod); % correct transform

>> norm(audioFileFft - audioFileDftMat )

ans =

   1.2985e+03

Why this difference exist? I would expect that those two operation(fft, dftmtx) would give the same result, as mentioned by MATLAB dftmtx's documentation.
How can I use the fft and still get the same results I got using the dftmtx?

Your resulting audioFile in the line audioFile= length(audioFile) / distortionPeriod; is a single number and not a vector signal.
And for this reason your audioFileFft and audioFileDftMat are computed from the processing of that number and not your initial signal vector in audioFile; which is not your intention. However if that line is for newLength, and this variable is an integer (i.e. initial audioFile length is multiple of distortionPeriod, here 512) then you can do this
audioFileFft = fft(audioFile, distortionPeriod, 2);
audioFileDftMat = audioFile*dftmtx(distortionPeriod);

And the resulting norm(audioFileFft - audioFileDftMat) should be very small.


How left out down-sampling in 3D-DWT just on the Z direction I am working on image processing by a  3D wavelet transform.

I have a problem with the classification size of the wavelet coefficients. As you know, when we apply 3D_DWT on the image (for example: Hyperion of the area by 256*256*128 size) the sub-band size would be decreased (for LLL in 1l level the size is 128*128*64).

How I could apply down sampling just for the Z direction (spectral direction)?

It means this size of the sub-band should be 256*256*62 after one level wavelet applying. I want to do it in the Matlab script (dwt3.m & wavedec3.m), but when I removed the down sampling, it was applied in all three directions.

You are applying a 3D DWT scheme whose down-sampling is coherent across the three  directions. Nothing forbids you from decomposing bands the way you want, and the literature is full of so-called isotropic, Tensor, Separable, Square or Standard implementations (vocabulary: separable/non-separable,
standard/non-standard, S-form/NS-form, rectangular/square, anisotropic/isotropic, tensor/Mallat, hyperbolic/isotropic, separated/combined wavelet transform)

To do that, it is more appropriate to perform 1D wavedec passes your own way on each direction separately.


frequency folding due to aliasing in FFT I am having problem understanding why negative frequencies are showing up to the right of positive frequencies in discrete complex FFT.

I am not sure how to apply aliasing concept to explain such phenomenon even after reading about modulo-2pi angular frequency concept at MIT OCW signal processing course.

Anyone ?

The Fourier transform $X(e^{j\Omega})$ of a sequence $x[n]$ is by definition periodic of period $2\pi$ in $\Omega$. The $N$-point discrete Fourier transform (DFT) of $x[n]$ is simply defined as a sampled version of the Fourier transform $X(e^{j\Omega})$, that is
$$X[k] = X(e^{j\frac{2k\pi}{N}}).$$
The periodicity of $2\pi$ in $\Omega$ for the Fourier transform becomes a periodicity of $N$ in $k$ for the DFT
$$X[k+N] = X[k].$$
The negative frequencies showing up to the right of the positive frequencies is simply a consequence of this periodicity. Or if you prefer, you get $X[k]$ for $k \in [0, N-1]$ instead of $k \in [-\frac{N-1}{2}, \frac{N-1}{2}]$. If you want the vector to contain first negative frequencies, then the DC and finally the positive frequencies instead, you can for example use MATLAB fftshift.

Those negative frequencies are not shown for the "real" DFT because if $x[n]$ is real, then $X[k] = X^*[-k]$ (conjugate symmetry). Hence, you only need $X[k]$ for $k \in [0, N/2]$.


Fourier Transform/Series DFT/DFS textbook problem (simple?) Suppose $x_c(t)$ is a periodic continuous time signal with period 1 ms and for which the Fourier series is

\begin{align*}
  x_c(t) &= \sum\limits_{k=-9}^9 a_k e^{j(2000 \pi k t)} \\
\end{align*}

The Fourier series coefficients $a_k$ are zero for $|k| > 9$. $x_c(t)$ is discretely sampled such that:

\begin{align*}
  x[n] &= x_c\left(\frac{n}{6000}\right) \\
  &= \sum\limits_{k=-9}^9 a_k e^{j(\pi k n/3)} \\
\end{align*}

$x[n]$ is periodic with $N=6$

Question: Find the DFS coefficients of $x[n]$ in terms of $a_k$.

My Work: The DFS coefficients of a periodic signal are:

\begin{align*}
  W_N &= e^{-j(2\pi/N)} \\
  X[k] &= \sum\limits_{n=0}^{N-1} x[n] W_N^{kn} \\
\end{align*}

Changing the variable $k$ in $x[n]$ to $m$ to avoid conflict and combining yields:

\begin{align*}
  X[k] &= \sum\limits_{n=0}^{N-1}
    \sum\limits_{m=-9}^9 a_m e^{j(\pi m n/3)}
 W_N^{kn} \\
  X[k] &= \sum\limits_{n=0}^{N-1}
  \sum\limits_{m=-9}^9 a_m e^{j(\pi n/3 (m - k))} \\
\end{align*}

I'm stumped on how to simplify or process this further.

I suspect this is the wrong approach. The problem gives Fourier series coefficients of the continuous function, there should be a direct way to convert them to the discrete Fourier series coefficients.

Textbook Answer: The answer given by the textbook is as follows. I am trying to figure out how to get to this answer.

\begin{align*}
  X[k] &= 2\pi \begin{cases}
    a_0 + a_6 + a_{-6} & k = 0 \\
    a_1 + a_7 + a_{-5} & k = 1 \\
    a_2 + a_8 + a_{-4} & k = 2 \\
    a_3 + a_9 + a_{-3} + a_{-9} & k = 3 \\
    a_4 + a_{-2} + a_{-8} & k = 4 \\
    a_5 + a_{-1} + a_{-7} & k = 5 \\
  \end{cases}
\end{align*}

Your result

$$x[n]=\sum_{k=-9}^9a_ke^{j2\pi kn/N}\tag{1}$$

with $N=6$ is correct. If you compare $(1)$ to the IDFT (or DFS)

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

you'll notice that in $(2)$ each term $e^{j2\pi kn/N}$ occurs only once for each value of $k$ (for fixed $n$), whereas in $(1)$, due to the periodicity of the complex exponential, each term $e^{j2\pi kn/N}$ occurs several times. E.g., for $k=0$ in $(2)$, in $(1)$ you get contributions for $k=0$, for $k=6$, and for $k=-6$ (because $N=6$). Consequently, comparing $(1)$ and $(2)$ gives

$$X[0]=N(a_0+a_6+a_{-6})$$

For all other values of $k$, the procedure is completely analogous: search for indices $k$ in $(1)$ which have the same complex exponential term $e^{j2\pi kn/N}$. You basically have to keep adding and subtracting the number $N=6$ to the indices, as long you remain in the range $|k|\le 9$, so e.g. for $X[2]$ you get contributions for $k=2$, $k=2+6=8$, and $k=2-6=-4$, because they all have the same exponential term $e^{j4\pi n/N}$.

PS: The term $2\pi$ in your textbook result doesn't make sense to me, it should be $N$, as already pointed out in a comment.


Detecting a single frequency in noisy data I am not very experienced in signal processing, and it seems like this should be a simple problem. I have two computers, one plays a tone a 19 kHz at a sampling rate of 44.1 kHz, and the other listens with a sampling rate of 44.1 kHz and tries to find if the 19kHz tone is present. The transmission of the signal is in a noisy environment, and I have tried to find the amplitude of the 19 kHz wave using an FFT as well as Goertzel, but neither has given me a conclusive value (the amplitudes are very small).

For the FFT, I am using JTransform, running the forwardFFT, then finding the bin in which 19 kHz lies and getting the amplitude by (imaginary * imaginary + real * real) for the components of that frequency bin.

For Goertzel, I followed the pseudo-code on this site: http://www.mstarlabs.com/dsp/goertzel/goertzel.html. I feel like both the FFT and Goertzel should give the same amplitude if the sample size is the same, but they give different values for that frequency bin.

I expect the amplitude to be greater than 10, but I do not believe this is working. I believe I need to at least filter out the noise below 19 kHz as a first step, but am unsure. Maybe some highpass filtering should be done first? Any suggestions on how to detect the 19 kHz sine wave? Thank you.
Updated:

Here is a plot of the raw data I am getting from the microphone:

I feel like this is just noise, but I am unsure.

If your goal is to find only the frequency then its an easy task .
Magnitude = sqrt(real*real + img*img)
This number will always be positive. Now , find the maximum value of magnitude in the array.

also,
                   frequency = max magnitude*(Fs/N)

Fs = sampling Frequency(44100 in your case) , N = size of FFT
Now you have the frequency. But the numerical value depends upon the fft size that you are using i.e 512 size will result in accuracy of 86Hz (+- 86) , 1024 => 43Hz and so on. So , the frequency will be 189__ or 190__ something :)


Matrix cross correlation in python I'm currently performing matrix cross correlation in python using :

C = scipy.signal.correlate2d(A,A)

where A is a 2D matrix, typically a picture. As you can imagine, it gets very long for large array, such that I'm looking for something faster.

Any hint ?
I also heard about scipy.signal.fftconvolve. Despite that it's much more rapid, it rarely give me the same result at the end.

scipy.signal.correlate2d(A,A) and scipy.signal.fftconvolve(A,A[::-1, ::-1]) will give the same result. The second one is much, much faster.


Can I ignore the Nyquist criterion? If I have a signal with frequencies up to 2000 Hz but I'm only interested in the frequency range up to 500 Hz, is it correct to sample at 1000 Hz $(f_\mathrm s \geq 2 \cdot f_\mathrm{max})$? Or do I have to apply an analog lowpass filter with cut-off frequency 500 Hz first?

I read this similar post, but didn't understand it well:
Ignoring Nyquist

You can only "ignore Nyquist" if you a-priori know that there is no higher frequency content (above your desired or required noise floor) that will alias down (or up!) to your band of interest.  But since your signal has content at and above half your sample rate, you will need to low pass filter that signal before sampling.

You will also need a sampling rate ABOVE 1000 Hz, if you are interested in a finite length of signal at 500 Hz.


Extraction of the phase lag form cross-correlation function I would like to find out the phase offset that my signal contains at the receiver side. For this, I use the cross-correlation function, where I cross-correlate signal before and after matched filter. The resulting cross-correlation function is represented in the figure below:

In the figure we can see that there is a pronounced peak at approximately 0.3 * 10^6 samples. I want to use this information to remove the phase error and my approach is following:

calculate the time lag
calculate the phase lag

this I code in the following way
time_lag = peak[0]/Fs # Fs - sampling frequency
phase = 2*np.pi*f_symb*time_lag # f_symb is symbol frequncy
print(phase)

out: 549160.4489439873

then I want to use this phase lag to shift the constellation, however even after the shift, the constellation is still tilted, which indicates phase error.
Can you point out, where I make a mistake?
Thanks!

The rotation of the constellation is due to carrier phase offset, which is not the same as a time delay. What the OP is seeing with the additional peaks is evidence that portions of the waveform itself repeat every 0.3x10^6 samples. Since the entire correlation is going down toward the edge of the plot, I assume the OP has done a linear cross-correlation and we are just seeing the effect from less overlap between the two waveforms that are being correlated.
Regardless of that, the carrier phase offset (constellation rotation) can be determined from the cross-correlation with zero time lag (the direct correlation function as a sum of complex conjugate products) as demonstrated below:
Consider such a correlation between the received $y[n]$ and known $x[n]$ signal where the only difference between the two signals is a phase rotation with
$x[n] = e^{j\theta[n] n}$
$y[n] = e^{j\theta[n] n}e^{j\phi}$
Where $\theta[n]$ is the known modulation of the carrier, while $\phi$ is an unknown phase offset (rotation).
The generalized correlation $r_{xy}[n]$ is a sum of complex conjugate products (the cross correlation function re-computes this at all delay offsets, but delay is NOT phase):
$$r_{xy}[n] = \sum_{n=0}^{N-1}x[n]y^*[n]$$
$$ = \sum_{n=0}^{N-1}e^{j\theta[n] n}e^{j\phi}Ae^{-j\theta[n] n}$$
$$ = Ne^{j\theta}$$
As we see, the phase of $r_{xy}[n]$ is the carrier phase offset we are trying to determine. Further, we have the advantage of processing gain in the correlation to optimally (in the presence of white noise) increase the SNR and therefore our estimate of the phase. The process is to use the cross-correlation function to align the received and transmitted waveforms, and once aligned (as evidenced by the largest peak) we can continue to use the generalized correlation function at that offset to determine the phase offset. Since the OP's cross-correlation function is maximum at offset = 0, the transmit and receive waveforms are already aligned, so the generalized correlation function can be used directly to determine phase (which would match the phase of that peak correlation at zero, since these would be the same).
I demonstrate this in this graphic pasted below showing correlation versus phase and frequency offset, with the added contributions of noise.


Maximum Likelihood Estimation in Presence of Colored Noise I am trying to test system identification in presence of measurement noise (1) A white Gaussian noise (2) Colored noise - pink, violet. When we are estimating parameters we do so in presence of iid, zero mean uncorrelated noise.

Q1: I would like to know if the colored noise are correlated or not. I think they have different distribution but I could not find any information if the samples will be correlated or not.

Q2: In estimation we assume that noise is additive white gaussian noise which is uncorrelated iid. What happens when the noise is not gaussian then how do we estimate theta?
For ex: $x = s(\theta) + Colored noise$
where we are trying to estimate $\theta$. Will the performance i.e, MSE vary with different levels of colored and non colored noise?

Samples of colored noise (taken at different times) generally
are correlated random variables because the
autocorrelation function of the noise process is not a delta function
as it is in the case of white noise. Thus, if we assume a zero-mean
process (noise is generally assumed to be regardless of its color),
then the covariance of two signals separated in time by $\tau$ seconds is
$R(\tau)$ where $R(t) = \mathcal F^{-1}(S(f)$ is the autocorrelation
function of the process (inverse
Fourier transform of power spectral density). Note that is is possible
for $R(t)$ to be zero for some values of $t$ (e.g. $R(t) = \operatorname{sinc}(t)$ is a valid autocorrelation function),
but it cannot be zero
for all nonzero $t$.

As far as the density function of any sample, if the process is
Gaussian, the sample is Gaussian even if the process has been filtered
with a linear filter before sampling. But if the process is not
Gaussian (it is,
let us say, LaPlacian), then while each sample will be LaPlacian, the
same cannot be said generally of samples of the process after filtering
of any kind. In other words, Gaussianity
survives linear filtering, LaPlacism generally does not.

So, how does maximum-likelihood estimation work when samples have
correlated noise? Consider the case when we wish to estimate the
unknown mean of a $\mathcal N(\mu, 1)$ random variable, and we
have two observations $x$ and $y$. In the standard case of
independent observations, the likelihood function is
$$L(\mu) = \frac{1}{2\pi}\exp\left(-\frac{1}{2}\left[(x-\mu)^2+(y-\mu)^2\right]\right).$$
The _maximum-likelihood estimator for $\mu$ is the number $\hat{\mu}$
that maximizes $L(\mu)$, which works out to be the number $\hat{\mu}$
that minimizes $(x-\mu)^2+(y-\mu)^2$. This is a quadratic in $\mu$
and the maximum-likelihood estimate turns out to be
$\hat{\mu}=\frac{x+y}{2}$. When the observations are correlated with
correlation coefficient $\rho$, then
$$L(\mu) = \frac{1}{2\pi\sqrt{1-\rho^2}}
\exp\left(-\frac{1}{2\sqrt{1-\rho^2}}
\left[(x-\mu)^2-2\rho(x-\mu)(y-\mu)+(y-\mu)^2\right]\right).$$
Once again we need to find the $\hat{\mu}$ where
$(x-\mu)^2-2\rho(x-\mu)(y-\mu)+(y-\mu)^2$ has a minimum.
We still have a quadratic in $\mu$ but now we get terms
like $xy$ in the coefficients. What $\hat{\mu}$ works out to
be is left for you to work out.

And what if we have $n$ observations where $n > 2$? All the above
still applies. For independent identically distributed Gaussian noise
in the samples, the sample mean
$n^{-1}\sum_i x_i$ is the maximum-likelihood estimate of $\mu$ but
in the case of correlated Gaussian random variables we get very
messy minimization problems because the quadratic that we are trying
to minimize depends on the inverse of the covariance matrix and
the result is a nonlinear function of the data instead of a simple
easy-to-remember result like the sample mean.

What if the noise is not Gaussian? The same principles apply -- set up the
likelihood function and find where it attains its maximum value --  but
the calculations are quite a bit different, all depending on
what you assume or know is the joint density of the observations.


elementary one and multidimensional real examples of continuous/discrete LTI systems By an LTI system, we mean a time-invariant linear map on continuous/discrete-time signals.
What (elementary one and multidimensional) real examples of continuous/discrete LTI systems do you suggest to provide a nice motivation for pure mathematicians?

If you care for multidimensional signals, let us use "shift-invariant" maps. To analyze such systems, Fourier transforms are tools (for engineers) of choice. I often say that Fourier transforms are to convolutions what logarithms are to multiplications: a tool that helps explaining what is less than obvious (diagonalization of the convolution operators), and that can become computationally efficient (FFT), and sometimes a game changer: mp3 and JPEG formats boosted the use of Internet.
Meanwhile, invariance is a very useful concept, and Fourier is a central operator in mathematics, and still subject to many research and open questions (for instance on series convergence). I had the opportunity to hire a skilled mathematics student (Jérôme Gauthier) for a PhD in signal processing on the synthesis in oversampled complex filter banks. I remember that during the overview, as a mere DSP engineer, I asked him about his knowledge on the Fourier transform, and he was (then) pretty confident. Later, he discovered another side of it (by the way, I rediscovered the least-squares parabolic extrapolation recently, which I though I mastered, from the computer science and real-time computing point of view).
So, the interplay between engineers and mathematicians is quite important, because one can shed a different light on how the other deals with it. For instance, research on the mathematics behind invariance is helpful to design novel algorithms. For instance, finding invariant transformations of graphs is of current interest.


how to remove color of light source? I have an image source emitting yellow light. Of course, when taking a picture with the light, the image has yellowness overall. How can I remove that color of light from the image? I have another picture of only background with the yellow light. Is there any way to subtract the image from the former one? Thanks

Do you want to remove the Yellow backgorund or make the Yellow light white?

If you want to color correct the image so that Yellow appears white then simply measure the R,G,B values of 'yellow', calculate the ratio between them. Then scale all the pixels in your image by that ratio.

Typically you would assume that green = brightness and just calculate red/green and blue/green then scale the red/blue pixels and just scale green to leave the overall brightness unchanged.

See http://en.wikipedia.org/wiki/Color_balance


Significance of outer product in structure tensor / second moment matrix How does one interpret the outer product $\mathbf{M} = \nabla I \otimes \nabla I = \begin{bmatrix}\left(\frac{\partial I}{\partial x}(x)\right)^2 & \frac{\partial I}{\partial x}(x) \cdot \frac{\partial I}{\partial y}(y)\\\frac{\partial I}{\partial y}(y) \cdot \frac{\partial I}{\partial x}(x) & \left(\frac{\partial I}{\partial y}(y)\right)^2\end{bmatrix}$ where $\otimes$ is the outer/tensor product and $I : \mathbb{R}^2 \to \mathbb{R}$ is an image? This matrix is called the structure tensor / second moment matrix.
Background:
The structure tensor tells us how much the pixels change between two image patches (infinitesimal small change $\mathbf{v}$).
\begin{align*}
    f(v_1, v_2) &= \sum_{x,y \in N} \left(I(x, y) - I(x + v_1, y + v_2)\right)^2\\
\end{align*}
where $N$ is a neighborhood. The first order Taylor approximation is
\begin{align*}
    I(x + v_1, y + v_2) &= I(x, y) + \mathbf{v}^T\nabla I(x, y)
\end{align*}
Then
\begin{align*}
    f(v_1, v_2) &\approx \sum_{x,y \in N} \left(\mathbf{v}^T\nabla I(x, y)\right)^2\\
              &= \sum_{x,y \in N} \left(v_1\frac{\partial f}{\partial x}(x) + v_2\frac{\partial f}{\partial y}(y)\right)^2\\
              &= \sum_{x,y \in N} \left(v_1\frac{\partial f}{\partial x}(x)\right)^2 + \left(v_2\frac{\partial f}{\partial y}(y)\right)^2 + 2v_1v_2\frac{\partial f}{\partial y}(y)\frac{\partial f}{\partial x}(x)\\
 &= \sum_{x,y \in N} \mathbf{v}^T\mathbf{M}\mathbf{v}
\end{align*}
I know this is an ellipsoid and there is a relationship to PCA and the covariance matrix. But I am trying to figure out the significance of computing the outer product of the same gradient. It looks a bit like the Hessian matrix $\mathbf{H}(I(\mathbf{x})) = \mathbf{J}(\nabla I(\mathbf{x})) = \begin{bmatrix}\frac{\partial^2 I}{\partial x^2} & \frac{\partial^2 I}{\partial x \partial y}\\\frac{\partial^2 I}{\partial y\partial x} & \frac{\partial^2 I}{\partial^2 y}\end{bmatrix}$ but without the 2nd derivative.

I have been looking again into this question. The tensor product doesn't really have such a special meaning here.
According to Nonlinear Structure Tensors:
"Although this tensor product contains no more information than the gradient itself, it has the advantage that it can be smoothed without cancellation effects in areas where gradients have opposite signs."
I think it is then just similar to a covariance matrix i.e. $M_{11}, M_{22}$ is the variance of the $x$ and $y$ direction. $M_{12}$ and $M_{21}$ is the covariance of $x,y$ and $y,x$. It is just important to see that we sum each component of the matrix (see Wikipedia). The outer product is then $E[X]E[X]^T$. So we consider an image $I : \Omega \to \mathbb{R}^n$ as a vector field and look for local variations $||dI||^2 = dI^TdI$. The outer product appears because all $4$ directions are compared.
Another explanation comes from diffusion equations (i.e. movement from a region of higher concentration to a region of lower concentration). Without the outer product, we would only look at the gradient at a point. The diffusion is the same in every direction (isotropic). This is called (nonlinear) isotropic diffusion:
$$\frac{\partial u}{\partial t} = \text{div}\left(g\left(|\nabla u|^2\right)\nabla u\right)$$
When $g(\cdot) = 1$ ("linear"), this is the regular heat equation because  $\text{div}\left(\nabla u\right) = \nabla^2 u$. The solution would be a Gauss filter, which blurs everything. Wikipedia calls nonlinear isotropic diffusion anisotropic, even though the scalar product $|\nabla u|^2$ is applied (instead of tensor product).
The second possibility is anisotropic diffusion. The PDE is
$$\frac{\partial u}{\partial t} = \text{div}\left(g\left(\nabla u\nabla u^T\right)\nabla u\right)$$
$\nabla u\nabla u^T$ is the structure/diffusion tensor. The diffusion is different in every direction. Diffusion across edges is prevented. When we only look at the diagonal of the matrix, the diffusion is isotropic again. Some more information can be found here and here.


Welch's Power Spectral Density - Time-Averaging Explanation I am just branching out to taking the power spectrum of short-term audio frames ($20$ ms) in order to extract useful audio features.  I have been reading about Welsh's method which states that after computing the squared magnitude of the result, the individual periodograms are time-averaged which reduced the variance of the individual power measurements.

I was hoping somebody could explain what and how time-averaging is achieved?  Is it done by looking at the individual squared bins within a given frame, or does it compare a previous frame to a present frame like a Spectral Flux?  Also, I take it that a periodogram refers to the periodicity of a given audio frame?  Thanks.

Periodogram averaging means that the squared DFTs of consecutive frames are averaged (per frequency bin). If $k$ denotes the frequency index and $l$ denotes the frame index, the averaged periodogram is computed as

$$\overline{P}[k]=\frac{1}{L}\sum_{l=0}^{L-1}P[k,l]\tag{1}$$

where $P[k,l]$ is the $k^{th}$ element of the squared magnitude of the DFT of the frame with index $l$. In (1) the averaging is done over $L$ frames. This is only useful if the signal can be assumed to be stationary during these $L$ time frames.


Check if a filter is IIR or FIR I have the following equation:

$$H(z) = \frac{z^3+5z^2+3z+1}{10z^3}$$

I want to find out whether if its IIR or FIR. Here is the steps I have gone so far:
$$H(z) = \frac{1}{10}\frac{z^3+5z^2+3z+1}{z^3}$$
Dividing nominator with denominator to find reminder gave me(I think its wrong):
$$z^3+5z^2+3z$$
After this step I have a $1$ as remainder and I do not know what to do with it.

Can someone please tell me the right procedure step by step?

In the general case you have

$$H(z)=\frac{P(z)}{Q(z)}$$

where $P(z)$ and $Q(z)$ are polynomials in $z$. If - as is the case in your example - $Q(z)$ just has one single term, $H(z)$ is definitely FIR, because you can simply divide each term of $P(z)$ by that respective power of $z$, and the number of terms of $H(z)$ equals the number of terms of $P(z)$.

If $Q(z)$ is a polynomial with finite zeros away from $z=0$, then $H(z)$ is generally IIR, except if all zeros of $Q(z)$ are cancelled by zeros of $P(z)$. One example of such a case is a recursive moving average (where $N$ is the window length):

$$H(z)=\frac{1}{N}\frac{1-z^{-N}}{1-z^{-1}}\tag{1}$$

The zero of the denominator is cancelled by one of the zeros of the numerator, so the transfer function $(1)$ can equialenty be written as

$$H(z)=\frac{1}{N}(1+z^{-1}+\ldots+z^{-N+1})\tag{2}$$

which is obviously FIR.


Detecting moved and dirty areas between two snapshots/images I'm new to image processing and trying to figure out best way to identify 'moved' and 'dirty' rectangles between two desktop screenshots.Basically, lets say there is a round circle moving on a screen, there is a text that is changing e.g. millisecond counter, and there is a image that is being updated. All this is happening at 60fps.
The algorithm should identify the moving circle as 'moved' region (and provide bounding rectangle of the moved object in source/initial and second/final image). Anything that changed in-place e.g. color/image change, text change etc should be identified as a 'dirty' region (with associated bounding rectangle). There could be a multiple of dirty & moved regions so the technique should be very fast on typical computer hardware (i5 processor) and be able to hit 60+ fps.
I'm open to using GPU/iGPU for speeding up the process, but at this point not sure about the right image processing algorithm.

I love abstract art.  It speaks to something inside of me that is outside the library, it speaks in words but not words.  One of my favorite artists is Wassily Kandinsky.

An image of his that could meet the requirements of the problem, both with "text" and translation of a circle, is his "Circle within a circle".  The bouncing globe animation is pristine, not textured, and doesn't have the text.

Also, the real world is messy and complex, so outside of explicitly specifying "this must work only on computer-generated content", any viable method has to hand a few real-world curve-balls.

Here is a copy of it:

So using MS Paint, I moved one of the circles and hashed up the lines in a way that should be consistent with local text changes.

Here it is:

So lets open up the files, convert to grayscale, and subtract them.

Here is some code:

import cv2
import numpy as np
from matplotlib import pyplot as plt
from PIL import Image

#read in images
#name of input file
fname1 = './wassily-kandinsky(A).png'
fname2 = './wassily-kandinsky(B).png'

#read the files, convert to grayscale
img1 = Image.open(fname1).convert("L")
img2 = Image.open(fname2).convert("L")

#convert to numpy array for math
arr1 = np.asarray(img1)
arr2 = np.asarray(img2)

#smooth (to remove some of the manual-editing real-world images artifacts)
kernel_width = 11
arr1b = cv2.bilateralFilter(arr1, kernel_width,75,75)
arr2b = cv2.bilateralFilter(arr2, kernel_width,75,75)

#subtract images
arr_diff = np.abs(arr2b-arr1b)

#show them
plt.imshow(arr1, cmap='gray')
plt.savefig("figure1.png", dpi=300)
plt.show()

plt.imshow(arr2, cmap='gray')
plt.savefig("figure2.png", dpi=300)
plt.show()

plt.imshow(arr_diff, cmap='gray')
plt.savefig("figure_diff.png", dpi=300)
plt.show()

And this gives slightly different initial images, and allows a decent difference image.

New (smoothed, grayscale) "unchanged" image:

New (smoothed, grayscale) "changed" image:

The difference between the two:

Here is a zoom in of the "signal area":

Next steps (to do when I get time):

clean up circle signal, show how edge histograms clearly give locations
block location from difference image, show locations of "salt-and-pepper" which would come from alternate text in same location.


How can I find the transfer function of the following block diagram? I've the following image and I want to find the transfer function from input $x(t)$ to output $y(t)$.
I know that I have to apply Laplace Transform, so the integrator becomes $\dfrac{1}{s}$, but I don't know what to do with the numbers $a$ and $b$. Should it still be the same?
If a and b still the same, then I found that $H(s) = \dfrac{Y(s)}{X(s)} = s^2-as-b$.

In such cases it can be useful to introduce auxiliary variables at the input of the integrators. For the given diagram you could use a signal $u(t)$ at the input of the second integrator. The equation for its Laplace transform $U(s)$ becomes
$$U(s)=\frac{1}{s}\big[X(s)-bY(s)\big]-aY(s)\tag{1}$$
You need another equation relating $U(s)$ to $Y(s)$, but that one is trivial. From those two equations you can express $Y(s)$ in terms of $X(s)$, and in this way you can obtain the transfer function $H(s)=Y(s)/X(s)$. The solution in your question is not correct. The correct solution must be a rational function, not a polynomial.


Multiplication property DTFT I was truing to solve an example of DTFT which is following multiplication property. The problem is
$$ a^n \sin(\omega_0 n) u[n]$$ we know that the definition of DTFT is
$$ X(j \omega) = \sum _ {n=-\infty} ^ {{+\infty}} x[n]e^{-j \omega n}$$
Multiplication in Time domain will be convolution in DTFT.

If we take the DTFT of $a^n u[n]$ we have $$\frac {1}{1-ae^{-j \omega}}$$ and DTFT of $\sin(\omega_0 n)\,u[n]$ will be  $$\frac{\pi}{j}\sum _ {l=-\infty} ^ {{+\infty}}  \delta(\omega + \omega_0 - 2\pi l) - \delta(\omega - \omega_0 - 2\pi l)$$

I have confusion how can I write it in the form of multiplication property.

Your sequence is

$$x[n]=a^n\sin(n\omega_0)u[n]\tag{1}$$

First of all, note that the Discrete-Time Fourier Transform (DTFT) of $(1)$ only exists if $|a|<1$. (The case $|a|=1$ can be handled by using Delta impulses). Anyway, for $|a|>1$ the DTFT of $(1)$ does not exist. Second, you can write $x[n]$ as the multiplication of $x_1[n]=a^nu[n]$ and $x_2[n]=\sin(n\omega_0)$. The DTFTs of these two sequences are given by

$$\begin{align}X_1(e^{j\omega})&=\frac{1}{1-ae^{-j\omega}},\quad |a|<1\\
X_2(e^{j\omega})&=\frac{\pi}{j}\left[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)\right],\quad -\pi <\omega <\pi\end{align}\tag{2}$$

Of course, $X_2(e^{j\omega})$ is also $2\pi$-periodic, but in $(2)$ only the interval $-\pi<\omega<\pi$ is considered (assuming $0<\omega_0<\pi$).
Note that in your expression for $X_2(e^{j\omega})$ you have a sign error. Furthermore, note that $X_2(e^{j\omega})$ is the DTFT of $x_2[n]=\sin(n\omega_0)$, and not of $x_2[n]=\sin(n\omega_0)u[n]$, as claimed in your question.

Now we have

$$X(e^{j\omega})=\frac{1}{2\pi}X_1(e^{j\omega})\star X_2(e^{j\omega})=\frac{1}{2\pi}\int_{-\pi}^{\pi}X_1(e^{j\theta})X_2(e^{j(\omega-\theta)})d\theta\tag{3}$$

The only thing you need to know to compute $(3)$ is the property

$$F(\omega)\star \delta(\omega-\omega_0)=F(\omega-\omega_0)\tag{4}$$

for any function $F(\omega)$. With $(4)$, evaluating $(3)$ with the DTFTs in $(2)$ results in

$$\begin{align}X(e^{j\omega})&=\frac{1}{2\pi}\frac{\pi}{j}\left[X_1\left(e^{j(\omega-\omega_0)}\right)-X_1\left(e^{j(\omega+\omega_0)}\right)\right]\\&=\frac{1}{2j}\left[\frac{1}{1-ae^{-j(\omega-\omega_0)}}-\frac{1}{1-ae^{-j(\omega+\omega_0)}}\right]\tag{5}\end{align}$$

The two terms in $(5)$ can be combined resulting in

$$X(e^{j\omega})=\frac{ae^{-j\omega}\sin(\omega_0)}{1-2a\cos(\omega_0)e^{-j\omega}+a^2e^{-2j\omega}}\tag{6}$$


Goodness of fit for complex valued curves (i.e. frequency responses in frequency domain) My apologies for perhaps the stupidity of this question.

Presume that one has the 'frequency response' $Y_{data}(k)$ of a system and also has an estimated model $Y_{syn}(k)$ that fits the data.

How can I quantify the 'goodness of fit' of $Y_{syn}(k)$ to $Y_{data}(k)$ ? Notice that both curves are complex-valued at every instant k (i.e. frequency domain).

For example, assume that I have two estimated models: $Y_{syn}^a(k)$ and $Y_{syn}^b(k)$. My goal (besides visualizing the residual) is to evaluate in a objective way the fit of both models to $Y_{data}(k)$ (the original data set). Is there a single metric to do this?

You could use the $p$-norm of the complex error $E(k)=Y_{data}(k)-Y_{syn}(k)$:

$$\epsilon_p=\sqrt[p]{\sum_{k=1}^K|E(k)|^p}\tag{1}$$

The most common values for $p$ are $p=2$ (Euclidean norm, $l_2$-norm) and $p=\infty$ (maximum norm or Chebyshev norm):

$$\epsilon_2=\sqrt{\sum_{k=1}^K|E(k)|^2}\tag{2}$$
$$\epsilon_{\infty}=\max_{1\le k\le K}\{|E(k)|\}\tag{3}$$

It eventually depends on your application which error measure best reflects the quality of the approximation. Note that for the maximum norm $(3)$, a single index $k$ determines the value of the error measure $\epsilon_{\infty}$, whereas the error measure $\epsilon_2$ in $(2)$ depends on all values of the complex error $E(k)$, and it is proportional to the square root of the mean squared error.


Will the resulting spectrogram from this algorithm be 1D or 2D? I'm implementing a spectrogram algorithm and I get the basic concepts, however, I'm struggling with something.

Let's say I have a 1D vector:

data = [0.36, 0.58, ...., 0.25]

I then compute a STFT of this data (with an overlap of say 128) this then produces a 2D array. I then compute the log of these values in order to compute the spectrogram.  Would the resulting vector be 1D or 2D? Since the STFT is 2D (as it is split into blocks), I'm guessing that the spectrogram will be a 1D representation since I'm calculating the logs.

Spectrograms are two-dimensional, usually calculated from a one dimensional input vector.

As an example say that you want a spectrogram with k-pts of frequency domain span and M-pts of time-domain span.   With NO overlap an input vector of length N = k x M will fill the spectrogram.  Every k samples a k-pt DFT is performed and makes one of the M time-axis vectors on the spectrogram.

For a running spectrogram with an overlap of k/2, then N/2 input samples will be needed to create a k x M spectrogram, and the time scale will be cut in half.

Usually the logarithm only comes in to scale the amplitude of the output of the DFT.  i.e., each of the k complex samples in the DFT output becomes a scalar value computed from the log of the complex sample magnitude.


Efficient frequency determination Is there any algorithm that tells you the most dominant frequency component, i.e. does argmax(fft(x)) or similar, that's more computationally efficient than argmax(fft(x))?

I need to continuously find which fft bin contains a tone, of (up to) 4096 bins in the full signal bandwidth. This is part of an RF detection system implemented in FPGA. As the input is taken from the air, the input signal is very noisy. Calculating argmax(fft(x)) with a length 4096 fft does what I need, but I was wondering if there's a way to do it without calculating all the 'unnecessary' fft outputs, to save fpga resources.

For example could there be some fft-type algorithm that bisects the spectrum to find the peak without calculating many of the outputs? Though I suspect this wouldn't really work without high SNR...

Is there any algorithm that tells you the most dominant frequency component, i.e. does argmax(fft(x)) or similar, that's more computationally efficient than argmax(fft(x))?

Mostly no.

The "trouble" is not the argmax, but fft. One alternative to this is to use autocorrelation and then determine whether that function is periodic and where does it crosses zero. The other alternative is to not evaluate the FFT in all 2048 possible points via the use of the Goertzel algorithm.

This covers why the answer is "mostly" no.

But here is a case where the answer might be yes.

Instead of trying to estimate the frequency via the typical "sinusoidal" concept, you can pass your signal through a sign function which will make it look like a square wave and then count its zero-crossings. (For complex $x$, you would have to do it twice). This is an inexpensive way of getting some estimate of "dominant" frequency.

Hope this helps.


Transform a data set by exploting the vectorfield I am somewhat new in the field of Digital Signal Processing / Image processing.
As shown in the figure, I have 4 straight lines $f_i(x)$ with $i = 1,\dots, 4$ that pass through $g(x)$.  Similiarly there are $4$ sample points $\mathbf{x}_i$ that lie on $g(x)$. The vectors $\mathbf{v}_i$ are the vectors obtained by joining the sample points $\mathbf{x}_i$ and the intersecting points between the lines $f_i(x)$ and $g(x)$.
The values of $\mathbf{x}_i$, $f_i(x)$ and $g(x)$ and hence the vectors $\mathbf{v}_i$ points are known.
If so, is it possible to find $g'(x)$ just by exploiting the information (direction and magnitude) contained in $\mathbf{v}_i$ such that the sum of the square of the magnitude of the vectors $\mathbf{v}_i$ is minimized?
My goal is to exploit the 1d-vectorfield and find a transformation matrix that maps $g(x)$ to $g'(x)$.
In the example above since $\mathbf{v}_1$ is pointing downward and $\mathbf{v}_4$ is pointing upwards (see $g(x)$), then intuitively, $g(x)$ must be translated about positive $x-$axis followed by some rotation. What is the best way to solve this problem?

The issue with considering anything as one dimensional in this problem is that they cannot be manipulated as required by the solution.
This can be solved by the optimsation of a function that modifies the $x_i$s towards the direction that minimises some expression of "error". The error here seems to be the distance of each point from the nearest intersection with the known lines.
The minimisation of that error alone is not enough to produce the required result here, but we will come back to this later.
The key "tool" here is the manipulation of the $x_i$s via modifications to the line they reside upon. This makes it a bit clearer why $g(x)$ is not a useful way of seeing things here.
For example, how are we going to rotate a $g(x)$? How are we going to "move it" towards the left or towards the right? If the answers to these questions are "we will just add some offset" or "send it through some rotation matrix", then $g(x)$ cannot stay as it is because it would be unaware of the space it is moving in.
Instead, suppose that all $x_i$s rest on a line running through and parallel to some vector $\vec{a}$. This, being a vector allows us to move it around (to translate the vector, simply add the required $x,y$ offset) and rotate it (to rotate the vector, simply scale its components by the appropriate coefficients). Consequently, we now cannot leave the description of the points as $x_i$ but instead something like $p_i$ taken from the set of points $P$ with each point being a pair of $(x,y)$ coordinates.
To complete the loop we need a way by which to change this line the points are resting on by some factor that is the result of the minimisation. We can do this by combining a translation and a rotation matrix and use the minimisation to inform us about the final $\text{offset}_x, \text{offset}_y, \theta$ parameters required to make the vector $\vec{a}$ coincide with the right points on the "fan" that is created by the $f_i$s.
In terms of practically doing this, you can establish a very simple function for manipulating the $p_i$s as $\text{mod_points}(P,\text{offset}_x, \text{offset}_y, \theta)=p_{i,n-1} \cdot M_{translate}(\text{offset}_x, \text{offset}_y,) \cdot M_{rotate}(\theta)$ where $P$ is the current position of the points, $p \in P$ and $n \in \mathbb{N}$ is the iteration step. Other platforms (e.g. MATLAB, Scilab, Octave, Julia and others) will have similar optimisation functionality of course.
At each iteration ($n$) use this to derive the error function you are trying to minimise. The error function takes the known lines and the positions of the transformed points and returns how far away from optimal is the solution. For example, $\text{fit_error}(f_i, \text{modified_P})$.
This $\text{fit_error}()$ will then be used in place of fun in scipy.optimise.minimise() and your initial $p_i$s (or $x_i$s as you call them now) will be the initial positions x0.
After a number of iterations, minimise() will have moved the line that the $P$s rest upon on the right intersection point with the "fan" that is defined by the (now known as) $x_i$s.
But, if you try to do this, you will notice that the minimisation of the error is not a sufficient condition because it allows the following solution:

Is this a valid solution according to the criterion of minimal error? Yes it is, this line produces points that rest on the known lines. Is it an acceptable solution?....No.
What else should we do to produce a valid and acceptable solution? Add a constraint to the calculation function so that the $P$ set gets closer to the intersections with the $f_i$ lines while at the same time the distances between the points themselves are also preserved. You can model this as two expressions for the error and add them with weight coefficients so that you can tune what sort of priority you want to give over either of the constraints.
Hope this helps.
EDIT:
You have two options regarding preserving the distance of the points:

Calculate the distance from the points, but infer where should the $\vec{a}$ go next. Then, since you do have the line that the points are on and their "inter-distances" are fixed, you can use the orientation of $\vec{a}$ to calculate where each point should be, knowing just the first one. (Subsequent points will be positioned at known distance $d$ from the previous known point towards the direction of $\vec{a}$).

Literally include all distances that need to be preserved in the expression for the error. In this case it would be distances $d(p_n,p_{n+1}), d(p_n,p_{n+2}), d(p_{n+1}, p_{n+2}), d(p_n,p_{n+3}), d(p_{n+1}, p_{n+3}), d(p_{n+2}, p_{n+3})$ where $d()$ is a function that returns the Euclidean distance and $n = 0$ is the index in $P$.


Algorithms for deformable image registration I need to perform a deformable image registration (DIR) for delineation of the different brain parts on the brain MRI 3D image. I faced with a problem that there are a lot of different approaches for solving this problem while I don't have enough experience to choose the right approach. First off I tried a lot of different algorithms from the itk library. After that I made several conclusions:

1) Affine registration is not appropriate for general case of this problem (I mean I always can find such examples that require not affine transformation for proper registration).

2) All itk cpu based algorithms for DIR takes at least several minutes to get the appropriate results. The only way to reduce this time is to use gpu based algorithms. (I use 256*256 and 512*512 T1 images)

3) Demons algorithm is really good approach for solving the general purpose registration, I just need to find the right parameters and number of iterations.

Are this conclusions right?
What other fast algorithms for DRI can you suggest? (I have found such algorithms as 6 types of Demon's, SG-LDDMM, IC-LDDMM, BSpline and fluid registration, but I can't find out do they really efficient for solving my problem)

P.S.
I would appreciate any links for gpu based source code for DRI.

Like Paul suggested, you may want to ask this question at the signal processing stackexchange as well. They may be able to provide you with more information. I am at best a talented amateur in the field.

Let's start with the biggest problem in image registration: there is no catch-all algorithm that just works for everything. Results are heavily dependent on image modality, site and also what you wish to use it for.

Luckily for you, brain mri is a fairly standard application. I would have a look at elastix, which is ITK based and used quite a bit at my department. It also has a nice database of parameter files for different applications, including several for brain mri. A quick glance reveals that most use B-spline registration for the brain.

I have mainly used registration for lung, where my personal experience was that Demons performed best. Note that Demons is NOT one single algorithm, although it is often presented as such. Rather it's a plethora of variations of the same basic idea.

Finally, some commercial software has registration algorithms that perform very well for specific applications. If you're associated with a hospital, ask researchers at your radiology or oncology department if they have access to something that works well.


Hilbert transform with a sliding window I have a input voltage signal which has disturbances present in it. I am actually discretizing the input signal into say $1300$ samples and taking Hilbert transform of the signal by taking $N$ samples at a time and then taking next $N$ samples at a time leaving behind the 1st sample..
Input Signal

While the output keeps oscillating instead of being steady for $N=50$ samples, (via Scilab) my time step is $0.4\ \rm{ms}$ and input signal frequency is $20\ \rm{ms}$ and that's why i am taking $50$ samples,
Hilbert transform output for $N=50$ samples.

If I increase or decrease the samples, the response has much more ripples than obtained here.

Cant i get output something similar to this??

 this is nothing but true rms of the signal obtained in excel.

My actual objective is to obtain the instantaneous amplitude and phase of the signal using Hilbert transform and I want to know what must be the length of the window (at present i have taken $N=50$ samples) and how to move that? Blocks or point by point?

The entire signal is not available and it just getting recorded then, and within the full signal getting recorded i need to obtain instantaneous frequency and phase.

Maybe including some code containing an minimum working example would help. I didn't see any problem.

I tried this because I am interested in related questions and recently learned about the Hilbert transform. I also added some noise because noise is always there, but often neglected in answers to questions like this (a pet peeve of mine), even though it sometimes completely defeats certain methods. Here's what I get:

import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

# Parameters
N = 1300
t = np.arange(0,N)
f = 1.0/20

# Discontinuous amplitude envelope
A = np.ones(t.size)
A[1*t.size/5:2*t.size/5] = .2

# Time signal with noise
y = A*np.sin(2*np.pi*f*t)
y += (np.random.rand(N)*2-1)*.2

# Filter noise
filt = True
if filt :
    cutoff = int(N*f*1.5)
    yrms1 = (y**2).mean()**.5
    yF = np.fft.fft(y)
    yF[cutoff:]=0               # Filter
    y = np.real(np.fft.ifft(yF))
    yrms2 = (y**2).mean()**.5   # Keep signal energy
    y *= yrms1/yrms2

# Calculate envelope
envelope = np.abs(sig.hilbert(y));

# Plots
plt.figure(1, figsize=[10,7], dpi=100)
plt.clf()
plt.plot(A, label='Actual envelope')
plt.plot(y, label='Time signal')
plt.plot(envelope, '.', label='Calculated envelope')
plt.legend(loc='lower right')
plt.show()
plt.draw()
plt.pause(.1)

EDIT:
The demonstrates performance in a sliding window. Just replace the line with the envelope calculation above with this code. Note that the envelope at each location is only ever calculated using N_win surrounding points.

envelope = np.ones(N)*np.nan
N_win = 50      # needs at least twice the period or so
lag = 20        # needs to be >= period.
for i in range(N-N_win) :
    envelope[N_win+i-lag] = np.abs(sig.hilbert(y[i:i+N_win]))[-lag]

Note, that for my implementation of this to work well, you need to be looking ahead of where you are calculating the envelope by about 1 period. For real time processing, that means your 'envelope' signal will have to lag behind your main signal by at least one period. This is because of edge effects of the Hilbert transform. I don't think there is really any way around this.


SDR and amplitude modulated signals Can I view every signal that has been transmitted using software defined radio as a amplitude modulated signal? For example I could make a frequency modulated baseband signal in GNU Radio and connect it to Adalm Pluto sink. Since analog part of Pluto SDR involves frequency mixing does it mean on top of FM I get AM?

Every signal transmitted by an SDR is a quadrature signal: the sum of two amplitude-modulated signals, with carriers that are in quadrature (one is at 90 degrees from the other).
A bit of trigonometry shows that this is equivalent to an amplitude and phase modulated signal; in other words, a quadrature signal is equivalent to $A(t)\cos(2\pi f_c t + \phi(t))$.
As long as you can express your signal with one of these two (equivalent) expressions, you can do anything! You can certainly do FM on top of AM if you want to. For many of these "exotic" signals, the challenge will be in figuring out a demodulation process, but don't let that stop you from experimenting.


Question about a diagram of PWM detection I have this diagram in a book(taken from Slideshare for clarity)
as described in the book PWM is put on the differential to obtain A(as starting of pulse) and B(as the end of pulse)
The A is used to start the integrator and B to stop it, at the same time A work as a clock for the sample and hold, to obtain the PAM signal of PWM
my question is, why A is used as the clock for the sample and hold?! I think B must be used for the clock to capture the amplitude! I don't understand this point!

From the diagrams, the integrator is reset to zero (sometimes this will be called "dumped", as in "integrate and dump") right after its value is sent to the hold.

The reason that I can see for doing this is that it appears that the leading edge of the PWM happens on a steady period, with the trailing edge varying to establish the width.  By sampling on the "A" pulse you insure that the timing of the PAM pulses is consistent; sampling on the "B" pulses would mean that the timing of the output would depend on the    amplitude; this would induce a nonlinearity.


(Rough) camera calibration of manipulated images Coming from a photogrammetry background, I am quite familiar with camera calibration, and the common assumptions that are used for crowd-sourced images (i.e. center point in the middle of the image, pixel resolution from the EXIF data etc.).

But I have been wondering a maybe naive question lately: Can I calibrate a manipulated image (e.g. downloaded from twitter or facebook), that has possibly been cropped or maybe even resampled and has no EXIF data? The nice thing is: More than 4 ground control points in the desired object coordinate system can be asummed available.

Is there any papers available on this problem?

1) It seems to me that you don't care about solving for correspondences or obtaining any high level paradigm in that sense. So I would treat this problem as a minimal problem and redirect you to check:

http://cmp.felk.cvut.cz/minimal/

Maybe specifically, 4-point method with unknown focal length:
http://cmp.felk.cvut.cz/minimal/p4pfr.php

You could find many variants in that page.

2) Well, this problem involves many issues of multiple view geometry theory. However, in this case, an uncalibrated approach would be more suitable.

You could obtain point correspondences from SIFT-like features. As your images are downloaded from social platforms, it is likely that they will contain a lot of such feature points. Then, you could match the descriptors in a RANSAC-like fashion and retain the correct matches. From there it is possible to compute $F$ the fundamental matrix. With further assumptions, constraints or knowledge about the scene you could compute the scale as well as  the 3D point coordintates. Some examples here:

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.31.4741&rep=rep1&type=pdf

Under certain constraints it is also possible to recover the camera matrices:

http://users.cecs.anu.edu.au/~hartley/Papers/eccv92/Higgins/higgins.pdf

After all, I would recommend you to read:

http://www.robots.ox.ac.uk/~vgg/hzbook/


FFT of triangular wave I am new in dsp.. so please bear with me.

I want to plot FFT of a triangular wave. I am using the code below to first generate the triangular wave and then take its FFT

Fs=50e6;
Ts=1/Fs;
NFFT=2^14;
Runtime=(NFFT-1)*Ts;
t=0:Ts:Runtime;

fin=5*1e6;

factor=Fs/NFFT;
bin=round(fin/factor);
fin=bin*factor;

fin_MHz=fin/1e6;

y_in=sawtooth(2*pi*fin*t,0.5);

% Calculate Spectrum
Y = fft(y_in,NFFT)/NFFT;
f = (Fs/1e6)*linspace(0,1,NFFT);

Ydb = 20*log10(abs(Y(1:NFFT))) + 6.02;

% % Plot amplitude spectrum.

plot(f,Ydb)
hold on
grid on
ylim([-150 5])
xlim([-5 55])
xlabel('Frequency (MHz)')
ylabel('dB')
set(gca,'YTick',-150:25:0)
refline(0,-78)

The response when input is 5 MHz is as follows :

The response when input is 5.5 MHz is as follows :

My questions are :

Do i have any mistakes in my code ?
The response should have been sinc^2 .. i dont see it
Why the difference between two plots.. i mean more spurs when ip is 5.5 mhz
The amp peak should be in 1 in fft .. it -1.825 dbfs .. ? why ?

Thankyou so much.

Do i have any mistakes in my code ?

Code seems clean. I would instead replace

% Calculate Spectrum
Y = fft(y_in,NFFT)/NFFT;
f = (Fs/1e6)*linspace(-0.5,0.5,NFFT); % it is better to center your spectrum

Ydb = 20*log10(abs(Y(1:NFFT))) + 6.02;

% % Plot amplitude spectrum.

plot(f,Ydb)
hold on
grid on

xlabel('Frequency (MHz)')
ylabel('dB')
set(gca,'YTick',-150:25:0)

  The response should have been $\text{sinc}^2$ .. i dont see it

It's not a sinc, per say. According to Fourier series, you can write a sawtooth wave as
$$x(t) = - \frac{A}{\pi}\sum_{n=1}^{\infty} \frac{1}{n} \sin n \omega t$$
To verify plot

plot(real(Y(1:NFFT)))

  Why the difference between two plots.. i mean more spurs when ip is 5.5 mhz

Note that the spectral information is in $[-f_{in},f_{in}]$. The FFT will display you information in $[-\frac{F_s}{2},\frac{F_s}{2}]$ in discritized bins. Therefore, anything outside the band of interest [-5,5] MHz serves as a replica.

  The amp peak should be in 1 in fft .. it -1.825 dbfs .. ? why ?

I don't quite get what you mean here.


Discrete and Continuous Signals I am new to the signals area and have been reading through a lot and have some questions.

Suppose I have a saved audio file on my computer, .wav file. I can view the time domain of the signal by identifying the sampling rate and bit depth. The resulted domain signal is a discrete signal, since everything stored in a computer is numbers.Now, I can convert this discrete signal into the frequency domain using the Fourier Transform algorithm. Since it is discrete, I must be using DFT equation below

Now, I can view the spectrogram of either continuous or discrete using the STFT, which computes the Fourier Transform over segmented windows of the whole signal (where we use FFT for DFT to reduce the complexity). The STFT equations of the continuous and discrete signals are:

My questions are:

1- Is my understanding correct?

2- Can we apply Hann Windowing on discrete signals?

Your understanding is correct and you can apply Hann windowing on discrete signals. The choice of your window depends on how easily you can visually resolve two sinusoids that are closely spaced in the spectrogram. $$STFT = X(m,\omega) = \sum_{n=-\infty}^{n=+\infty}x_nw_{n-mL}e^{-j\omega n}, $$ where $L$ is the step size of spectrogram (interval after which you place next window). and $W$ is the size of window, which decides resolution of frequency (row-wise separation of frequencies). Higher the $W$, more finer the resolution. You can see there is overlap between adjacent windows, the overlap being $W-L$.


Relationship between the carrier frequency and the propagation distance? For electromagnetic waves, Is there a relationship between the carrier frequency and  the propagation distance?

When transmission distance increases the signal tends to lose power, so carrier signal is added along with the message signal ( to strengthen the original message signal). Is that true?

Even in vacuum, there exists an attenuation that is inversely proportional to the frequency squared of the signal. It is called the free-space loss.

The ratio of received power vs. transmitted power (for polarization-matched antennas aligned for maximum directional radiation and reception) reduces to:

$$
P_r / P_t = (c/(4\pi ~d ~f))^2 ~G_{0t}~G_{0r}         ~~~~~~~~~~~~~[1]
$$

where $c$ is the speed of light in vacuum ($\approx3\times10^8 ~$m/s), $d$ is the separation between antennas (in meters), $f$ is the central frequency (in Hz) of the signal, and $G_{0t}$, $G_{0r}$ are the antenna gains under the conditions stated above.

As can be seen from Eq. [1], for a specific power ratio (attenuation), there is a fixed relation between the frequency and the propagation distance.

Equation [1] is a particular case of the well known Friis transmission equation. The term $(c/(4\pi ~d ~f))^2$ is called the free-space loss factor, and it takes into account the losses due to the spherical spreading of the energy by the antenna.

See, for example, the following reference for further details:

C. A. Balanis, "Antenna Theory. Analysis and design", 2nd edition, 1997. (Section 2.17, Chapter 2)

The inverse relationship between distance and frequency is not due to the propagation phenomenon itself, but due to the physical dimensions of the receiver antenna (which does depend on frequency, in something called the effective area of the antenna). In practice we always have an antenna matched to the desired frequency, and thus we have a propagation distance limited by this restriction.

Regarding the second part of your question, normally a carrier signal is added (for example, in AM modulation) for simpler detection circuits. Besides, the carrier signal is used sometimes to perform synchronization at the receiver. I have not seen before the use of a carrier to "strengthen" the message signal.

Actually, instead of using an extra power to transmit a carrier, it can be used to increase the power of the message signal, and thus let the message travel a longer distance with the same attenuation.


Interpretation of Histogram in Statistical Image Processing I am learning statistical image processing by myself. In papers and books, it always show the histogram of original images and gradients as the following image shows. The histograms of images vary significantly while histograms of image gradients show some similarity. Does it assume that each pixel in images obey the same probability distribution for the histograms of images? Does the histogram of any image gradient obey the same probability distribution?

In the paper Image Denoising Using Scale Mixtures of Gaussians in the Wavelet Domain by Javier Portilla, Vasily Strela, Martin J. Wainwright, and Eero P. Simoncelli there is one paragraph

Contemporary models of image statistics are rooted in the
television engineering of the 1950s (see [4] for review), which
relied on a characterization of the autocovariance function for
purposes of optimal signal representation and transmission. This
work, and nearly all work since, assumes that image statistics are
spatially homogeneous (i.e., strict-sense stationary). Another
common assumption in image modeling is that the statistics
are invariant, when suitably normalized, to changes in spatial
scale. The translation- and scale-invariance assumptions, coupled with an assumption of Gaussianity, provides the baseline
model found throughout the engineering literature: images are
samples of a Gaussian random field, with variance falling as
in the frequency domain. In the context of denoising, if one
assumes the noise is additive and independent of the signal, and
is also a Gaussian sample, then the optimal estimator is linear.

image statistics are spatially homogeneous What does it mean? Does image statistics means the histogram?

an assumption of Gaussianity What is Gaussian?

images are samples of a Gaussian random field If one image is considered as a random field, can histograms be used? The assumption that each pixel obeys the same probability distribution will not hold.

Does it assume that each pixel in images obey the same probability distribution for the histograms of images?

Images of different scenes will definitely not obey the same probability distribution of the pixel values.

Histograms are one way that people use to do dimensionality reduction: move from a 2D image to a 1D signal.

  Does the histogram of any image gradient obey the same probability distribution?

What you are seeing in the image gradient is the "diffs" in the image. Because images are generally low-pass in nature,  this means you are picking out the places where they change.  There will be (at least) two components to this change: how the scene being imaged changes and how the sensor capturing the image perturbs the "true" pixel values.

For the same camera taking the images, this second component should be very similar across all images.

  image statistics are spatially homogeneous  What does it mean? Does image statistics means the histogram?

means that the image statistics are very similar regardless of where in the image you look. One way the statistic show up would be in the histogram, yes.

  an assumption of Gaussianity What is Gaussian?

Gaussian means that the noise (random fluctuations in the image) follows a normal or Gaussian distribution.

  images are samples of a Gaussian random field If one image is considered as a random field, can histograms be used? The assumption that each pixel obeys the same probability distribution will not hold.

If the images are not random, then they will follow some well-defined deterministic rule.

Certainly, histograms can be used.


Combining CDMA with differential Manchester-like encoding? I wonder to which extend it is possible to combine a self-clocking code (like Differential Manchester encoding) with synchronous CDMA, and how it would be done?

Suppose I want to receive some bits from devices where each of them has sufficiently synchronized clocks but the receiver of the bits does not (and hence suffers from clock drift, delays etc).

Right now (without CDMA) I do not even need to care about the clock ... I just detect the transitions at the receiver and can infer the bits (because there is a guaranteed transition every bit). I wonder if something similar can be done while still using CDMA?

Take for example 3 IoT devices which send data in parallel via 4-order Walsh-Haddamard code.

The code matrix is:

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

and the three devices send bits $001011010101$, $001010110011$ and $001010001111$. The waveform as seen at the receiver looks like:

(as per assumption above they sum up nicely in the channel and preserve orthogonality of the data bits but the receiver does not know the clock). We can see that there is not necessarily a transition (or at least change of sample value) every bit/symbol. I tried replacing each 1 with $[1, 0]$ and each 0 with $[0, 1]$ in the code matrix but this doesn't change anything significantly.

What would be the best way to implement this?

(I hope my question is clear enough. If not, please leave a comment and I will clarify)

With CDMA everything you would traditionally do with a symbol such as timing and carrier recover still applies- but in CDMA the symbol is the series of many “chips”.

So yes you could certainly add any higher level coding to the symbols where in CDMA the symbol may represent a data bit such as +1 and -1 only (perhaps), but each of these data bits after such coding would then be expanded to the CDMA sequence.

There is no motivation to change the sequence itself other than its cross correlation properties to the other codes used and the overall processing gain required. Timing and carrier recovery is done after despreading.


What is homogeneous noise field? What is the formal definition of homogeneous noise field? Is there any example of it?

Thanks.

Think of being at the center of sphere with a very large radius that has a very large number of uniformly distributed noise sources on its interior surface. What you receive at the center ls spherically isotropic noise.  One can envision a similar cylinder for cylindrical isotopic noise.

These are ambient noise models.  The interesting thing is that the noise is spatially correlated between two points and has a sin(x)/x form for spherically isotropic noise.

This is a model typically used in SONAR where the electronic thermal noise of each sensor is lower than the ambient noise.

This model is derived in a number of books or at least referenced in just about every SONAR or Underwater Acoustics book.  I like the presentation in  Underwater Acoustic System Analysis by William S. Burdic (first edition) because he also derives the half sphere and infinite plain models. The only criticism of Burdic's analysis is that it assumes straight line propagation. Van Trees, Detection Estimation and Modulation Theory, Vol 4, Array processing has a section on ambient noise fields with many references.

The really interesting thing of spherically isotropic noise is that for a given temporal frequency, the zeros of the spatial correlation function are separations of $\lambda/2$ which is also the spatial Nyquist criteria.


How is a Particle Filter used to Estimate Parameters of a State Transition Function? In the tutorials on particle filters that I have seen, it seems that the state transition function is already known. This example from Matlab, for example, states the the Plant or State Transition function is provided ahead of time.
This video, also from Matlab, talks about predicting future state to move particles forward based on the ded reckoning  done by a robot. In that case, the robot knows what it's motion vector is (with process noise) at each step, but it is solving for its position and orientation. So the particle filter isn't involved in predicting this delta that results in future state.
I noticed the same thing in Matlab's description of the Unscented Kalman filter, where the state transition function is given a priori, when the filter is instantiated.
I'm looking for an explanation on how to use the particle filter (and the UKF) to estimate the parameters of the state transition function, but I'm not finding anything about that.
I suppose I can treat the parameters of the state transition function as the space in which I want to generate the particles, and transate those into the same space as the measurements. But for a simpel case of constant velocity in 3 dimensions, that would mean generating particles in 6 dimensions.  Is this how it's done with a particle filter? If not, what am I missing and where can I get a simple explanation of how this would be done?

I suppose I can treat the parameters of the state transition function as the space in which I want to generate the particles ... Is this how it's done with a particle filter?

More or less, yes.  You'll find the richest sources of information on this if you look for books on system identification and on adaptive control (of which one important branch is doing system identification on-line, then adapting one's controller to match the identified system).
I can't recommend a good system identification book, but Karl Åström and Björn Wittenmark's "Adaptive Control" is a pretty good book on that topic.


Denoise Image with Gaussian Noise Using MATLAB / Octave I want to remove a noise for an image using MATLAB, when the observed image is
$$f=u+v$$
where $u$ is the restored image (is the image i want recovered) and $v$ is the gaussian noise.

To restore $u$, I solve the following minimization problem:
$$\min_{u \in H^1(\Omega)} \int_\Omega \gamma|\nabla u(x)|^2dx+  \int_\Omega (f(x)-u(x))^2dx,$$
where $\gamma$ is the regularization coefficient and $\Omega=[0,n]\times[0,m]$ and $[n,m]=$size($u$).

I want to solve the PDE (Euler-Lagrange) using MATLAB:
\begin{eqnarray}
div(\gamma \nabla u) + u = f \;in \;\Omega \\
\frac{\partial u}{\partial n}=0 \; in \; \partial \Omega
\end{eqnarray}

Can anyone help me to solve this problem? Thank you!

I tried the following code :

clear all, close all,clc

uor=imread('gourd.bmp'); % the original image
u0 = imnoise(uor,'gaussian',0,0.01);
u0=double(u0);
[m n]=size(u0);
uor=double(uor);
u=u0;
c=0.028;
h=1;
for Iter=1:50,
    for i=2:m-1,
      for j=2:n-1,
          Lap=0.003*(u(i+1,j)+u(i-1,j)-4*u(i,j)+u(i,j+1)+u(i,j-1));
          u(i,j)=(u0(i,j)+(1/(2*c*h*h))*Lap);
     end
    end
    for i=2:m-1,
          u(i,1)=u(i,2);
          u(i,n)=u(i,n-1);
        end

    for j=2:n-1,
          u(1,j)=u(2,j);
          u(m,j)=u(m-1,j);
        end

        u(1,1)=u(2,2);
        u(1,n)=u(2,n-1);
        u(m,1)=u(m-1,2);
        u(m,n)=u(m-1,n-1);

en=0.0;
    for i=2:m-1,
      for j=2:n-1,
      ux=(u(i+1,j)-u(i,j))/h;
      uy=(u(i,j+1)-u(i,j))/h;
      fidelity=(u0(i,j)-u(i,j))*(u0(i,j)-u(i,j));

      en=en+c*fidelity;
      end
    end

Energy(Iter)=en;

%  Error between uor and u0
 ur=reshape(u,m*n,1);
 uori=reshape(uor,m*n,1);
 residu=norm(ur-uori)/norm(uori);

 [peaksnr, snr] = psnr(uor, u);

disp(['    iter ' num2str(Iter), ' :     Error = ' num2str(residu), ...
    ' ,    Peak-snr ' num2str(-peaksnr), ' ,    SNR ' num2str(snr)]);

  end

% show the structural similarity index for measuring image quality
[ssimval, mapssim] = ssim(u,uor);
disp([' the structural similarity index is ' num2str(ssimval)]);
figure,imshow(mapssim,[]); axis square;

figure,imagesc(u); axis image; axis off; colormap(gray);

The original image is here :
https://www.dropbox.com/s/4bccby1f4lxp4j9/gourd.rar?dl=0

Best regards

You're trying to solve what's called Perona Malik Non Linear Diffusion Problem (Sometimes people call it, by mistake, Anisotropic Diffusion).

Anyhow, the simplest code for that is Anisotropic Diffusion (Perona & Malik) on The MATLAB File Exchange.

There is a more advanced (Anisotropic for real) algorithm in Fast Anisotropic Curvature Preserving Smoothing (Also on the File Exchange - Fast Anisotropic Curvature Preserving Smoothing).

Enjoy...


Can I use standard computer vision techniques for images taken in the NIR spectral range? I'm relatively new in Computer Vision and I need some advice with respect to the general possibilities and shortcomings of images taken in the near infrared (NIR) range. I passed a great deal of time researching about that matter but I rather found specific papers in areas like pedestrian detection and face recognition which not completely answered my question.

The background is that I have to chose camera and lens for a stereo-vision system which will be used for object detection tasks.To make the system robust against variations in illumination I thought about the feasibility of employing active illumination with NIR-light and a NIR-bandpass filter. Thus, working exclusively in the NIR spectral range.

I found publications about stereo matching in the far-infrared range, which come to the conclusion that it is not possible to obtain a dense depth field. However, I could not find anything related to the NIR range. It may be because it is too obvious but unfortunately not to me at this point.

Thats why I would like to know in particular:

if it is possible to employ standard stereo matching algorithms like block-matching to images taken in the NIR spectrum and if so would it be possible (of course depending on the conditions) to obtain a dense depth map.
if standard feature detection and description techniques like SIFT, HOG,etc produces similar results in the NIR spectrum.
if a neural network which was trained with visible-light images would be also able to detect an object in an NIR image. (I am aware that this is somehow connected to the prior question)

I would be grateful for any advice, personal experience, book/paper references you may could share with me.

Thank you!

As far as I know, it is a common practice to use regular image processing techniques in this field. In fact I have recently finished a project in this field.

  if it is possible to employ standard stereo matching algorithms like
  block-matching to images taken in the NIR spectrum and if so would it
  be possible (of course depending on the conditions) to obtain a dense
  depth map.

Yes, it is, I think, since NIR images as opposed to Far-Infra red Images has a lot of texture which you probability can (It is very very similar to RED band of the image), as an example people doing these:

repository.upenn.edu/cgi/viewcontent.cgi?article=1115&context=cis_reports

  if standard feature detection and description techniques like SIFT, HOG,etc produces similar results in the NIR spectrum.

I think so, since again it has lots of texture and wherever there is texture is corners, edges and all types of  feature points,which you can extract for further porcessing, check this out:

enter link description here

I think you are better to also take a look at reference books such as:

Machine Vision and Advanced Image Processing in Remote Sensing:

https://books.google.com/books?isbn=3642601057

Remote Sensing, Models, and Methods for Image Processing:

https://books.google.com/books?isbn=0126289816

Hope you've got your answer.


Digital 16-bit PCM audio signal scaling and STFT spectogram amplitude I'm sampling a 16-bit PCM audio signal at 44100 Hz and I'm trying to understand the amplitude level of both the digital audio signal and its corresponding spectrum. The final goal is creating audio fingerprints for visualization and sound comparison. By now, I'm conscious I'm not calibrating my hardware so that dB scale is just a relative measurement to itself

First, for the 16-bit PCM audio signal I do understand the signal (Audio ADC lectures) ranges from -32768 to 32767, scaled to [-1 to 1) sig_scaled = sig_int16 / 32768 so that both -32768 and 32767 correspond to the maximum amplitude in db (0 db):

max_db = 20*np.log10(abs(+-1))= 0 dB

*Here, I do also understand it could be scaled to [0, 1) sig_scaled=(sig_int16-+32768)/65536

Meanwhile, the minimum one (~96db) corresponds to:

min_db = 20*np.log10(1/65536) ~= -96 dB

Which scaling is better [-1 to 1) or [0 to 1)?
How to avoid the "RuntimeWarning: divide by zero encountered in log10" when the ADC lecture or sig_int16 value is 0?

I think [-1,1) makes more sense since a sine wave goes from -1 to 1 but the obtained spectograms for both scales are similar.

Second, regarding the spectrum I've performed an stft on the [-1, 1) scaling, obtaining the following image (stft_db = 20*np.log10(stft/np.max(stft))):

What's the meaning of the -159.93 dB minimum and 0 dB maximum value? / Shouldn't the minimum value be -96 dB or bigger?
Which scaling ([-1, 1), [0, 1) or dB amplitud) should the STFT be performed on?

Thank you folks, let me know if you need more explanations or code!

1: Scaling

Audio corresponds to a pressure wave, so it is mostly symmetrical in amplitude. Therefore audio is normally represented as a signed value (with 0 average). [-1,1) is the standard representation.

You may offset your signal by 1/2 and represent the wave in the range [0,1), but this adds nothing but a huge component at 0 Hz in your spectrogram.

2: Warining when log(0)

You cannot calculate log(0). If that bothers you, create a function that produces a minimum value. Something like mylog(x) = log10(max(abs(x),1E-5)).

3: Small values in spectrogram

Values in spectrogram indicate energy per frequency bin. The more bins, the smaller values you may get there. And you may well have a 0 as result (see (2) above).

4: Signs again

Stick to signed values. See above.


Mathematical Approach to Detect If a 2D Signal Is Separable In the Special 2-D Sequences page the following examples are demonstrated,

2D dirac
2D diagonals
2D unit step function

Is there a more defined method or series of steps for determining if a function is separable or not? Other than:

$$x(n_1, n_2) = f(n_1)g(n_2) $$

Note: is there any difference between this and saying that:

$x(n_1, n_2) = f(n_1)\times g(n_2)$, where $\times$ is the orthogonal product (See SIMG-716 Linear Imaging Mathematics I section 2.1)? Is this simply referring to the inner product?

To me the solutions are highly intuitive and the requirement to solve such problems solely relies on the understanding of properties associated with the function.

Let's assume our data is in finite dimension.
So $ x \left[ m, n \right] \in \mathbb{R}^{M \times N} $. So it can be written as a matrix $ X \in \mathbb{R}^{M \times N} $.

Using the SVD Decomposition the matrix can be written as:

$$ X = \sum_{i = 1}^{n} {\sigma}_{i} {u}_{i} {v}_{i}^{T} $$

Seprable 2D Matrix is defined as a rank 1 Matrix which can be composed by Outer Product of 2 vectors.
Looking at the SVD Decomposition of $ X $ we can conclude that $ X $ is separable operator if and only if $ \forall i > 1 \; {\sigma}_{i} = 0 $ and it is given by:

$$ X = {\sigma}_{1} {u}_{1} {v}_{1}^{T} $$

So, here is a simple recipe:

Write you data in Matrix Form.
Apply the SVD on the matrix.
In MATLAB it can be done using svd().
Look for the rank of the matrix. Namely how many non vanishing singular values it has.
If the number is 1 the matrix is considered separable. It can be separable by the its first left eigen vector and first right eigen vector.

For example on Separable Operators have a look at - How to Prove a 2D Filter Is Separable?.


Image Processing- ornament detection So I am working on a project where we are supposed to detect and count the red and gold ornaments in the picture provided in the link below. I am having issues detecting and counting each ornament and cant seem to figure out why. I will outline my process:

1) input image and convert to grayscale then binary
2) using a disk structuring element perform a closing operation to remove the areas around the ornaments.
3) use a method to detect the circles (No problems here I think the main issue is in 2).

Is there a better way for me to isolate the red and gold ornaments? I just need a pointer in the right direction as I have tried a variety of things and none of them seem to work so I am starting to panic. I am working on this project in MATLAB and have included the current output below.

http://www.deviantart.com/art/Xmas-tree-png-3-HQ-large-398885414

The problem lies in step 1. Since you need the color information, you should not convert it to grayscale. Convert the image to HSV color space and work on the HUE(H) to filter the objects with color you require.


matched filter position I have seen that matched filter should be used in TED (gardner loop) after
resampler

resampler --------->   matched filter  ---------> down sample ---------> loop filter -> interpolation control for resampler
     (4 samples/symbol)            (4 samples/symbol)     (2 samples/symbol)

is this right?

Can matched filter be used with carrier recovery through costas loop?

The resampler would typically directly use an appropriately sampled matched filter; that makes the most sense, since that would actually give you the best crosscorrelation for timing correction.


Detection of a Pulse pattern with Matched filter I am working on a hobby transmitter-receiver project. Presence of 4 pulses in T seconds is regarded as 1(not any number of pulses, should match below):

lack of pulses during that T seconds is regarded as bit 0. This is how transmitter block sends out information.
This signal is transmitted over a lossy channel so what I receive is r = alpha.s + noise , which is then digitized and passed to a DSP. Noise is not necessarily Gaussian (not sure of it's distribution).
1- Can I use matched filter(done in DSP) to detect if s was the Bit 1 shown above or Bit 0?
2- When should I sample the middle of match filter response? Does this need to be synched to start of every bit?
I am not sure what other method could be robust and reliable for decoding this.

IF a 1 is always represented by the waveform in the figure, you can use that waveform as the matched filter response.
However, it sounds like the waveform could be anything as long as there are some pulses in it. In that case, the matched filter should be a single pulse, and you'll have to count how many pulses were detected afterwards.
The non-Gaussian noise can be a problem; if you don't know anything about its distribution, you'll probably just have to build a prototype and test it. You may end up having to design a filter tailored to the specific noise distribution in your system.