Power Spectral Density estimate (PSD) I'm trying to understand PSD, and, it's functions.
So far I have carried out the following:

Segmented a signal with NFFT = 256 and hop_size = 128
Windowed each of these segments (Hamming)
DFT of each Window
Taken the first half of the resulting matrix (n/2+1)

From what I have read so far, the algorithm can be best described as the following:
We take the squared magnitude of the LENGTH NFFT, we calculate this by: mag = (bin.re * bin.re + bin.im * bin.im). This gives the magnitudes of each of the points in the DFT bin.
Q1: These are "averaged" to form Pxx Pxx is a bin of size NFFT/2+1 which I assume store the averages, but, do these store the averages for each of the bins? Because I don't quite understand:
If I have a resulting 2D matrix containing values:
a = [[1, 2, 3, 4, 5, 6, 9, 8],
     [2, 5, 8, 9, 7, 8, 9, 7],
     [4, 5, 6, 4, 7, 8, 9, 8]

I then take the "average" of these.. I.e. (1+2....., n)/NFFT/2+1. This is only going to give one value, for each of the bins?
OK, so now onto Scaling. Can we therefore scale to the Frequency? If so, do we scale by Fs * window_size
EDIT:
I don't think my last question made that much sense..
So far I have a 2D vector containing 146x128..
block1 = [1, 2, 3, ...... n] (where `n` is 128)

block2 = [0, 5, 6, ...... n]

.....

block146 = [1, 4, 5, .... n]

Now, if I calculate the average magnitude of each of these blocks, i.e. sum(sqrt(re*re+im*im))/128 this will give me the averages for each block (a vector containing 146 samples).
For each of my samples in the blocks do I therefore divide each of these by the average magnitude for that particular block? I.e.
block1 = [0, 1, 2, ...... n]
average = (0+1+2)/128 = 0.58 (for example)

result = [0/0.58, 1/0.58, 2/0.58.... n]

Therefor my resulting vector will still be a 2D vector.
Does this make sense?
EDIT:

1 - What type of signal are you trying to process? I.e: What type of application is it? Is it a unidimensional signal, like Audio sampled at 44.1Khz? Or is it a multidimensional signal like an image or video?

I'm trying to process a Audio signal.
What I have done so far is:
Taken a 1D input vector, this vector is then split into blocks (of size 256 with an overlap/hop of 128 this gives me a total of 146 blocks with each containing 128 samples. These samples have been multiplied by the Hamming Window and then a DFT has been passed through.
These blocks contain the resulting values after Windowing and DFT has taken place.
2 - Where are you performing these calculations? Matlab, Octave, python, or some other?

I am performing these calculations using C++
3 - When you say 2D vector (maybe could be less confusing to call this a Matrix), What does the numbers inside mean? what are the columns and what are the rows?

The columns and rows contain the DFT results of the Hanning windowing of each block at 256/2 = 128
Essentially, what I'm attempting to do is create a spectrogram. Currently, I can plot the spectrogram with magnitudes of the values, i.e. (re*re+im*im) however, instead of plotting the magnitudes, I want to plot the overall power spectrum (PSD) so I can determine where the power is mainly within the signal.
If I do this using matplotlib I can get the following results:

So in theory, instead of showing the magnitudes which is what I'm currently doing, I want to show the PSD (like in the graph). I looked at the matplotlib source code and this is what was given:
result, windowVals = apply_window(result, window, axis=0,
                                  return_window=True)
result = np.fft.fft(result, n=pad_to, axis=0)[:numFreqs, :]

result = np.conjugate(result) * result

# for PSD
result /= (np.abs(windowVals)**2).sum()

result[1:-1] *= scaling_factor

# where scaling_factor is either 1 or 2

To me, the above code looks like they are taking the squared sum of the WindowValues (Which is Hamming) and dividing each value from the result (which is the resulting DFT) and finally multiplying by the scaling factor.
Does this look correct?

First of all, what do you mean when you say "Taken the first half of the resulting matrix (n/2+1)"?

The FFT will result in a vector, but maybe you mean taking only the first half of the FFT samples because the others are the negative frequencies and hence redundant?

Anyway, the process you describe is actually correct, and there is no need to do anything else. In your example with the resulting matrix a, you aparently took an FFT of size 8 (Or 16 if you then choped the second half of it). But what you do is simply take the average of all the first numbers of these vectors, and that would be an estimate of the Power spectral Density at f = 1*fs/NFFT.
Taking the average of all the second number will be an estimate of PSD at f = 2*fs/NFFT, etc.

In the example the PDS would be

PSD = 1/3*[(1+2+4) , (2+5+5) ,(3+8+6) , ... ];

When you are analizing in real time as you carry out the computation, if the signal is Wide sense stationary, then it should converge or hover around some value, but if it is not, it will deviate over time to adjust to the new coming statistics of the signal.

If i understood you correctly, that should be it.

Edit:  oops i forgot about the important scaling factor. The scaling factor should be introduced to preserve the parseval theorem between the PSD and its transform (The autocorrelation). If you windowed your signal first, i think you should apply the scaling factor 1/sum(W*fs/NFFT), where W is your window vector. Since you are merely dividing by the integral of the window function in the frequency domain This is what you add to the estimate of the frequency component  when you do the FFT by windowing first. Because the FFT alone already preserves the energy.

Edit 2:

I am even more confused with the 2-D thing. Let me try an ask you a few questions:

1 - What type of signal are you trying to process? I.e: What type of application is it? Is it a unidimensional signal, like Audio sampled at 44.1Khz? Or is it a multidimensional signal like an image or video?

2 - Where are you performing these calculations? Matlab, Octave, python, or some other?

3 - When you say 2D vector (maybe could be less confusing to call this a Matrix), What does the numbers inside mean? what are the columns and what are the rows?

4 - Did you check the link to mathworks? i think you will find what you are trying to do, is well explained there.

Also, just to make sure you understand the backgroun well enough, you could refer to some classical books that have the topic very well explained:

Digital Signal Processing, Principles, Algorithms And Applications -Proakis & Manolakis - 3rd Ed. Ch. 12

And the Father of all books in spectral Estimation:

Spectral Analysis of Signals - Stoica & Moses

Edit 3:

Ahh, ok, now i get what you are trying to do. So, here's the thing. I gave a quick read to the sections on FFT and PSD in the books by Stoica and Proakis, and still can't quite understand why they propose a scaling factor of simply 1/N to the FFT, while the matlab function for PSD using periodogram method scales the FFT by 2/(Fs*N), being N the length of the FFt, in your case 256.

So i don't know a lot of C++ but i would say do the following:

1 - Take first frame X, which consists of 256 samples of audio signal.

2 - apply 256 FFT to it. you obtain Xfft.

3 - To each complex value Xfft(i), for i = 0...255, do xfft(i)*conj(xfft(i))
which is equal to abs(xfft(i))^2. this you already did,so up to here nothing new.

4 - Multiply the whole vector Xfft by 2/(Fs*N).

This will give you the Power spectral density for your frame X. Now take the next frame and do the same, And when you reach a number of frames that you would like to plot, you stop and plot your vectors.

now since you are taking frames every 128 samples, your estimates at each time step will be sort of crammed up together, so a trick that youwould use is the fact that a Speech/audio signal is approximately wide sense statinary (WSS) in an interval of 20-50 msec (See Speech Enhancement - Jacob Benesty).
So what you could do, is assuming you have an Fs = 44100, that gives you a period of 50msec of around 2205 samples.

With the dimensions of the frames you took you could easily average 15 or 16 overlaped frames and get a Smoother estimate of the PSD. But the choice is yours, experiment with it and have fun!


Expectation of deterministic signals Maybe this a stupid question, but I came along with this point and I would like to confirm that I am not wrong.

I found some similar posts in StackExchange, like:

Expectation and Auto-correlation of an Output from a Matched Filter.

In one of the answer, it is stated "as S(t) is determistic[...]" and then the deterministic signal is taken out of the expectation.

and we know that for instance if I have a constant value, I can take it out of the Expectation operator, e.g.:

Where lambda is a random variable and "a" is a constant. Ok, so, if "a" is a constant, then fair enough, but what happens when "a" is a deterministic signal, but it is not a constant, for instance, what happens if "a" is a squared subcarrier?

Can I still safely say that the expectation will be the product of the deterministic signal multiplied by the expectation of the random variable?

Thanks, guys!

If $x(t)$ is a deterministic signal, then you have $E[x(t)]=x(t)$. Furthermore, if $Y(t)$ is a random process you have $E[x(t)Y(t)]=x(t)E[Y(t)]$.


Power Spectral Density (PSD) from autocorrelation I have a 1D fractional Brownian motion (fBm) signal, $u(x)$, of size $N$,  generated through a random number generator with a gaussian mean, $\mu=0$, and standard deviation, $\sigma=\sqrt{N}$. I want to find the power spectral density (PSD) of said signal. Signal data, MATLAB-R2018a code, and citation is included in the Github repository (https://github.com/dustinma324/Buruglence) .

From "Numerical Recipes" by Press 2007,
$$ Corr(g,h)_{j} = G_{k}H_{k}^* $$
where $G_{k}$ and $H_{k}$ are the discrete Fourier transforms (DFT) of $g_{j}$ and $h_{j}$, and the asterisk represents the complex conjugate. Normalizing the above equation by $N^2$ returns the PSD.

I am having trouble understanding the difference between the order of performing DFT. I have two road maps in mind to get the PSD. The first being,

$$ PSD = fft\Big(u(x)\Big) \cdot fft\Big(u(x)\Big)^{*} $$

where the PSD is found through breaking up the auto-correlation function into the product of the Fourier transform and the complex conjugate of the Fourier transform (convolution) of the two signals. The second,

$$ R_{11}(r) = u(x)*u(-x)$$
$$ PSD = \bigg|fft\Big(R_{11}(r=0) \Big)\bigg| $$

where the PSD is found by finding the auto-correlation function of the signal first, then taking the magnitude of the Fourier transform of the auto-correlation function.

The dark blue line is the second method, and red is the first. By taking the $FFT$ of the auto-correlation function, the general shape is captured, but the latter method shows small amounts of energy accumulation at the higher frequencies. This violates the Schwartz inequality, since the inner product of to vectors cannot be larger than the product of their magnitudes.

Is there a specific mathematical reason as to why the latter method appears to carry more energy at the higher frequencies?

Here's my MATLAB code below for the two different methods.

Method 1 (Convolution)

n = length(x); %Nx
pwr = [];
meanU = mean(u(:,tavg_start:end),2);
for i = tavg_start:1:nt
   U_fft = fft(u(:,i)-meanU,n)/n;
   pwr = [pwr U_fft.*conj(U_fft)];
end
meanpwR = mean(pwr,2);

Method 2 (Auto-correlation function)

pwR = [];
n = length(x);
meanU = mean(u(:,tavg_start:end),2);
for i = tavg_start:step:nt
    U11 = double(u(:,i)-meanU);
    Ucorr = xcorr(U11); % Output will be 2*length(Ucorr)
    N = 2^nextpow2(length(Ucorr));
    R_fft = fft(Ucorr,N)/(N^2);
    pwR = [pwR abs(R_fft)];
end
meanpwR = mean(pwR,2);

I have not been able to find anything that explains my particular question. Any inputs/suggestions/discussions would be greatly appreciated.

UPDATE Feb 4th, 2020

Auto-correlation function is analogous to a low-pass filter, meaning low frequencies are getting through, and high frequencies are attenuated out (gradual loss) at a certain cutoff frequency. The attenuated frequencies (unresolved energy) is being aliased into lower level frequencies when performing the Fourier transform of the ACF, resulting in the slight increase.

A transition band was described in this link (https://www.ni.com/en-us/innovations/white-papers/18/anti-aliasing-filters-and-their-usage-explained.html). I believe for me this band ranged from $k = N/2$, Nyquist frequency, to $k=2N/3$, where Basu dealiased the same signal using the 3/2 rule in the Fourier domain. This band is said to cause aliasing in the signal, thus showing up in my graph above.

Technically both methods shows some form of aliasing. So are there different kinds of aliasing? One is observed to aliased over a range of wavenumbers, while the other shows a spike.

Thanks,

Dustin

I am not sure I follow your exact questions.

  Is there a specific mathematical reason as to why the latter method
  appears to carry more energy at the higher frequencies?

No. The two methods are equivalent and will match to the numerical precision you are using.

  Auto-correlation function is analogous to a low-pass filter, meaning
  low frequencies are getting through, and high frequencies are
  attenuated out (gradual loss) at a certain cutoff frequency.

I do not understand this claim. How is auto-correlation a low-pass filter? Auto-Correlation is auto-correlation.

  Technically both methods shows some form of aliasing. So are there
  different kinds of aliasing? One is observed to aliased over a range
  of wavenumbers, while the other shows a spike.

Your input signal is already sampled and you are not changing the sampling rate, so I do not understand this either.

If the PSDs as estimated from two methods do not match, then there is a mistake in your code. Below is an example using octave/matlab code:

N = 1000;
x = rand(N,1);
Nconv = 2*N - 1;

% PSD from autocorrelation.
a = conv(x,flipud(conj(x)));
A = abs( fftshift(fft(a)) );

% PSD from FFT.
X = fftshift(fft(x,Nconv));
A2 = X .* conj(X);

% The residual between the two estimates.
err = 10*log10( norm(A - A2) / norm(A) );


Covariance matrix $P$ and $Q$ for complex valued data The equations for the Kalman filtering are

  IM = H*X;
  IS = (R + H*P*H');
  K = P*H'/IS;
  X = X + K * (y-IM);
  P = P - K*IS*K';

The covariance matrix in my application is of dimension 3 by 3 size because the there are 3 state variables that need to be estimated.

For data in the real domain, I initialized the covariance matrices $Q$ (process covaraince) and $P$ (prediction covariace) in Matlab as

P= 5*(eye(3,3))
Q  = 0.001*(eye(3,3));

For the complex domain, I initialized $P$ and $Q$ as

    P_real = 5*(eye(d,d));
    P_imag = 5*(eye(d,d));
    Q_real = eye(d,d)* 10^(-2);
    Q_imag = eye(d,d)* 10^(-2);

    P  = P_real + sqrt(-1)*P_imag;

    Q  = Q_real + sqrt(-1)*Q_imag;

R is a scalar real valued number.

Problem  : The default operation in Matlab is complex conjugate represented by the operator (.) '. I am using a toolbox which applies the above expression for Kalman filtering and there is no term conj() present especially when evaluating the expression IS , K, P. However, the toolbox never uses any (.).' or (.).* operators as well. When in complex domain, would the equations change to

IM = H*X  % no conjugate

IS = (R + conj(H)*P*ctranspose(H))

K = P*ctranspose(H)/IS;

X = X + K *(y - IM)  % here * denotes the multiplication operation

P = P - conj(K)*IS*ctraspose(K)

where H, X, y, P are all complex valued.

Are these okay?

In general, the equations from your first code block work equally for real and complex values.

Note that the ()' operation is the Hermitian Conjugate (i.e. transpose + conjugate). If you are in the real domain, it only becomes transpose (because conjugation on reals doesn't change them). Hence, you would not need to change anything in the code.

Regarding your covariance matrices: Note that a covariance matrix is always positive (semi-)definite and Hermitian symmetric, i.e. $P^H=P$. This ensures that each single variance on the diagonal is a real value. Also, how would you imagine a e.g. Gaussian with complex variance?

However, there is one thing which might be important for you here: When considering white Gaussian noise in the real domain, you generate random values according to

noise_real = sqrt(sigma2) * randn(N,M);

where sigma2 is the variance of this random variable. If you are in the complex domain, you have to generate them via

noise_complex = sqrt(sigma2/2) * (randn(N,M) + 1j*randn(N,M))

i.e. the variance becomes halfed for real and imaginary domain. Then, when you calculate var(noise_complex) you will get exactly sigma2.

I see from your post that your P and Q denote uncorrelated complex variables. In this case, if they are supposed to have a variance of sigma2 on each element, this is the corresponding code:

P = sigma2*eye(3);


How can I measure noise and sharpness of an image? I'm working on a project which ranked some images based on quality.

For this project, I want to figure out the noise and the sharpness from a image.

To calculate the  noise (from a CCD-Sensor) I think I can calculate the average noise from a image. How could I do this? I've read something about MSE, PSNR and SSIM.  However, all these require a second image. Or is there a modification, so as I can use these algorithms for one image?
To calculate a value which represent the sharpness I have no idea. Maybe I can calculate the peak with the sharpest point of the image (because there are different kinds like out-of-focus or motion-blur)

Is there an available reference paper or an implementation in any language?

Update
 Can I use (slanted-edge) MTF to measure sharpness? I need some normalized metric to compare sharpness between images. Or is there a better solution?  A stand alone Laplacian is not normalized, I think.

Simply, you can not 'calculate' the noise from single observation. Rather, you can estimate the target signal (denoised form), and then by using this estimation, you can determine the noise level. If the observation $y = x + n$ where the $x$ is target and $n$ is noise, after estimation of $x$ (via standart denoising methods), you can calculate the MSE, the SNR or other measures for noise level by using estimation $x$ and observation $y$. Of course these values will depend on the denoising method you use to estimate $x$, which can critically change the result.


Best way to approximate a curve? I have a curve like this one:
I need a function to approximate this curve, but I need that the function be a low order function (less than 5). What is a good way to obtain what I expect? Thanks!

So, you have a bunch of datapoints of the form (x,y), and considering all of those datapoints together, you have the vectors $x$ and $y$.

To do a curve fit, you would like to solve the equation $Aw=y$ for the column vector $w$, which holds the coefficients of your curve fit polynomial.  These coefficients are also the weights of the basis vectors in the matrix $A$.

The matrix $A$, for an $n^{th}$ order approximation, will hold $n+1$ columns and have a number of rows equal to the length of $x$ and/or $y$ (the number of datapoints).

The first column of $A$ is the vector $x$ raised element-wise to the power of 0. So, it's just a vector of 1's.  This is the basis associated with the offset or y-intercept.  The second column of $A$ is the vector $x$ raised element-wise to the power of 1. So, it's just the vector $x$.  The third column of $A$ is the vector $x$ raised element-wise to the power of 2.  You see the pattern?

To solve for the weight vector, $w=pinv(A)y$, the pseudoinverse of $A$ times $y$. Note, this doesn't actually solve the equation $Aw=y$, but performs an optimal projection of your data onto the vector space described by $A$.

EDIT:  $pinv(A)=(A^TA)^{-1}A^T$


Real-time BPSK demodulation I am currently studying BPSK demodulation. Most probably my question  very basic but i am having difficulties to understand the case of real-time bpsk demodulation. Basically, all the examples (matlab, theory books etc) assume i have the whole signal at the receiver with frequency and phase offset in which they apply frequency and phase offset estimation and later symbol and timing recovery. I find it difficult to realize this scheme in a real-time fashion, in which for example I constantly receive samples from antenna in a buffer.

How can i start to perform bpsk demodulation without isolating the signal since the the timing recovery is performed latter in the examples i have checked?

For sure I am missing something or confusing the concepts.
Can you please provide a toy example of such methodology or point me to a usually technique to perform in a real-time implementation.

Kind Regards,
J

In the implementation of a real receiver, you have to do what is called "streaming", or handing of samples from one block to the next in real time.

Matlab does not lend itself to working with streaming samples; you need a different tool.

For instance, GnuRadio is designed to work in real time with streaming samples. You can even design your demodulator using a visual tool with drag-and-drop of different modules.

It should also be possible to use Python, by setting up a pipe to read samples from your front-end and processing them in sequence.

If you want to use a DSP, I can recommend this book: http://amzn.com/0387748857 which describes the entire design of a streaming digital communications system using a DSP evaluation board.


Shelving filter (2nd order) I'm trying to find the equation of a shelving filter of the second order. I easily found the equation for a first order filter on wikipedia with a transition region which have a 6db per octave slope. The equations are :

coef^2 * (1 + (f / (coef * f0))^2) / (1 + (f/f0)^2) for the lowshelf
(1 + (f / (coef * f0))^2) / (1 + (f/f0)^2) for the highshelf

With :

coef = multiplicator coef (gain = 10 log10 (coef))

f = frequency

f0 = cutoff frequency

But now, how can I use a kind of Q factor to modify the slope of the transition region ?

Of course, as I'm using already a frequency signal and not a temporal signal, I can't use the equations made for the z transform...

Thank you for all your answers

Ok, I finally found a solution after few search about sigmoïd charts and few try.
What I wanted is to modify the slope factor in the transition area without modifying the high and low border values.
So, the equations are becoming :

coef^2 * (1 + (f / f0)^(SFactor+2) * (1/coef)^2) / (1 + (f/f0)^(SFactor+2)) for the lowshelf

(1 + (f / f0)^(SFactor+2)) * coef^2 / (1 + (f/f0)^(SFactor+2)) for the highshelf

With :

coef = multiplicator coef (gain = 10 log10 (coef))

f = frequency

f0 = cutoff frequency

SFactor = slope factor. If SFactor = 0, it acts like a classical first order shelving filter, if SFactor increases, the slope in the transition area is steeper.

Thank you for all the previous comments and the Audio EQ CookBook is also a very interesting source I'm going to study...


Time domain representation of a modulated signal when symbol rate >> carrier frequency? I need to understand the relationship between carrier vs. bits with PSK modulation, and what the signal in the time domain would look like if my symbol rate is far greater than my carrier frequency. How can I put more symbols on a carrier which is at a set frequency? For example how would I put 1000000000000000 bits on a 1 Hz carrier wave, is this possible and what would that look like?

tl;dr: you can't. That's basics!
longer story:

You ask about symbols.
When the symbol rate is much larger than the carrier frequency, you've basically lost, and cannot transport that many symbols per second; you don't have the bandwidth; you get at most twice the carrier frequency in bandwidth, namely from 0 Hz to twice the carrier frequency.

However, you're asking how to modulate 10¹⁵ bits on a 1 Hz carrier, and that's a different question: theoretically, you can take forever to transport these 10¹⁵ bits, so there's no problem here. Realistically, we'll all be dead by the time you've finished your transfer, but whether or not that is  a problem is of mere philosophical interest.

If you, however, mean you want to transport 10¹⁵ bits per second on a 1 Hz carrier wave, you have less than two channel access per second, so you'll need enough SNR that you get more than  $2^{5·10^{14}}$ different constellation points.
If you reasonably try this, I hope someone stops you in time: as physics dictates the difference between two detectable constellation points needs to be at least the Planck constant, in fairly good approximation, the action needed to get from the all-0 bitword to the all-1 bitword would be at the very least,
$$\frac{2^{\left(5\cdot 10^{14}\right)}\cdot h}{\pi} =\frac{2^{\left(5\cdot 10^{14}\right)}\cdot 2^{-110.2}\frac{\text{kg·m²}}{\text{s}}}{\pi} \approx 2^{\left(4\cdot 10^{12}\right)} \frac{\text{kg·m²}}{\text{s}}\approx 10^{4\cdot10^{11}} \text{J·s}.$$
That is incredibly much action. This is nuts. Your photon of 1 Hz would have to have an energy of $10^{4\cdot10^{11}}\,\text{J}$ to be able to walk a path of that action; that's as much energy in a single photon as the sun emits in 487 years. This photon would end the universe as we know it. There'd be only heat and chaos. And then, who the hell will care about your 10¹⁵ bits?


Fourier Transform within a certain Limit I want to evaluate fourier transform within a certain limit in MATLAB,the expression of which is
$$X(f) = \int_{1}^{4}{x(t)e^{-i2\pi ft}}\,dt$$
I have to find value of the above expression within limits which are definite in nature.
I came across this post on MATLAB discussion forum which says to multiply the function by dirac(x-lowerbound) * dirac(upperbound-x) and then use fourier to calculate the desired results,but i am not getting the desired results.
What else can I try?
https://in.mathworks.com/matlabcentral/answers/16049-need-command-for-continuous-time-fourier-transform

So, your integral with bounds

$$X(f) = \int_{a}^b x(t) e^{-i2\pi ft}\,dt$$

could really be written as

$$X(f) = \int_{-\infty}^\infty (r_{a,b}(t)\cdot x(t)) e^{-i2\pi ft}\,dt$$

with $r_{a,b}(t) = \begin{cases}1 & a<t\le b\\0&\text{else,}\end{cases}$  i.e. as a "windowed" view at $x(t)$.

Luckily, the convolution theorem of the Fourier transform allows us to say that

$$\mathcal F\{r_{a,b}\cdot x\} = \mathcal F\{r_{a,b}\}*F\{x\}$$
with $*$ being the convolution operator.

Thus, you'd just need to numerically determine an estimate for the Fourier transform of $x$, $F\{x\}$, and convolve that with the Fourier transform of your window $r$. Luckily, the Fourier transform of a unity-wide rectangle is known (it's a sinc function), and you just need to apply the scale and shift theorems of the Fourier transform to convert that standard sinc to the Fourier transform $F\{r_{a,b}\}$.

Then you'd need to estimate the convolution of the two.

Assuming your signal is well-behaved (especially: band-limited), a discretized signal $x[n]$ represents your $x(t)$, and hence, you could estimate the continuous time Fourier Transform within the signal bandwidth by doing a DFT (matlab: fft()) of sufficient samples of a period of $x(t)$.

As said, finding $F\{r_{a,b}\}$ is easy to do analytically. Simply generate the same amount of samples of $F\{r_{a,b}\}$ as you did for $x(t)$, and numerically convolve the result (matlab: convolve, filter).

This might be much more effort or much less elegant or much less accurate than simply doing numeric integration, which matlab can also do out of the box. But I'm assuming you're not interested in that, because a) it's trivial to find "matlab numerical integration" if you just look for it, so I presume you already did that and b) I'm assuming you're asking how to solve this with signal processing methods, considering the place you ask.


Why dft of moving average impulse response has 0 angle I wanted to visualise angle and magnitude of moving average DFT. I used python3 for that:

import numpy
irs = [1 / 5] * 5
dft_irs = numpy.fft.fft(irs)
print("dft", dft_irs)
print("magnitude", numpy.absolute(dft_irs))
print("angle", numpy.angle(dft_irs))

I got results that I do not understand

dft [  1.00000000e+00+0.j   0.00000000e+00+0.j   1.38777878e-17+0.j 1.38777878e-17+0.j   0.00000000e+00+0.j]
magnitude [  1.00000000e+00   0.00000000e+00   1.38777878e-17   1.38777878e-17 0.00000000e+00]
angle [ 0.  0.  0.  0.  0.]

Here phase shift is 0 for every frequency while I know it is not for moving average filter. Please help me understand why I got this strange result.

You are looking at the frequency response at multiples of $2 \pi /5$, where the response is either 1 or 0.

If you add several 0 values to the end of irs, you'll get the frequency response at more intermediate frequencies, where you can see that phase is not zero.


What is the difference between the Weiner-Hopf equation and the normal equation? In adaptive filters,  the development of LMS algorithm typically starts from the Weiner-Hopf equation, while the development of RLS algorithm starts from the normal equation. As I understand, these two equations are the same, and both their solutions is the optimal coefficients that the adaptive filter has to find. Is this true or am I missing something?

Weiner-Hopf equation leads to Wiener filter that is optimal filter. For the case of stationarity in some time span it's the only filter minimizing MSE at its output. LMS and RLS algorithms are the adaptive approaches and they converge to Wiener optimal solution (as you can see from their convegence curves). While LMS is the simple steepest descent method of finding such an optimum, RLS solves Least Squares estimation problem at every step recursively.

LS estimation problem can be treated as a practical approach of Wiener filter. Look:

$w = \hat{R}^{-1} \cdot \hat{r}$

is the normal equation in matrix form for the LS problem, $\hat{R}$ and $\hat{r}$ are estimations of autocorrelation matrix and cross-correlation vector respectively. So if these estimations is properly done this solution leads us to the Wiener-Hopf eq., i.e. to Wiener filter. But it's only true if estimations are valid. So in general LS solution is not a Wiener solution, it's an approximation.

RLS approach is derived from LS (as it can be viewed from abbreviation):

$w(k) = S(k) \cdot \hat{r}(k)$,

where $S(k)$ is recursive estimation of $R^{-1}$ and $\hat{r}(k)$ is a recursive estimation of $r$. From the equation above one can see RLS try to solve LS problem at every step.

So the conclusion might be both LMS and RLS try to converge to Wiener solution and Wiener filter theory is in the core of these algorithms. But ways they do this are different.

This is my own understanding, maybe some discrepancies are possible. There is very good reference about adaptive filters: http://books.google.ru/books/about/Adaptive_Filter_Theory.html?id=MdDi_PF7gMsC&redir_esc=y


How to compute fundamental frequency from a list of overtones? Given a list of overtones (F1, F2, F3, etc), how do I compute the fundamental frequency? Can I do something like F2/F1=F1/F0? Is it the correct method to use?

Nope. Difference between overtones is a good point to start, i,e F3-F2, F2-F1. The differences should be all the same or multiple of each other. The smallest one is often the fundamental. It gets more tricky of the spectrum is "sparse", i.e. a lot of the harmonics are missing. Then you need to find a largest possible divisor that turns all frequencies into integers or, to be precise, so that the ratio of frequency to fundamental is within the measurement accuracy of the nearest integer.


DFT and multiplication/convolution equivalence Is there a simple or potentially intuitive explanation for, with the DFT, vector multiplication in one domain being equivalent to circular convolution of the transforms of the vectors in the other domain?

Since a DFT is just multiplication by a (special) square matrix, what about this matrix and matrix multiply allows the above duality?

As you correctly say, the DFT can be represented by a matrix multiplication, namely the Fourier matrix $\mathbf{F}$. On the other hand the DFT "transforms" a cyclic convolution in a multiplication (as all Fourier transform variant as DFT, DTFT, FT have a similar property of transforming convolution to multiplication) and vice versa.

To understand this in the matrix picture, note that also (circular) convolution with a certain sequence can be represented by a matrix multiplication. More specifically this is a circulant matrix, a special kind of a Toeplitz matrix.

so $\mathbf{y} = \mathbf{c} * \mathbf{x}$ with $*$ the cyclic convolution can be written as
$\mathbf{y} = \mathbf{C}(\mathbf{c}) \mathbf{x}$ with $\mathbf{C}$ denoting the circulant matrix formed from entries of vector $\mathbf{c}$.

If we "transform" this equation with the DFT (i.e. multiplication by $\mathbf{F}$) we obtain

$\widehat{\mathbf{y}} = \mathbf{F} \, \mathbf{C}(\mathbf{c}) \,\mathbf{F}^H\, \widehat{\mathbf{x}}$

with $\widehat{\mathbf{y}}  =\mathbf{F}\mathbf{y}$ and $\widehat{\mathbf{x}}  =\mathbf{F}\mathbf{x}$ the respective DFTs (note $\mathbf{F}^H$ represents the IDFT).

The point is now that $\mathbf{F} \mathbf{C}(\mathbf{c}) \,\mathbf{F}^H$ is always a diagonal matrix, because all circulant matrices are diagonalized by the Fourier matrix. This means that the eigenvectors of circulant matrices are just given by the rows of the Fourier matrix.

This is of course consistent with the convolution picture,
because the DFT transforms the convolution to an elementwide multiplication. Moreover the diagonal elements of this matrix are just the DFT of $\mathbf{c}$, or, eqivalently, the eigenvalues of the circulant matrix formed from $\mathbf{c}$.


What is impulse response if i know input and output? I would like to find the impulse response, h[n], of an LTI system given the input

x[n]=[1,−3,2]

and the output

y[n]=[1,−1,−4,4]

I know that y[t]=x[t]∗h[t], but I am having hard time to figure out the right way to calculate the impulse response. I know very little about signal processing, so if you don't mind giving an easy explanation, then I appreciate it. Or, if it's possible to do an example, that's better.

If $x[n]$ is non-zero for $n=0\ldots2$, and $y[n]$ is non-zero for $n=0\ldots3$, then $h[n]$ must have non-zero components for $n=0$ and $n=1$. This is a property of the convolution.

So now you have to solve $y[n] = h[0]x[n] + h[1]x[n-1]$ for the relevant values of $n$, and there you go.

Start with the edges ($n=0$ and $n=3$).


What are some typical lossless compression ratios? A client was trying to send me 250 GB worth of files. After attempting various ways of sharing the data, he sent me a zipped folder only 4 GB in size. That sounds like too much compression to me--I don't think when I've zipped things I've ever reduced more than 20% of the size.

What are some typical lossless compression ratios that one would see in practice? (Or alternatively, a range.)

UPDATE: I realise it's impossible to say without being able to guess at the actual information content, so perhaps this was an ill-worded question. I can't share the client's data. But looking at the XML files, there are a lot of repeated phrases e.g.

<thing>
    <property="1" value="2" />
    <property="3" value="4" />
    <property="5" value="6" />
    <property="7" value="8" />
    <property="9" value="10" />
    <property="11" value="12" />
    <property="13" value="14" />
</thing>

which does seem quite compressible.

It really depends on how much redundancy there is in the data. If all 250GB was just '0', then you could get fabulous levels of compression.

This page shows some results for compression of English text. It compresses 2,988,578 bytes worth of text using various techniques. The top three are: 330,571 (88.94%), 333,759 (88.83%), and 352,722 (88.20%).

Relating this back to your case: that means 250GB would compress (at best) to 27.65GB.

So perhaps there are lots of blank, repetitive patches in the data?

UPDATE:

Based on your edit, here is another paper that looks at XML compression rates.
The upshot (Figure 8; screenshot below) is that they appear to see compression rates ranging from 1 (none) to about 50 at a maximum.  That suggests that you are correct to be suspicious of your customer's 62.5 to 1 (250 to 4) compression rate.


Difference Between Dirac delta and unit impulse function The definition of Dirac delta function states that it gives a value of $\infty$ at t=0 and 0 elsewhere.

But, the definition of unit impulse function states that it gives a value of 1 at t=0 and 0 elsewhere.

Many textbooks state that Dirac delta and unit impulse are the same function.

So What is the difference between Dirac delta and unit impulse function?

The Dirac delta $\delta(x)$ is a continuous argument (generalized) function where $x$ is a continuous variable and can be time, space, frequency etc. Unfortunately, a simple definition like

$$\delta(x) =
\begin{cases}
\infty ~~~&,~~~ x=0 \\
0 ~~~&,~~~ x \neq 0
\end{cases}
$$

does or can not properly describe its behaviour. You need to define its properties under the integral sign, and that requires advanced mathematical stuff for a formal verification. But  a simplified Riemann definition of $\delta(x)$ is given as:

$$\int_{-\infty}^{\infty}f(x)\delta(x-a) dx = f(a)$$ for a sufficiently smooth test function $f(x)$. All properties of $\delta(x)$ follows from this Riemann (seeming) integral operation.

On the other hand, the discrete counterpart, the unit-impulse or the Kronecker delta is a sequence of integer index $m$ with very well defined behaviour as:

$$\delta[m] = \begin{cases}
1 ~~~&,~~~ m=0 \\
0 ~~~&,~~~ m \neq 0
\end{cases}
$$
These two functions serve the same purpose for discrete and continuous systems, hence share the same name impulse. Yet they are two different, distinct entities; the Dirac impulse is a very problematic concept, whereas the discrete  impulse is a mere simplicity.


Hartley Transform vs Fourier Transform Can you explain to me in what way Hartley Transform differs from Fourier Transform? Is it even used today or is it some mostly forgotten, obsolete archaic thing?

Please don't use equations - I dont understand them, but I understand Fourier Transform.

You are probably interested in the real discrete transforms. Here are the basis functions for a transform of length 6. The inverse transforms are weighted sums of these basis functions, and the forward transform figures out the weights given the time domain sequence. A reconstruction of the time domain sequence by the inverse transform is exact. The discrete transforms don't actually involve the continuous functions shown below but rather their samples at integer time indices (here 0, 1, 2, 3, 4, 5).

As a service to those who prefer to see equations, the discrete Hartley transform $H_k$ and discrete Fourier transform $X_k$ of $x_k$, where $k = 0, \dots, N-1$ are defined by:

$$H_k = \sum_{n=0}^{N-1} x_n \left[ \cos \left( \frac{2 \pi}{N} n k \right) + \sin \left( \frac{2 \pi}{N} n k \right) \right]\\
X_k = \sum_{n=0}^{N-1} x_n \left[ \cos \left( \frac{2 \pi}{N} n k \right) - i \sin \left( \frac{2 \pi}{N} n k \right) \right],$$

where $i$ is the imaginary unit.

Hartley

Figure 1. Hartley transform basis functions are harmonic sinusoids. Each basis function is a sum of sine and cosine.

Being harmonic means that the frequency of each is an integer multiple of the lowest non-zero frequency. Some of the sinusoids have frequency greater than the Nyquist frequency.

Or as discrete sequences: [1, 1, 1, 1, 1, 1] (yellow), [1, 1.366, 0.366, -1, -1.366, -0.366] (black), [1, 0.366, -1.366, 1, 0.366, -1.366] (green), and so on.

Figure 2. Hartley basis functions, constructed from the discrete basis functions by perfect bandlimited interpolation assuming that the sequences repeat forever. The discrete Hartley transform could alternatively be expressed directly with these functions, colored the same as in Fig 1.

Except for the special cases of frequency 0 and π, the aliased Hartley basis functions come in pairs of harmonic equal frequency sinusoids that are 90 degree apart, just like in the real discrete Fourier transform below. Except for scaling and phase shifts (with one scaling constant and one phase shift constant for each pair), the basis functions are identical to those of the real discrete Fourier transform.

Fourier

Figure 3. Real discrete Fourier transform basis functions are harmonic cosines and sines, here colored the same as the equivalent scaled and phase-shifted Hartley basis functions in Fig. 2.

For frequencies 0 and π there are no sine basis functions because their samples would be zero-valued and would never contribute to the weighted sum in the inverse tramsform. Unlike with the Fig. 1 formulation of Hartley transform, the discrete Fourier transform basis functions stay within the Nyquist frequency limit.

Popularity of Hartley

The fast Hartley transform algorithm came out in 1984, but there has been no significant increase in its popularity since.

Figure 4. Google Ngram of Hartley transform counts the relative frequency of the word pair in literature as it developed over the years.


Finding linear convolution of two time series PROBLEM

Two real, causal time series

$$f_k = a^k \quad \text{and} \quad g_k = b^k$$

where $a \neq b$, $|a| < 1$, and $|b| < 1$, are given for $k = 0,1,2,...$

Find the linear convolution $h_n$ of the time series $f_k$ and $g_k$.

ATTEMPTED SOLUTION

I think perhaps my attempted solution here can be improved upon. We have from the definition:

$$h_n = \sum_{k=0}^{n} f_k g_{n-k} = \sum_{k=0}^{n} a^k b^{n-k} = b^n \sum_{k=0}^{n} \left(\frac{a}{b} \right)^k$$

If we let $R = \left(\frac{a}{b} \right)$, the last term can be written as:

$$h_n = b^n \cdot \frac{1 - R^{n+1}}{1 - R}$$

or

$$h_n = \frac{b^n - \frac{a^{n+1}}{b}}{1 - \frac{a}{b}}$$

$$h_n = \frac{b^{n+1} - a^{n+1}}{b-a}$$

I don't see how I can improve upon this answer though, but since we are given certain conditions in the problem, I have a feeling that it is possible to get a more elegant solution.  I see that as $n \to \infty$, $h_n \to 0$, but is there any way I can "improve" on the answer for any $n$?

There is a small mistake in your solution. Here $f_k$ and $g_k$ are infinite duration signal, ie. $n \to \infty$. You made a mistake that you took  n in $h_n$ and number of coefficients in $f_k$ or $g_k$ as same, but it is different.

$$h_n = \sum_{k=0}^{\infty} f_k g_{n-k} = \sum_{k=0}^{\infty} a^k b^{n-k} = b^n \sum_{k=0}^{\infty} \left(\frac{a}{b} \right)^k$$

$$h_n = b^n \cdot \frac{1}{1 - \frac{a}{b}}$$

$$h_n = \frac{b^{n+1}}{b-a}$$

If you want to finite the number of coefficients in $f_k$ and $g_k$, use other variables other than n.

Then n in the RHS of your equation will replace with the new variable.


Is the Channel shortening similar to Time domain Equalizer (TEQ) I used to read about Time domain equalizer (TEQ) in many lectures, i.e here , But in many other papers,  that's called  Channel shortening i.e here .

However, in the a two above links, the cost functions is almost similar, but in first paper it's called Time domain equalizer, and in the second, it's channel shortening.

What's the difference between them (I mean between Time domain equalizer TEQ and Channel Shortening)?  And why do we call it sometimes TEQ and another time channel shortening?

Thank you.

Most likely people use different words to describe the same concept. This happens quite often. The TEQ is the more general term and also used to equalize other methods that OFDM, while CS is related to OFDM and tries to make the channel shorter so that it fits within the CP duration.

You can check online, you will find that almost all articles related to CS are in OFDM whereas the TEQ is a general term

Good luck


dBFS scaling and spectrum I am trying to replicate the spetrum as given in Audacity and other commercial software... I am getting crazy with the correct scaling needed in order to have dBFS (with reference to a full scale sine wave).

Suppose $x_t$ is the PCM sample at $k$ bit precision, first of all I rescale the PCM samples

$$
s_t =  \frac{x_t}{  \frac{ 2^{k-1} }{ \sqrt{2}  }  }
$$

Suppose now $w_t$ is some window function, I compute the FFT of of size N of the product $s_t w_t$.

Let $P_1, P_2, \ldots, P_N$ be the FFT complex values. Since I am only interested in one-side spectrum, I just retain the fist half of the sequence. Moreover suppose for simplicity that the time series matches the window size $N$.

Taking

$$
p_i = 2\left| \frac{P[i]}{N}\right|^2
$$

for $i=1,2,\ldots, N/2-1$ should give the power spectrum values. Then one would take $10\log_{10}(p_i)$ to get dBFS.

But I think that in order to express dBFS one needs to scale the $p_i$ in some way taking into account the effect of the window. In fact  I produced a pure 1KHz sine wave then I computed my spectrum and compared it with that obtained with Audacity,  Adobe Audition etc. My FFT stays below by few dBFS so I think it's really a problem of scaling.

I hope somebody knows if there is a standard ways of doing that!

Regards
Pierre

Here's a paper with a pretty thorough explanation of proper scaling in order to obtain the power spectrum or power spectral density. Essentially N should be the ENBW or the sum of the window squared.


What type of images are required for voxel classification? In an algorithm, voxel classification approach is used to perform image segmentation. The aim was classifying voxels either as cartilage or background in MRI image. The project has to be completed in matlab.

MRI images has .dcm extension, it may have .avi and other types too, I am not sure.  Is it possible to perform voxel classification on .dcm images or it require another type of images?

In general it doesn't matter what the format is, because you will load it as array/matrix anyway.

That being said, I never saw preprocessing being done on DICOM images, usually they are transformed to nifty or something else before anything is being done to it. However, ignorance is not an argument.


Can the inverse system of a stable system be unstable? Can the inverse system of a stable system be unstable?

For the class of LTI systems, the criteria for stability of a system with impulse response $h(t)$ and systems function $H(s)$ are:

$h(t)$ be absolutely integrable.

or in s-domain:

$H(s)$ have the $j\omega$ axis (or $s=0$) in its convergence region, and do not have any derivatives of impulse in it (which means for rational system functions like $P(s)/Q(s)$ the order of the nominator be less than or equal to the denominator).

So I thought that if we consider the following system
$$H(s)=\frac{s-2}{s^2-1}$$
with $ROC:s>1$, it is stable, but its inverse $H_i(s)=1/H(s)=\frac{s^2-1}{s-2}$ with $ROC:s>2$ has the first derivative of delta in it and thus is not stable.

Is my reasoning correct? and is there any way to answer the question for all of the systems (not only LTI)?

The ROC is determined by a condition on the real part of $s$, so for the given transfer function $H(s)$ you could define the ROC as $\Re\{s\}>1$, but then the system wouldn't be stable because the ROC does not include the $j\omega$ axis. You could as well have the ROC $-1<\Re\{s\}<1$. In this case the system would be stable but not causal.

And now for your main question: yes, the inverse of a stable system can be unstable. If we consider systems described by rational transfer functions (i.e. lumped elements) then the location of the system's poles determine its stability. The poles of the inverse system are the zeros of the current system, so there's no relation between the stability of a system and the stability of its inverse, because the zeros do not at all influence the system's stability, but they do determine the stability of the inverse system.

For your example, the inverse system has a pole at $s=2$, and there are two possible ROCs: $\Re\{s\}>2$, which corresponds to a causal but unstable system, and $\Re\{s\}<2$, corresponding to a stable but anti-causal system.


Why the energy is a good feature extraction for detect disturbances in signal processing? To detect disturbances in signal processing a common step is to extract signal characteristics to analyze them, among these characteristics it is recommended to use the signal energy.
Why energy is recommended over other features that can be extracted, such as statistical parameters of the signal or frequency components?

In short, energy is not a good feature per se. Disturbance enhancement sometimes requires preprocessing, and might not be resistant to outliers.
The notion of disturbance in your  question could be made more  precise.

However, since the energy squares amplitudes of signal values, and is generally integrated over several samples, it has the ability to emphasize disturbances whose amplitude dynamics vary from the "normal" ones.

Aside, within the scope of automated detection, energies are quantities that are relatively easy to minimize: least-squares are common techniques in signal processing.

Yet I would not consider energy as a good feature in general. Depending on the need, other nonlinear operators can be used, as can be seen in the family of Teager Kaiser operators, see for instance: Signal processing using the Teager Energy Operator and other nonlinear operators.


Shape recognition, hough transform I'm trying to match a picture saying "this is ducks, this is towers" by giving votes  from the hough transform of a picture.

I've found this article http://ttic.uchicago.edu/~smaji/papers/hough-cvpr09.pdf, but I didn't find a implemented code (and I'm trying several things).

Do you know anything similar to that, for example included on OpenCV?

Thanks.

Hough forests for object recognition are initially proposed by Jurgen Gal. The paper was named "Class Specific Hough Forests for Object Recognition".

Here is his homepage. Here are the source code and the details. I think you could make use of those.

Certainly better methods exist and one could well talk about them. However I am sticking to the Hough framework, as you mentioned that in the first place.


Envelope of I/Q data In finding the envelope of the I/Q data set (which is a complex data set, Inphase /Quadrature data) I cannot use envelop(x) in MATLAB, as it requires the real-valued "X" in MATLAB, so got an error. Does anyone have a suggestion for solving this, finding envelope from I/Q data? Thanks,

If your signal is an analytical signal, then envelope(real(X)) would simply be your answer. This would be identical to using the Hilbert transform to extract the Quadrature part of the signal from the Inphase part. If your signal is not analytical, your envelope would simply be:
$$G = \sqrt{I^2 + Q^2} $$


What is the ideal type of barcode for paper tickets? I can use the following type of barcodes in my project:

Code 39
Code 128
QR Code
Data Matrix
PDF-417

They will be printed with standard office printers and with boca printers. It should also be possible to show them on smartphone screens.

They will be read with webcams or Windows tablets. It would also be helpful I could buy an USB barcode reader or something similar to use with this codes.

The data in the code will be text and numbers, at least 10 characters.

Which barcode type is most suitable for this situation? The most critical factor is the fault toleration, this has to be as low as possible on the mentioned media.

First, about the USB barcode reader: at my previous job (QA at a label printer company), we had one from Scanology. It worked as a keyboard device, the scanned data is sent to the pc as if it were typed. So you don't need any special software.

Code 39 has no error correction, but even a single misread line won't result in wrong data, because it simply won't produce a valid character. So you either get valid data, or you get nothing.
A major advantage of code 39 is that you don't need barcode software, you can simply use a barcode font.
Code 128 is slightly better than code 39 in terms of error detection, but still no fault tolerance. And you can't use a barcode font. In your case, you would use 128A (0-9, A-Z and a few special chars).
If your most critical factor is fault tolerance, then I would suggest QR, with high error correction setting, which gives up to 30% redundancy. This means that up to 30% of your QR code can be unreadable before the scanner gives up. The Wikipedia page on QR shows an example of a torn QR code that is still readable: http://en.wikipedia.org/wiki/QR_code#Error_correction
Another advantage is that it can be scanned in any direction, and it's a quite popular format.
Data Matrix and PDF-417 are similar to QR in terms of fault tolerance, but they are not widely supported for smartphones. I scanned thousands of QR and PDF-417 codes during my time as QA engineer at that label printer company, and I found a lot of high quality QR scanning software for Android and iPhone, but hardly any for PDF-417 and Data Matrix. I assume the situation will be the same (or worse) for Windows Phone.

My advice:

Easy implementation: Code 39. Use a Code 39 barcode font. Works with cheaper barcode scanners.
Very high fault tolerance: QR at level H (High). Use the ZXing library. Handheld USB barcode scanners may cost a tiny bit more.


Combining two separate radar signals? I have two small radar sensors that use 25 MHz. From each sensor I get the sent signal and received signal (Doppler effect). I know that I can use the FFT, use a window, or FIR filter, take the inverse, and get the frequency within that window. Thus I can analyze the signals separately. I'll have two plots. The radars can be positioned in 3D space at different points, but they can be pointed at the same point in space or at different points.

Is there any way to combine the two signals to analyze motion? For example if one sensor is pointed at the cheek and the second sensor is pointed at the chin, can I somehow average the signals to get one signal that represents the observation of the cheek and the chin? Can I average the two separate plots together?

Or say we had a spinning top. One sensor is on the right and another on the left. If the top leans right the sensor on the right gets less of a time delay and the sensor on the left gets a longer time delay for the received signal. Thus using the information from both I can tell at what angle the top is spinning hypothetically because speed is distance over time.

Any help, math equations or references would be appreciated.

Yes, there is a way to combine radar sensor data.  There are instances where you can correlate the two signals with each other, but it depends on what you are trying to do.  A couple of points about your examples:

Chin/cheek.  At 25 MHz, any reasonable aperture (antenna) will likely have a large enough beamwidth that both sensors would see the same thing: both cheek and chin.  That said, you may see a Doppler shift from both objects.

Top.  Conjugate one of them, provided that all other parameters are exactly the same.  That said, one sensor is already measuring the exact event that you describe.  If you need higher SNR, there are other things you can do.

And don't forget that each radar can pollute each other's data


I cannot find reference (paper) of this relation $u(t)t^{α} ↔^{FT} f^{-(α+1)}$ I am struggling a lot in finding the reference (book or, even better, academic paper) that states that:
$u(t)$     $t^{α} ↔^{FT} f^{-(α+1)}$
as I found here An Interesting Fourier Transform - 1/f Noise in the formula in the first box of the link.
I searched a lot, but I did not find any reference, and they are very important to me.

You will find a "proof" of principle in math-stackexchange: Fourier transform of power function $t^\alpha$ (the title lacks the mention of the unit step function $u(t)$). At least formally, the Fourier transform of one-sided power law turns into a Laplace transform, and looks a lot like the Gamma function (the one related to the factorial).
You can find a similar description on page 276 of the chapter  Appendix A. Mathematical Background, Power law.
What I find troubling in this "simple proof" is that the convergence and well-behaved properties of the power-law functions are not satisfied. They are not integrable in the traditional sense, so it should be treated with tempered distributions or generalized functions, and Schwartz functions decaying faster than powers.
You will find a clean demonstration in Terence Tao blog: 245C, Notes 3: Distributions. 3. Tempered distributions, around Equation 9. In dimension $d$ ($d=1$ for you), the formula reads:
$$\widehat{u(t)t^\alpha}(f) = \frac{\pi^{-(d-\alpha)/2}\Gamma((d-\alpha)/2)}{\pi^{-\alpha/2}\Gamma(\alpha/2)}|f|^{-(d-\alpha)}$$
[Nota: this is a follow-up of question: Can I simplify $x=\frac{\ln(|fft((ifft(\sqrt{(f^{-5/3})}))^{2})|)}{\ln(f)}$]


If the impulse response of a causal LTIC system is $h(t) = \delta(t) + \sin(t)u(t)$, is it marginally stable or unstable (BIBO)? If you take the $H(s) = \mathcal{L}[h(t)]$ the poles are on the imaginary axis so the system should be marginally stable, but is it?

The Laplace transform of the given impulse response, i.e., the transfer function of the corresponding system, has two poles on the imaginary axis. According to the definition of marginal stability, that system is marginally stable.
It is, however, not BIBO stable, because given an input with a frequency equal to the frequency of the poles will result in an unbounded response. For a causal continuous-time LTI system to be BIBO stable, all its poles must be strictly in the left half-plane, i.e., their real parts must be negative.


how gabor filter is used on image texture edge detection as I searched gabor filter is most used tool on texture analysis and edge detection in digital image. its said that it searches a certain frequency as a texure around a pixel or point to find the texture. I cant really understand it well. as I know frequency in digital image is sign of change in grayscale change. high frequency has high change in grayscale so a highpass filter finds the edges but I cant understand a filter like gabor can find a special texture.

thank you for helping

A periodic geometrical texture may appear as a localized wiggling pattern, with dominant orientations and 2D frequencies. Its correlation with Gabor atoms is bigger when the atom has orientation and frequencies close to the texture's, thus yielding feature vectors allowing texture detection or classification.


Why is DTFT of $e^{jn\omega_0}$ an impulse train? update : After asking the question, I figured out that DTFT result is an impulse train. Now my question evolved to, how it is derived in this way?

Using the DTFT formula seems not to be working, because of "sum does not converge".

$$e^{j\Omega_{0}n}=2\pi \sum_{k=-\infty}^{\infty} \delta\left(\Omega-\Omega_0-2\pi k\right)$$

However, that DTFT is defines as above in most places.

Why cannot I find the transform using the forward DTFT formula?

First of all, note that the discrete-time Fourier transform (DTFT) of a sequence is always periodic with period $2\pi$:

$$X(\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\tag{1}$$

This is obvious because $e^{-jn\omega}$ is $2\pi$-periodic in $\omega$.

Next, consider the inverse DTFT:

$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)e^{jn\omega}d\omega\tag{2}$$

For $X(\omega)=\delta(\omega-\omega_0)$, we get from $(2)$

$$x[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}\delta(\omega-\omega_0)e^{jn\omega}d\omega=e^{jn\omega_0}\tag{3}$$

[Note that without loss of generality we assume that $-\pi\le\omega_0\le\pi$ because $X(\omega)$ is $2\pi$-periodic. $X(\omega)=\delta(\omega-\omega_0)$ is a shorthand for $X(\omega)=\sum_k\delta(\omega-\omega_0-2k\pi)$. That shorthand can be used here because we know that $X(\omega)$ is the DTFT of some sequence, and hence $2\pi$-periodicity is understood].

Since $\textrm{IDFT}\{\delta(\omega-\omega_0)\}=e^{jn\omega_0}$ we're inclined to accept

$$\textrm{DTFT}\{e^{jn\omega_0}\}=\sum_{n=-\infty}^{\infty}e^{jn\omega_0}e^{-jn\omega}=\delta(\omega-\omega_0)\tag{4}$$

(where $2\pi$-periodic continuation is understood), even though
the sum

$$\sum_{n=-\infty}^{\infty}e^{jn\omega_0}e^{-jn\omega}\tag{5}$$

doesn't converge in the conventional sense. However, this can be made mathematically sound by interpreting the sum $(5)$ as a distribution.


Subtracting DC Offset a Complex Signal? If I take the magnitude of a complex signal, and then subtract the mean amplitude, I remove DC offset but lose the real and imaginary components. How can I properly subtract the DC offset of a complex signal? Subtracting the complex amplitude results in a non-zero magnitude so I don't believe that's the way to go...

Using the term "DC offset" is a little misleading in this context.  I think what you are asking is how do your "center" your signal.  That is simple:  Take the average value (sum of N samples divided by N) and subtract it from each sample.  A "DC offset" is usually considered in terms of real valued signals.  The DC stands for "Direct Current".


Coherent Sampling And The Distribution Of Quantization Noise I have a question concerning the distribution of noise in the frequency domain in case of Coherent Sampling. I read up about Coherent Sampling and understood that in order for a frequency $f_{in}$ to be located on one single bin of the DFT, the period or a multiple of it must fit exactly into the data window of the DFT, meaning that the following must be true:

$$
\frac{f_{in}}{f_s}=\frac{N_c}{N}
$$
$N$ ... Record Length, $N_c$ ... Number of cycles in data window

Furthermore in the slides of a course about Mixed-Signal Processing I am attending it is mentioned that $N_c$ and $N$ need to be coprime in order for the quantization noise to be distributed uniformaly. Otherwise the quantization noise would show a discrete spectrum. Why is that the case?

This is most easily understood by examining the error between the true sinusoid and the quantized sinusoid. Let's start by defining a simple sinusoid. For simplicity, I will assume a cosine, but the same arguments apply for a sinusoid with an arbitrary phase offset. The true sinusoid is defined in continuous-time as

$$ x(t) = \cos\left(2 \pi f_{in} t\right). $$

Now, we need to sample this signal with a period of $T_s = \frac{1}{f_s}$ to obtain a discrete-valued (but unquantized) signal

$$ x_D(n T_s) = \cos\left(\frac{2 \pi f_{in} n}{f_s}\right).$$

Using the relationship that you provided in your question, we can rewrite this as

$$ x_D(n T_s) = \cos\left(\frac{2 \pi N_c n}{N}\right). $$

From this it should be clear that this signal lies in two bins, specifically the bins at $N_c$ and $N - N_c$ of an $N$ point FFT. However, this signal is not realizable with finite precision numbers. So, we need to quantize it. Let us assume that we are quantizing to $q+1$ bits using two's complement. A reasonable approximation to the quantized signal is given by

$$ x_Q(n T_s) = 2^{-q} \left\lfloor 2^{q} \cos\left(\frac{2 \pi N_c n}{N}\right) \right\rfloor.$$

Another way to represent this would be to represent the sinusoid added with an error term (i.e., the quantization noise $\epsilon_q(n T_s)$). This is captured in the following:

$$ x_{D}(n T_s) = x_Q(n T_s) + \epsilon_q(n T_s). $$

Now, let's consider what $\epsilon_q(n T_s)$ looks like for the two cases you mentioned. That is, the coprime case and the common factor case. Let's rearrange this slightly so that we can see things a little better,

$$ \epsilon_q(n T_s) = \cos\left(\frac{2 \pi N_c n}{N}\right) - 2^{-q} \left\lfloor 2^{q} \cos\left(\frac{2 \pi N_c n}{N}\right) \right\rfloor. $$

Now, consider the case where $N_c$ and $N$ share a common factor, $M$. Let's say $M$ is the greatest common factor. In this case, $N_c = M K_c$ and $N = M K$. Using this and the expression derived for $\epsilon_q(n T_s)$, we can show that

$$ \epsilon_q(n T_s) = \cos\left(\frac{2 \pi K_c n}{K}\right) - 2^{-q} \left\lfloor 2^{q} \cos\left(\frac{2 \pi K_c n}{K}\right) \right\rfloor, $$

and from this, we can see that

$$ \epsilon_q(n T_s) = \epsilon_q(n T_s + K T_s) $$

since

$$ \cos\left(\frac{2 \pi K_c n}{K}\right) = \cos\left(\frac{2 \pi K_c (n + K)}{K}\right).$$

In other words, the sequence is periodic with period $K < N$ since $K$ and $K_c$ are coprime. This means that all of the information in this signal can be captured with a $K$ point FFT. This, in turn, means that all of the energy from this signal is contained at the equivalent bins in the $N$ point FFT. These are the bins at $0, M, 2M, ... (K-1)M$. All of the other bins are zero-valued. This is why you see a discrete spectrum as shown in the slides you mentioned.

Now, if $N$ and $N_c$ are coprime, then the sequence is non-sinusoidal and aperiodic (or viewed another way, it is periodic with period $N$). There is no shorter period embedded in the sequence as there was with the common factor case. Therefore, all bins are required to represent the signal.


How to constrain the object detection problem? How to constrain the object detection problem?

It would be tempting to ask for precision such as: "It depends, does it have to be viewpoint invariant, is there going to be viewpoint changes, etc.?" Actually, it is exactly what I want to know. I have a specific application in mind and I need a structured way to help me constrain and define correctly the problem. I am pursuing a graduate degree but I am new to the field.

I mean, at first I read the literature about things such as binary keypoint descriptor, etc. but in the end, I am putting the cart before the horse, not having constrained my problem correctly, I cannot know or prove whether binary keypoint descriptors or else are a good starting point for my application...

Unfortunately my adviser is unavailable for a week or two and my colleagues are not specialized in that field.

A long time ago I have come across a very good advice: research should be problem-oriented, not approach oriented. That means find a problem you want to solve, and then look for appropriate methods.  Not the other way around. So, I would recommend starting with a specific problem, such as a real-world application.

For example, maybe you want to build a system with a camera looking at your bookshelf, which detects a particular book.  Or maybe you want to build a mobile robot with a camera, which can follow you around.  Or a video player program, which detects whether you are looking at your computer screen, and pauses the video when you look away (some cell phones do that now).  Or a program that can organize your photo collection by detecting and recognizing people, and maybe landmarks, like the Eiffel Tower.  The possibilities are endless.

So the algorithm is as follows:

Find a real-world vision problem that you find interesting.
Make as many restrictive assumptions as you can to scope it down and still keep it interesting.
Start reading papers to find applicable algorithms.
Understand the limitations of existing algorithms, and come up with your own ideas to address them.


Spectrum in real life with Spectrum in theory as you know it is possible to measure the power spectral density ($|X(f)|^2$) of a certain voltage signal $x(t)$ through a Spectrum Analyzer.
As far as I'm concerned, the simplest spectrum analyzer does not perform the Fourier transform X(f) of the input signal to evaluate the power spectral density. On the contrary, it simply selects each harmonic of the input signal through a tuned filter (or a tuned demodulator) and measures the power of the selected harmonics. The measured power is plotted against all the possible frequencies in the display.
This is what I've seen when I used a 2GHz sine wave as input signal.

It is fine because it is similar to an ideal dirac-delta what I would have expected from theory.
But, there is something I cannot understand. The Fourier transform $X(f)$ (and hence the power spectral density $|X(f)|^2$) is something that may be evaluated only by knowing the full signal behaviour on the infinite time interval $[-\infty;+\infty]$. And the resulting quantities $X(f)$ and $|X(f)|^2$ are obviously constant in time.

However, in real world we cannot know the future and we are sending to the Spectrum Analyzer a real time signal. And it does provide a real-time spectrum, which is not the result of a Fourier Transform evaluation (which wouldn't be possible as we do not know the future values of the signal) but of a simple power measurement.
So, my question is: why does this power measurement (which is an instant function of time: if I stop the input signal, the spectrum becomes 0) coincide with the theoretical power spectral density $|X(f)|^2$ (which is not function of time and which is a pure ideal quantity which requires we to know the signal behaviour also in the future)?

why does this power measurement (which is an instant function of time)

It's NOT a instant function of time.
For example instantaneous power of a sine wave is
$$ p(t) = x(t)^2 = A \cdot cos^2(t) = \frac{A^2}{2}(1+cos(2t))$$
You have to apply some amount of time averaging to get rid of the $cos(2t)$ component and that's what all spectral analyzer do.

: if I stop the input signal, the spectrum becomes 0) ...

Not instantaneously. Since there is time domain averaging involved the power measurement can only change at a finite rate.
In general there is an "uncertainty principle" at work here: the higher the spectral resolution, the slower the system can react or track. For your first implementation, that is determined by the bandwidth of your demodulation filter. The smaller the bandwidth, the slower the output will change.
If you use a Fourier Transform based analyzer, this will be determined by the length of the (finite) time domain window you apply the transform over.

why does this power measurement ... coincide with the theoretical power spectral density |X(f)|2

It doesn't. All physical measurements are approximate in nature. They all have a margin of error, tolerance, repeatability variance etc. You can't measure anything exactly, but only "good enough" for your specific application".
You also need to define what exactly you want to measure within the context of your specific application and requirements. You don't need the exact power spectrum of an infinitely long signal, since infinitely long signals do not exist in nature.


Log-Polar DFT Based Scale-Invariant Image Registration I'm trying to do image registration using phase correlation as described in the Reddy Chatterji paper. In my case, the images may be scaled and translated relative to each other.

The algorithm for finding the relative scale, as I understand it, is (see: the flowchart from the paper):

F1 = DFT(I1)
F2 = DFT(I2)
H1 = Highpass(F1)
H2 = Highpass(F2)
L1 = LogPolar(Magnitude(H1))
L2 = LogPolar(Magnitude(H2))
PC = PhaseCorrelate(L1,L2)
PM = norm(PC)
R = IDFT(PhaseCorr/PM)
P = Peak(R)
Scale = LogBase^P[1]

Scale gives me seemingly nonsensical values (wildly different from image to image and never correct).

But ignoring scale, the same phase correlation approach works fine for translation; and so I suspect I have a problem with my log-polar transform. Here's an example, where I've solved for translation--the left image is the original and the right has been cropped and translated--the solution is shown on top of the orignal:

For the log-polar transform, I first transform into polar space:
$$
\hat{I}(\rho,\theta) = I\left(r+\rho\cos\left(\frac{2\pi\theta}{N_{\theta}}\right),r-\rho\sin\left(\frac{2\pi\theta}{N_{\theta}}\right) \right)
$$
where $I$ is the original image, $r$ is the image radius (half-width) and $N_{\theta}$ is the number of samples in the $\theta$ direction. I then sample from this to transform into log polar space:
$$
\hat{I}_{log}(\rho,\theta) = \hat{I}\left(\log_{b}(\rho),\theta\right)
$$
where $b = (2r)^{-N_{\rho}}$ as described in 1 so that it spans the whole polar space.

Here're the example images in log-polar space with $\rho=\theta=256$ (in case there's something glaringly wrong):

Lastly, this shows the actual transformation the images go through before the phase correlation step (top is DFT magnitude post highpass filter, bottom is that in log polar space):

I'm using OpenCV, which has LogPolar and PhaseCorrelate methods. While the PhaseCorrelate, like my manual implementation, gives me the correct answer for translation, it's incorrect on scale. Since using the OpenCV LogPolar or my own doesn't affect the correctness, I must be missing something.

Any help would be appreciated.

If you want something really robust, but that might be more computationally expensive, you might want to check out the algorithm I implemented here. It implements the paper, "Robust Image Registration Using Log-Polar Transform" (pdf). It also has the advantage of being rotation-invariant, in addition to translation and scale invariant. In my application (art), it was able to register even even similar-looking images, not just transformed versions of the same image.


Generate audio data using VAE+GAN I am trying to train a VAE+GAN model to generate sounds produced by honeybees. I build my model by slightly modifying this tutorial, which aims to generate new MNIST images. Since my data are 1D signals, I apply STFT to the input to generate image-like features whose height is time and width is frequency. However, the total loss quickly goes to inf or nan just after a few steps(~10). Any idea why? What is the legit way of using VAE+GAN model to generate audio data?

If you are applying the Short Time Fourier Transform (STFT) , then your data are more likely to look like a spectrogram, where the "height" is frequency and the "width" is time.

The..."legit" way to train an autoencoder type of network so that it learns to generate periodic-like sounds is to "teach" a bank of oscillators to "fire" in a predefined pattern. This "pattern" is what the VAE will "learn" from the training dataset.

By applying the STFT, you are segmenting your time into $N$ bins, where $N$ depends on the length of your STFT "window". Consequently, the STFT "window" also determines your frequency resolution. If you choose a low $N$, you will get very good temporal resolution but very poor frequency resolution too and therefore, this is a design choice.

Once you apply the STFT, you get the equivalent of a spectrogram. Taken a bit more liberally, the spectrogram gives you this information. Only that, in the case of the spectrogram, the "keyboard" is rotated 90 degrees clockwise (in a spectrogram the vertical axis is frequency and high pitch is upwards, while the horizontal axis is time with positive time extending to the right).

Every one of the keys is an oscillator. If you scan the spectrum at a given "time window" $t_n$ within a sound file and you find a high amplitude in that frequency then the key is "pressed" (the oscillator is active). Of course, the other parameter is how hard is the key pressed (amplitude).

So, your auto-encoder has to learn to predict the sequence by which it has to activate the oscillators and at what amplitude, in order to generate the sounds your are after.

This is the bigger picture. The "detail" to look for here is that the output of the STFT is $\mathcal{C}$omplex and you really need it to be to be able to take into account the phase that these oscillators are firing with respect to each other. Otherwise, you will be getting a baseline error below which your model would not be able to optimise any more.

More generally, I think that you would benefit from a little bit of background on speech (or music) generation by looking into what a vocoder is and how does it work, Prony's method by which signals can be synthesized from damped sinusoids (This is the closest to the piano analogy), speech synthesis and of course, the STFT and spectrogram links from above.

Hope this helps


How to Calculate Local Spatial Variance of an Image? For those who work often on image processing, I'm facing an algorithm of Non-Uniformity Correction (NUC) that requires the calculation of a parameter called "Local Spatial Variance". I can't find a definition for it.
Does anyone have any idea about it ?
The article is called A Novel Infrared Focal Plane Non Uniformity Correction Method Based on Co Occurrence Filter and Adaptive Learning Rate, IEEE Access, 2019.
The "Local Spatial Variance" is on the 3rd page if you please.

Usually it means the variance within a window of pixels where the (i, j) pixel is in its center:

In the above you can see a 3 x 3 window (Radius of 1).
You'd calculate the variance in this window and set its value to the location of the Pij pixel.


how to make fft ignore phase shift when only interested in magnitude I have a signal that consists of a sine wave, on which I apply a rectangular window and then an FFT. I'm only interested in the magnitude of the sine wave, and I want it to stay the same even if the signal shifts a few samples to the right relative to the window (of course assuming the sine wave has the same amplitude from top to bottom). Is this possible and how can this be done best? The application is a kind of pitch detector where I have to detect a small group of different pitches.

You can't make an FFT ignore phase (for strictly real input), except for signals which have frequencies that are exactly periodic in the FFT window length.

Otherwise, you will have to interpolate to get any between bin magnitude peaks.  And an accurate magnitude interpolation of any non-bin-centered frequency will depend on phase (easiest if you reference phase to the center of the data window), especially near the DC or Fs/2 bins, but to a much smaller extent, thru the rest of the FFT result as well.


Algorithms for detecting the difference between two signals If I have two signals $X$ and $Y$, where $X$ is a clean signal and $Y$ is the same signal with linear white Gaussian noise and an amplitudes of $10^{-4}$. How could one use an algorithm to detect if there are similarities between the signals.

I have currently come across cross correlation, which I am still trying to understand but are there any other methods that could be use, the idea is that I would be able to code it in matlab, so no tool boxes can be used. Could there be a matrix approach where you read the signal into a matrix form and decompose the noise and amplitudes from the matrix to get back to the original signal

This is currently a project I am doing for one of my course, so I am not looking for someone to give a answer, I am looking for advice and where I should be look, maybe relevant documentation.

The two samples are 16-bit and sampled at 44.1Khz

My apologise I have just been informed that the Y signal was renormalized after the Gaussian white noise and amplitude modulation.

This is a pretty common problem, and you’ll see it all the time in radar/sonar/___ar processing. Cross-correlation is really the thing you’d want to use here in my opinion.

In layman’s terms, you can use cross-correlation to do “template matching”. In radar applications, were looking for the signal we sent out to return to us. Since we know what signal we sent out, we can use that as a template and use cross-correlation as a mathematical tool to “search” for that template signal in the received signal.

So in your problem statement, you have the signal X which and the signal Y. It seems that Y = X + N, where N is some noise signal. To use cross correlation in MATLAB, you can simply just do [xc, lags]=xcorr(Y, X). Now ideally, the peak of the magnitude of xc will tell you precisely where the signal X “starts” in your signal Y. It’s important to note that the output of the MATLAB xcorr function will include “negative” time lags, so when you plot it, use plot(lags, abs(xc)).

If you plotted the magnitude of xc, that is abs(xc), you’ll likely notice a discernible peak above all of the noise. This would tell you that indeed, the signal X is present in the signal Y.

Now that’s just a simple overview of cross-correlation; naturally there are a lot of fine details, but hopefully that should get you started!

Implementation note for cross-correlation:
If you wanted to do circular cross-correlation or not use any built ins, you could do it with discrete Fourier transforms (DFT/FFT):
xc = IFFT( FFT(Y) .* conj(FFT(X)) ). If you’d like to learn more about the mathematics of this, check out Wikipedia and the convolution theorem.


How to find a certain frequency in a signal that varies in time? I have an EEG signal. It is amplitude verses time. I want to find out how certain frequencies in that signal vary with time.

One way I think I can do this is to do a short-time-Fourier-transform on it (that is, split it in small pieces and do a Fourier transform on each piece), then track the frequency I am interested in across the slices.

Is there a different/standard way of doing this? I only ask because Its going to be a lot of processing overhead doing it my way and I'm doing this in python.

thanks in advance

cleerline

If I understood it right, you intend to figure out how some frequency content of a given signal evolves in time.

Next to short time Fast Fourier Transform, you have the opportunity to make use of a time-frequency analysis. I would suggest you to use either the Hilbert Transform (https://en.wikipedia.org/wiki/Hilbert_transform) or the Teager-Kaiser energy operator (http://folk.uio.no/eivindkv/ek-thesis-2003-05-12-final-2.pdf). Both methods are able to make estimation of istantaneous magnitude and frequency of a signal over time. I assume that your signal is multicomponent, therefore you should first isolate the frequency components you want to analyze by means e.g. of a filter-bank and then apply one of the two mentioned methods.
Otherwise, you could directly go for an Empirical Mode Decomposition (https://en.wikipedia.org/wiki/Hilbert%E2%80%93Huang_transform). Such a method will decompose the signal in terms of frequency components: after that you can again apply Hilbert or Teager-Kaiser methods to check how your target-frequency develops over time.


Doubt about wide sense stationary random process I have white Gaussian noise $F[n]$ with zero mean and autocorrelation $R_F[n_1,n_2]=\delta[n_1-n_2]$.

If now I consider the random process defined as

$$X[n]=u[n]e^{-kn}F[n]$$ Is $X[n]$ a wide-ense stationary random process for all values of $k$?

I tried to solve the autocorrelation function to check if it depended only on the lag but my result was that $X$ is WSS only for $k=0$, but I think that since it's a sort of a linear filtering it should be that it is WSS for all $k$.

The mean of $X(n)$ is always $\mu_X =0$, because the noise has zero mean. Thus we should check whether the autocorrelation corresponds to a WSS process.

$$\begin{align}
R_X(n_1,n_2)
&=\mathbb{E}[X(n_1)X(n_2)]\\
&=\mathbb{E}[u(n_1)e^{-kn_1}F(n_1)u(n_2)e^{-kn_2}F(n_2)]\\
&=u(n_1)u(n_2)e^{-k(n_1+n_2)}\mathbb{E}[F(n_1)F(n_2)]\\
&=u(n_1)u(n_2)e^{-k(n_1+n_2)}\delta(n_1-n_2)
\end{align}$$

We can evaluate the autocorrelation function in two different cases:

If $n_1-n_2\neq0$, then $$R_X(n_1,n_2)=0$$
If $n_1-n_2=0$, then $$R_X(n_1,n_2)=u(n_1)e^{-2kn_1}=u(n_2)e^{-2kn_2}$$

The second case shows that $R_X$ depends on the specific values of $n_1$ or $n_2$ and not just their difference, thus the process is not WSS.


How can I implement a binary chirp signal containing zeros and ones? I build a signal, $x$, which is linear in the frequency domain as (in MATLAB syntax)

x = 1:0.1:1e3;

Then I do:

y = (ifft(x));

and I would like to build a binary signal (zeros and ones) from $y$, but I can't define any threshold value such that

if y(index) > threshold
    y(index) = 1
else y(index) = 0

How can I implement a chirp of zeros and ones in the time domain?

What is a chirp signal?

A chirp signal is a sinusoidal signal of the form $\sin(bt+\frac{1}{2}at^2)$ whose instantaneous radian frequency (the derivative with respect to time of the argument of the sine, viz. $b+ at$, is a linear function of $t$. If the chirp signal is of
duration $T$, then the instantaneous frequency at the end is $B=b+aT$.
If $a > 0$, the frequency increases
with time from the initial value $b$ to its maximum value $B$ at $t = T$,
while if $a<0$, the frequency decreases from
$b$ downwards to the minimum value $B$ at $t = T$.

Is the Fourier transform of a chirp signal constant over
the frequency range $[b,B]$  (or $[B,b]$ for pessimists who choose $a < 0$)?

No, the Fourier transform is not constant over this
frequency range.  Nor does the DFT have equal amplitude
in all bins if we sample the chirp signal and compute its DFT.

What is a fake chirp signal?

A fake chirp signal is one whose instantaneous radian frequency
increases in steps of equal height $\Delta = (B-b)/n$ but decreasing
durations so that it consists of one full period of $\sin(bt)$
followed by one full period of $\sin((b+\Delta)t)$ followed by one
full period of $\sin((b+2\Delta)t)$, $\ldots$, followed by one full period
of $\sin((b+(n-1)\Delta)t) = \sin((B-\Delta)t)$.  Note that the total
duration is
$$T = 2\pi\left[\frac{1}{b} + \frac{1}{b+\Delta}
+ \cdots + \frac{1}{b+(n-1)\Delta}\right]$$
To maintain phase continuity, we insist that each full period start
and end at $0$ amplitude and $0$ phase so that, for example, between
$t=\frac{2\pi}{b}$ and $t=\frac{2\pi}{b}+\frac{2\pi}{b+\Delta}$ when
we have a sinusoid of radian frequency $b+\Delta$,
the exact mathematical expression
for the signal is $\sin((b+\Delta)t+\theta)$ where $\theta$ is chosen
so that $(b+\Delta)\frac{2\pi}{b}+\theta$ is an integer multiple of $2\pi$.

What is a binary chirp signal?

If you create a chirp signal, sample it, and hard-limit it to $\pm 1$
(or, as you prefer, shift the dc level so that the signal is $0$s and $1$s)
then you can call it a binary chirp signal, but you shouldn't
expect it to have the properties one expects of a chirp
signal.  In any case, it is necessary to have $b,B \ll f_s$
where $f_s$ is the sampling rate to get something reasonable.
If you don't want to go through the hassle of creating a sampled
chirp signal, etc. you can get something that can be called
a faker  binary chirp signal very simply.  Here are
two examples.

$$111110000011110000111000110010~\text{or}~1111100001110010$$

The first sequence is effectively one period of a sinusoid with $10$ samples
per period followed by one period of a sinusoid with $8$ samples
per period followed by one period of a sinusoid with $6$ samples per period,
etc.  The second sequence has half periods of sinusoids with decreasing
numbers of samples per period.  I call this a faker binary chirp
signal because instead of the frequency of the sinusoid
increasing in constant steps, it is the period of the sinusoid
that is decreasing in constant steps.  So we don't have the linear
frequency sweep that characterizes chirp signals.  But the signal is
trivial to synthesize, and if $b, B \ll f_s$, the difference between
a fake binary chirp and a faker binary chirp can be small.

I don't read MATLAB but if x = 1:0.1:1e3; means that it is desired
that the frequency sweep from $1$ Hz to $1000$ Hz (or by a factor of $1000$
in general), then the digital sequence (whether sampled from a true chirp
pr a fake chirp or synthesized as a faker chirp) is going to begin with
long strings of $1$s followed by long strings of $0$s.  Matters
will be even worse
if a frequency resolution or step size $\Delta$ of $0.1$ Hz is required.


Is Oversampling a Signal Same as Discretizing the Signal? I am confused about the new terminology 'oversampling' I have encountered very recently.
To explain it better, I am going to use numerical examples.
Now let's assume I am using a BPSK modulator and my symbol period is 1 second. To digitize and process this using MATLAB or Python, what I was doing was to discretize the signal with 100. (Using np.repeat() or repmat to information bit sequence).
My questions are:

1) Is this value (100) the oversampling factor?

2) Is this discretization value of mine, is the ratio between symbol duration times sampling frequency of the system ?

3) If the ratio that I mentioned earlier is not what oversampling is, can you explain or give me links for some kind of explanation? ( I found a lot about it but it made me even more confused due to the fact people are using upsampling and oversampling interchangeably, which I think is totally wrong)

Thanks in advance.

Oversampling is the case the rate the data is sampled is higher than required by the data Bandwidth and Nyquist Shannon Sampling Theorem.
It has many good properties regarding the processing yet it also generate some issues (Like coloring the quantization noise).

Regarding your questions:

In BPSK the symbol is the phase range of modulated signal (Sign). By using repmat() you indeed oversmapled the data (Symbols) yet not the modulated signal. What you needed to do is creating a grid where you sample the modulated signal. So if the the modulated signal has a frequency $ f $ than for 100 oversampling rate you should sample at $ {F}_{s} = 100 \cdot 2 \cdot f $.
What you did is replication of the symbol. Not higher sampling rate unless you imitate the effect of oversampling after Symbol inference.
In order to see the effect of oversmapling by factor of $ 100 $ generate the modulated signal over a grid where $ \Delta T = \frac{1}{100 \cdot 2 \cdot f} $ with $ f $ being the modulating signal frequency.

Regarding your question in the title.
The process of sampling data is basically discretizing data over 2 grids:

The Time Grid
Basically sampling the continuous data in specific times. Oversampling refers to the grid interval.
The Value Grid
Process we usually call quantization which means the value we sample isn't continuous but has discrete property of being a value withing a pre defined set of values (Quantization Grid).


A-GPS frequency (fd) search ranges A-GPS is able to reduce the fd search range from several kHz (GPS search range) to several hundred Hz (a-GPS search range). I know the GPS search range has a center at 0 Hz. What is the center frequency with the reduced A-GPS search range? I got an impression that this A-GPS frequency search center shouldn't be located at 0 Hz. Instead, it should be the assumed initial fd value.

A-GPS will give you ephemeris data. That will be used to update your initial guess for the frequency of the downlink signal, indeed. The relevant remaining uncertainty  is in the local oscillator of the GPS receiver itself. And that can be large compared to Dopplers!

So, yes, assuming the frequency probability density to be symmetrical convex around the nominal frequency, that would put an ML estimator at the doppler frequency.

But: your GPS receiver typically doesn't run the first time. It already has acquired estimates about its own frequency error in the past, and together with parameters like current temperature and the time that has passed since the last successful fix, your initial frequency error estimate doesn't have to be 0; for example, if I remember correctly (it's really a bit of a fuzzy memory, please check the source code) RTKlib uses an extended Kalman predictor for such purposes.

Also, the highest likelihood estimate for an initial frequency error is not necessarily one that makes for the best starting point in a search algorithm!

For example: Due to the ambiguities in multiples of (chip rate)/(chip sequence length), it might be cleverer to first test a frequency hypothesis that's intentionally a bit off, because it can give a drastically improved expected convergence rate when looking for a frequency error optimum – in the end, this is a relatively complex optimization problem that's not globally convex.


How unauthorized access is prevented in biometric identification systems biometric access control is a rapidly maturing technology that is used to police our ability to gain access to places, allow us access to software and hardware as well as guard against unauthorized access.
If the number of people allowed to access a certain a place is for, example, 100. Training samples from the 100 persons are used to construct the stored database. Features are extracted from these training samples, and are stored in the system database. K-nearest classifier can be used in the test stage to identify the person, by selecting the person with the minimum distance between the test features and the stored features.
I can understand that this identification process will correctly identify the person if he is actually one of the 100 persons allowed to enter the place. I can not understand how the system will prevent unauthorized access. If the person is not one of the 100 persons allowed, of course the system will identify him as the one with the closest distance to the features stored in the database.

It's not the closest distance alone which lets you make a decision.
Given a pattern recognition (or object classification) system that's trained with N objects, then for a given test sample you will do the following:
Compute the metrics based on the error produced by the chosen features of the test sample and the database.
Then seek through the classes for which this error-metric is a minimum.
Then further check if this minimum error is acceptable (i.e., below a suitable threshold).
If yes, then you decide that you have identified the test sample.
Otherwise, you were not able to decide; a.k.a. unidentified access.