Why do we get different imaginary parts of a zero centered Gaussian for the the same number of data points N? Suppose we have a total number N= 2048 points in a data and we wish to have zero centered Gaussian. There are two possibilities that we use the x-axis as

x1=[-1023:1:1024]; % x axis spans from -1023 to 1024 with 1 unit steps

x2=[-1024:1:1023]; % x axis spans from -1024 to 1023 with 1 unit steps

and if we make two zero centered Gaussians using these x values:

Gauss1=normpdf(x1,0, 10); % The syntax is normpdf(x, mean, standard deviation)

Gauss2=normpdf(x2,0,10);

and obtain their FTs as follows in MATLAB.

FFTGauss1=fft(Gauss1);
FFTGauss2=fft(Gauss2);

The real parts are identical and their magnitudes exactly match. For some reasons, the imaginary parts vary drastically. Why do we see large imaginary parts in one case and almost non-existent imaginary parts in the other? Thanks.

Answer : When $x_2 = [-1024:1:1023]$, then $x_2[n]$ satisfies the condition $x_2[n] = x_2[(N-n)\mod N]$. That is why when $x_2 = [-1024:1:1023]$, then the FFT is real and hence the imaginary part is $0$. If you see the scale of the $y$-axis for imaginary part of $x_2$ plot, it is of the order of $10^{-17}$ which is almost $0$ in MATLAB.

Detailed Explanation:

When $x \in \{-1024, -1023, -1022,..., 0, 1, 2, ..., 1023\}$, then you get the following mapping :
$$\begin{array}{lcl}x[0] &=& {\rm gaussian}(-1024)\\
x[1] &=& {\rm gaussian}(-1023)\\
x[2] &=& {\rm gaussian}(-1022) \\
& \vdots\\
x[1024] &=& {\rm gaussian}(0)\\
x[1025] &=& {\rm gaussian}(1) \\
& \vdots \\
x[2047] &=& {\rm gaussian}(1023)\end{array}$$ Observe that $x[1] = x[2047] = x[(2048 - 1)\mod\ 2048]$, $x[2] = x[2046] = x[(2048 - 2)\mod \ 2048]$ and so on. This makes $x[n]$ real and symmetric mod $N$, which will in turn make $X[k]$ real and that is why you see that the imaginary part as $0$. MATLAB shows $0$ as values of the order of $10^{-17}$.

Do the same mapping for $x = [-1023:1:1024]$, and you will see that $x[n] \ne x[(N-n)\mod \ N]$ and hence imaginary part is not $0$.


Speech Analysis - Zero-Crossing Basically, (I know I've asked a simular question, but I have been away and really thought about it this time)!

Let me try and explain this in the best way possible:

Basically, I'm attempting to create an algorithm that takes a voice sample (.wav)
of someone saying either "Yes", or "No".  Now, in order to do this I am attempting to split the sample into equal blocks in order to get the Phones.

Now, I have thought about this..

I don't need to technically analise all of Phones, I can identify if the the sample is either "Yes" or "No" by just looking for a "Y" Phone. Because if a "Y" doesn't exist then the sample is "No". So, I can handle this by implementing a decision tree.

Ok, so to my question. I know that this can be achieved using Zero-Crossing, but, I don't have any experience in this (Ok, people are going to tell me to stop doing it then) but, we all had to learn somewhere, right?

If someone could possibly link me to a website or a book where I can learn about the Zero-crossing method? Then it would be great.

The main thing I don't understand about it is, how can Zero-crossing determine whether the Phone is a "Y"? Or, will I need to have some samples of "Y" in order to match through Correlation?, or something.

I hope someone can advise me.

Thanks :)

Zero-crossing rate (ZCR) might be useful for voiced/unvoiced frame discrimination, speech/music discrimination, but it is of much lesser importance in speech recognition. One reason is that it is pitch-dependent and not robust to background noise or hum. It is not difficult to craft very different signals (say a female voice saying a phone and a male voice saying another phone) that have the exact same ZCR; or to add some noise to a signal in a way that will half or double its ZCR.

Using zero-crossing rate might work in the single-speaker + prior training case. Collect a handful of samples of each words, compute the zero crossing rate, and train a naive bayesian classifier for this two-class problem. Be prepared to poor performances, and don't expect your system to generalize to other speakers.

The only plausible reason one would use such a feature would be on embedded systems with very limited computational power. Otherwise, there's no reason not to use "proven" features like MFCC or AR coefficients.


What could be the discrepancies between theory and practice when implementing a proportional controller? I implemented the following proportional controller circuit in practice.

The transfer function of the plant I uses is
$\frac{20}{(s+2)(s+3)}$.
Before doing the experiment, I plotted the root locus plot of the system.
Then after doing the experiment from the plots, the damping ratio was found using the 2nd order equations. During the experiments the gain of the circuit was changed such that different responses where obtained such as an under damped response.  These were found by the finding the value of the peak overshoot from the plots.  After finding the damping ratio, the new poles of the system were found and they were compared to the root locus obtained before.
They did not match exactly.  I know that there could be practical errors such as mismatches in the resistor values in the ciruit, error when finding the overshoot amplitude from the graph etc,   but could there also be some theoretical errors when using this system.

Since I am using the general equation of the second order system, could this have an effect, I know that the plant is second order but when considering the whole closed loop system, does this still apply?

There are no theoretical errors as you are simply describing a mathematical equation exactly. Also there are no actual 2nd order systems in the analog world due to real-world bandwidth limitations ("higher frequency" poles that we can typically ignore due to the dominant effects of lower frequency poles resulting in a reasonable estimate), and especially as your loop bandwidth increases, the effects of physical delays ($e^{-sT}$) in your system (which if you do the Taylor series expansion on that $e^{-sT}$ equation you see that it has an infinite number of poles!).

I found the following graphic useful in helping to understand the impact of the pole locations on the overall response (in particular to the point made above about "lower frequency" poles being dominant; the poles aren't actually on the frequency axis so the more accurate explanation is poles closer to the frequency axis are dominant where the frequency axis is the unit circle in the z-plane, or the imaginary axis in the s-plane as drawn below). Just look at a simple polynomial equation with two poles and you will see why it is "dominant" and how close a pole that is further out can be where it starts to matter based on the kinds of differences you are actually seeing.

When we do have a system with two complex poles which are sufficiently closer to the frequency axis than any other ones, we can make the simplifying and well documented approximations of a second order system. These are approximations.

For the graphic above which is an example Root Locus (so the poles move on the red trajectories shown as we change loop gain) I have used the following approximations:

Case 1: When the pole on the real axis is greater than 10x (in the negative direction) the real value of the complex poles, the system is "second order dominant" and it can be reasonably approximated by the two poles alone.

Case 2: When the pole on the real axis is less than 10x (in the negative direction) the real value of the complex poles, the system is "first order dominant" and it can be reasonably approximated by the one real pole alone.

Case 3: For values in between I expect to have greater errors in either approximation and likely need to model it as a third order system (depending on the accuracy of the results needed).

Here's an example from James Kamman at Western University that also used a 10x criteria but shows the significance of that in more detail:

https://web.archive.org/web/20150513061828/http://homepages.wmich.edu/~kamman/ME360DominantPoles.pdf


How Last Edge Image Can Be Achieved from Law Masks as we know thechnique of law edge detection has 25 2d masks that is obtained form 5 different 1d filters. we can use each of those masks on image with convolution but how can we obtain final result? is that achieved from from some kind of fusion or it have to be chosen?

thank you

If I understand correctly, the question is, given many images which are result of different Edge Filter applied on the same image, how to actually mark edges.

Well, you basically created 25 tests for each pixel to decide whether or not it is an edge.

You could apply many approaches to decide:

Majority Votes - If more than half of the voters decided it is an edge, nark it as an edge.
Threshold - If more than $ x $ voters vote for edge the pixel will be declared as an edge.
Spatial Model - Instead of per pixel decision, look around it and other voters.
Weight of Votes - Don't mark votes as "Yes / No" but give it a scalar. If the sum of all scalars above a threshold value, declare an edge.

As you can see the options are endless.


Statistical significance of coherence values I have two signals sampled at $100\textrm{ Hz}$ during $3\textrm{ seconds}$. So I have $300$ sample points at each signal series.
I want to calculate coherence values between these two time series and would like to know whether this coherence value is statistically significant. I found two resources online.

In the first resource, there is a parameter $n$ (it is the sample size as I could interpret). A significant coherence value for that $n$ is:
$$c = \sqrt{1-\alpha^{1/(n-1)}}$$

and this formula does not include any frequency component, unless $n$ is related to frequency instead of sample size.
In the second resource, the function returns degree of freedom and significance level at each frequency. However I am not sure whether this function returns the correct computation or not. Besides I couldn't interpret the results of this function.

My question: How can I find the significance of my computed coherence value at a certain frequency (or frequency band) preferably using MATLAB in the light of the above resources?

Great question! I've been grappling to determine how to calculate statistical significance of coherence for several days. This question helped me find some key references I hadn't previously located!

One detail regards whether you're talking about "coherence" or "magnitude-squared coherence". The first resource uses the former and the second uses the latter. I generally prefer the latter:

$$C_{x,y}(\omega)= \frac{|\langle{X^*(\omega)Y(\omega)}\rangle|^2}{\langle|X(\omega)|\rangle^2\langle|Y(\omega)|\rangle^2}$$

Where $X(\omega)$ is the fourier transform of $x(t)$, and $\langle \rangle$ denotes an ensemble average of $n$ independent ensembles.  Obviously the plain "coherence" is the square root of what is defined above.

The significance level of this magnitude-squared coherence as defined by the two resources is:

$$L_1(\alpha, q) = 1-\alpha^{1/q}$$
$$L_2(\alpha, q) = F_{2,2q}(\alpha)/(F_{2,2q}(\alpha)+q) $$

Where $\mathrm{DOF} = 2n = 2 ( q + 1 )$, and $F$ is the inverse-CDF (percent point function) of the F-statistic. Also, just FYI, the second formulation can also be found in Priestley 1981 page 706.

The following Python code implements the two approaches for $\alpha=0.05$ (i.e. a 95% significance level):

import numpy as np
import scipy.stats as st

alpha = 0.05
p = 1 - alpha

dof = np.arange(5, 100)
n = dof / 2.

C_1 = 1 - alpha ** (1. / (n - 1))

fval = st.f.ppf(p, 2, dof - 2)
C_2 = fval / (n - 1. + fval)

Plotting the two functions shows that they are apparently identical:

import matplotlib.pyplot as plt
ax = plt.gca()

ax.plot(dof, C_1, 'k-', lw=10, label=r'$1-\alpha^{1/q}$')
ax.plot(dof, C_2, 'r-', lw=5, label=r'$F_{2,2q}(\alpha)/(F_{2,2q}(\alpha)+q)$')
ax.plot(dof, 6. / dof, 'b-', lw=2, label=r'$6/d$')

ax.set_xlabel('Degrees of Freedom')
ax.set_ylabel('{:0.0f}% Significance Level'.format(100 * p))
ax.legend()
ax.set_title(r"Coherence level for $\alpha=0.05$")

I compared $L_1$ (black) and $L_2$ (red) for several different values of $\alpha$ and they always agree (their difference was always $O(10^{-16})$). I didn't bother doing the algebra to show that $L_1=L_2=L$, but perhaps someone else will find that to be an interesting undertaking?

Lastly, I've included the approximation: $L(\alpha=0.05)\approx6/\mathrm{DOF}$ (blue). I got this from U. Washington Applied Physics Lab professor Jim Thomson, which he apparently got from legendary oceanographer Carl Wunsch. It appears to be a very good approximation when $\alpha=0.05$ and $\mathrm{DOF} > 20$ or so. I'd be interested to see the algebra that shows that:

$$\lim_{\mathrm{q}\rightarrow\infty}L(1/20,q) = \frac{6}{\mathrm{DOF}}$$


MLE parameter estimation -- confusion regarding some terms in the pdf of complex normal r.v (Part 2) This question is based on the application of the pdf which was an earlier question of mine asked here Confusion regarding pdf of circularly symmetric complex gaussian rv

If $v \sim CN(0,2\sigma^2_v)$ is a circularly complex Gaussian random variable which acts as the measurement noise in this model
$$y_n = A + v_n \tag{1} $$ where $y$ is the observation and $A$ is a scalar unknown value which needs to be estimated. I am having a slight confusion whether there will be a 2 in the denominator of Eq(3) and Eq(4) with the $\exp(.)$ term. Based on the answer in the link, there should be no sqrt term with $\pi$ in the denominator, if $v \sim CN(0,2\sigma^2_v)$. If $v \sim N(0,\sigma^2_v)$ then there is a sqrt term.

Can somebody please check if I have correctly written  out the log-likelihood? I think I am missing a 2 in the denominator of $\exp[.]$ term in Eq(3) but I am not quite sure.

Thank you for your time and help.

$$P_y(y_1,y_2,...,y_N) = \prod_{n=1}^N\frac{1}{2\pi \sigma^2_v} \exp \bigg(\frac{-{({y_n-A})}^H ({y_n-A})}{2\sigma^2_v} \bigg) \tag{2}$$

taking log
$$\ell   = -N\ln(2\pi\sigma^2_v) -  \frac{1}{\sigma^2_v} {\bigg[{[\sum_{n=1}^{N} {(y_n - A)}{(y_n - A)}^{\mathsf{H}}  ]}\bigg]}. \tag{3}$$
$$ = -N\ln(2\pi\sigma^2_v)- \frac{1}{2 \sigma^2_v}{\bigg[ \sum_{n=1}^{N}y_n y_n^\mathsf{H}
- 2 \sum_{n=1}^N y_n A\bigg]}
  - \frac{1}{2 \sigma^2_v}{\bigg[  \sum_{n=1}^N   {AA}^\mathsf{H} ] \bigg]} \tag{4}$$

Hi: It's worded more clearly now in that you're estimating A and there's only one RV which makes more sense. Consider the first link I sent in the other message and go to where it says the "likelihood is". Your likelihood is quite similar except, since your variance is multiplied by 2, this causes 2 things to happen in terms of how would you change the likelihood in that document accordingly for your problem.

A) The $2 \pi \sigma^2$ that is getting raised to the $-\frac{n}{2}$ instead becomes  $4 \pi \sigma^2_{v}$

B) The denominator of the exp instead of  being $2 \sigma^2$  becomes $4 \sigma^2_{v}$.

Make those 2 changes to that likelihood and you will have your expression. I get confused looking at your 2) because my memory of 1/2's and the squares ( in your case H's ) is blurry and it can be confusing. So better to rely on document where likelihood is correct so you can't go wrong. Keep in mind that this is under the assumption that complex RV's don't change much compared to real ones.


How Can I Resample a Signal with an Arbitrary Factor (For Example - 128000Hz to 16000.1Hz) in MATLAB? I need to simulate the sampling of a continuous (fsCtu = 128000Hz), acoustic signal with two microphones that have a slight offset in sampling rate (fsMic1 = 16000, fsMic2 = 16000.1) in Matlab. What is the best way to do this?

Things I tried:

The Matlab "resample" command only works for resampling to 16000, not to 16000.1
"interp1" doesn't seem to be an option because I think I need to use bandlimited interpolation for a correct simulation. (Is this assumption correct?)
I tried to write my signal to a wav file and resample it via a system call using this software, then load the processed file with wavread. I'm not sure if this is a good solution. A quick test revealed that this method doesn't give the same result as the "resample" command for resampling to 16000Hz, which I find strange.

Any ideas or suggestions?

There are a couple of ways that you can do it.  The first is with resample, but it is a multi-step process.  First, you have to figure out which interpolation and decimation factors will get you the sample rate you want.

[n, k] = rat(16000.1 / 128000);

That gets you an interpolation factor of 20000 and a decimation factor of 159999.  You factor those to break them up into smaller chunks.

nFactors = factor(n)
kFactors = factor(k)

It turns out that n factors to $2^5 * 5^4$ and k factors to $3 * 7 * 19 * 401$.  All of those are doable, though the decimation by 401 will not have great filtering properties.  Anyway, if you resample in stages you can get the final sample rate you want.

The other way to do it is polynomial interpolation.  Essentially you model your signal as a polynomial through curve fitting techniques, and then you can simply feed in the time values that you want and out will pop the signal values.  This technique can be very effective, but modelling the signal well is a bit of an art.  In particular you don't want to overfit.


How to modify zeros and poles in a delta-sigma modulator loop? In this specific case, I am modeling a delta-sigma analog-to-digital converter with a basic loop in which the output $V$ is fed back with unity gain to be subtracted from the previously sampled input $U$ (discrete time) and this difference is fed to an integrator modeled with the transfer function:
$$G(z) = \frac{z^{-1}}{1- z^{-1}}$$

Quantization noise is added in form of error $E$ to the output of the integrator, giving the output $V$. It follows:
$$V(z)=\frac{G(z)}{1+G(z)}U(z)+\frac{1}{1+G(z)}E(z)$$
The complete model, however, should be a third-order system in which the output $V$ is fed back with unity gain twice more so that it is once subtracted from the sampled input in the outer loop, this difference is then integrated and fed to the next loop, in which the fed-back output is subtracted again. This new delta enters the second integrator and the same process repeats once more in the inner loop, where the quantization error is added as described in the first lines above.
In total there will be 2 identical outer loops and a third inner loop where the quantization is added. In total there will be 3 integrators with transfer function $G$.

the minus sign is missing in the second picture for all three summation nodes. it was a mistake. the feedbacks should be indeed negative
Now I'm interested in the Noise Transfer Function: NTF = $V / E$. According to my computations:

NTF with gain unity in all three feedback paths (alpha, beta and gamma are 1) should be equal to:
$$\text{NTF}(z)=\frac{V(z)}{E(z)}=\frac{1}{G(z)^3+3G(z)^2+3G(z)+1}$$
It follows:
$$\text{NTF}(z)=\frac{(z-1)^3}{(z-1)^3+3(z-1)^2+3(z-1)+1}$$
From here, I would like to do two things:

Replace the poles $z1 = 1$, $z2 = 1$ and $z3 = 1$ of the NTF with different values. In order to do so, I simply introduced a gain factor in the feedback path of each of the three loops (alpha, beta and gamma as above). This will add three coefficients in the denominator that I can play with.

However:

I would also like to replace two of the zeros of the NTF, which are now all equal to 1, with a pair of complex conjugate zeros. This should provide "lower noise around the zeros, at the expense of a less effective shaping at $z = 1$, thus increasing the usable frequency range".

If I understand correctly, the numerator $N(z)$ of the NTF should be something like this in the end:
$$N(z)=(z-1)(z-e^{jω})(z-e^{-jω})$$
My question is: how do I achieve this? What should I add in the model to shift the zeros slightly away from DC?
Just for information: once I have also obtained a pair of complex conjugate zeros, I should solve an optimization problem to find the optimal position of poles and zeros in the NTF for a determinate quantization noise power. So my humble guess is that ideally I should add something as simple as a gain in the model for the zeros, too, so that can adjust the numerator of the NTF, too, and solve for the optimal zeros within certain boundaries.

i dunno, i might have posted this before.  but here is a 2nd-order sigma-delta quantizer that operates not as an oversampled rate.  so you can hear the quantization noise, but you can also hear the music beneath it.  and the output is binary: $y[n] = \pm A$.
your professor is mistaken if he/she thinks they can model this without the $z^{-1}$ in the feedback path nor leave out the gain of the comparator.
%
%
%
%
%
%
%   simulated 1 bit sigma-delta converter:
%
%
%            x(n)-y(n-1)    w(n)                v(n)                 ( mean(y^2) = A^2 )
%
%   x ---->(+)--->[1/(z-1)]--->(+)--->[1/(z-1)]--->[Quantizer]----.---> y = +/- A = quantized value
%           ^                   ^                                 |
%           |                   |                                 |
%           |                   '----[-fbg]<----.                 |
%           |                                   |                 |
%           '------[-1]<------------------------'------[1/z]<-----'
%
%
%
%
%
%   "linearized" model:
%                                                          .---- q = quantization noise  ( mean(q) = 0 )
%                                                          |
%                                                          |
%            x - y/z        w                   v          |         ( mean(y^2) = G^2*mean(v^2) + mean(q^2) )
%                                                          v
%   x ---->(+)--->[1/(z-1)]--->(+)--->[1/(z-1)]--->[G]--->(+)-----.---> y = G*v + q
%           ^                   ^                                 |
%           |                   |                                 |
%           |                   '----[-fbg]<----.                 |
%           |                                   |                 |
%           '------[-1]<------------------------'------[1/z]<-----'
%
%
%
%
%
%
%
%
%           W = 1/(z-1)*(X - Y/z)
%
%
%           V = 1/(z-1)*(W - fbg*Y/z)
%
%             = (X - Y/z - fbg*Y*(z-1)/z)/(z-1)^2
%
%             = (X*z - Y*(1+fbg*(z-1))) / (z*(z-1)^2)
%
%
%           Y = G*V + Q = G*(X*z - Y*(1+fbg*(z-1)))/(z*(z-1)^2) + Q
%
%             = G*X/(z-1)^2 - G*Y*(1+fbg*(z-1))/(z*(z-1)^2) + Q
%
%
%           Y + G*Y*(1-fbg + fbg*z)/(z*(z-1)^2) = G*X/(z-1)^2 + Q
%
%
%           Y = (G*X/(z-1)^2 + Q)/(1 + G*(1-fbg + fbg*z)/(z*(z-1)^2))
%
%             = (G*X/(z-1)^2 + Q)*(z*(z-1)^2)/((z*(z-1)^2) + G*(1-fbg + fbg*z))
%
%             = z*(G*X + Q*(z-1)^2)/(z^3 - 2*z^2 + (G*fbg+1)*z + G*(1-fbg))
%
%             = z*(G*X + Q*(z-1)^2)/(z*(z-1)^2 + G*fbg*z + G*(1-fbga))
%
%
%    as z -> 1  (DC)
%
%           Y  ->  z*X/(fbg*z + (1-fbg)) =  X/(fbg + (1-fbg)/z)  -->  X
%
%

if ~exist('mean_vv', 'var')
    linearized_model = 0                % run this with 0 the first time to define G and mean(q^2)
end

if ~exist('A', 'var')
    A = 1.0                             % comparator output magnitude
end

if ~exist('fbg', 'var')
    fbg = 2.0                           % feedback gain to internal integrator
end

%
%   if there is an input soundfile specified, use it.  else, create a sin wave
%

if exist('inputFile', 'var')

    [inputBuffer, Fs] = audioread(inputFile);

    fileSize = length(inputBuffer);

    numSamples = 2.^(ceil(log2(fileSize(1))));  % round up to nearest power of 2

    x = zeros(numSamples, 1);                   % zero pad if necessary

    x(1:fileSize) = inputBuffer(:,1);           % if multi-channel, use left channel only

    clear inputBuffer;                          % free this memory
    clear fileSize;

    t = linspace(0.0, (numSamples-1)/Fs, numSamples);   % time

else

    if ~exist('numSamples', 'var')
        numSamples = 65536                              % number of samples in simulation
    end

    if ~exist('Fs', 'var')
        Fs = 44100                                      % (oversampled) sample rate
    end

    if ~exist('f0', 'var')
        f0 = 261.6255653                                % input freq (middle C)
    end

    if ~exist('Amplitude', 'var')
        Amplitude = 0.25                                % input amplitude
    end

    t = linspace(0.0, (numSamples-1)/Fs, numSamples);   % time
    x = Amplitude*cos(2*pi*f0*t);                       % the input

end

sound(x, Fs);                                   % listen to input sound
pause;

y = zeros(1, numSamples);                       % the output (created and initialized for speed later)

if linearized_model
                                                % artificial quantization noise for linearized model
                                                % mean(q) = 0, var(q) = mean(q^2) = mean(y^2) - G^2*mean(v^2)
                                                % does not have to be uniform or triangle p.d.f.
    q = sqrt(6.0*(A^2 - G^2*mean_vv))*( rand(1, numSamples) - rand(1, numSamples) );
else
    q = zeros(1, numSamples);
end

sum_yv = 0.0;
sum_vv = 0.0;

w = 0;
v = 0;
for n = 1:numSamples

    if linearized_model

        y(n) = G*v + q(n);                      % here the comparator is modelled as a little gain with additive noise

    else

        if (v >= 0)                             % the comparator
            y(n) = +A;
        else
            y(n) = -A;
        end

        q(n) = y(n) - (sum_vv+1e-20)/(sum_yv+1e-20)*v;

    end

    sum_yv = sum_yv + y(n)*v;                   % collect some statistics on v
    sum_vv = sum_vv +    v*v;

    v = v + w  - fbg*y(n);                      % second integrator
    w = w + x(n) - y(n);                        % first integrator

end

if ~linearized_model                            % don't recalculate this if using the linearized model
    mean_yv = sum_yv/numSamples;
    mean_vv = sum_vv/numSamples;
    G = mean_yv/mean_vv;                        % the apparent comparator gain (assuming stationary input)
end

%
%
%
%     Y = ((G*z)*X + (z^3 - 2*z^2 + z)*Q) / (z^3 - 2*z^2 + (G*a+1)*z + G*(1-a))
%
%
%
Hx = freqz([0  0 G 0], [1 -2 G*fbg+1 G*(1-fbg)], numSamples/2);
Hq = freqz([1 -2 1 0], [1 -2 G*fbg+1 G*(1-fbg)], numSamples/2);

plot(t, y, 'b');
sound(y, Fs);                                   % this could sound pretty bad
pause;

Y = fft(fftshift(y .* kaiser(numSamples, 5.0)'));
Q = fft(fftshift(q .* kaiser(numSamples, 5.0)'));

f = linspace(0.0, (numSamples/2-1)/numSamples*Fs, numSamples/2);

plot(f, 20*log10(abs(Y(1:numSamples/2)) + 1e-10), 'b');
hold on;
plot(f, 20*log10(abs(Q(1:numSamples/2)) + 1e-10), 'r');
plot(f, 20*log10(abs(Hq) + 1e-10), 'g');
axis([0 Fs/2 -50 100]);
hold off;
pause;

semilogx(f(2:numSamples/2), 20*log10(abs(Y(2:numSamples/2)) + 1e-10), 'b');
hold on;
semilogx(f(2:numSamples/2), 20*log10(abs(Q(2:numSamples/2)) + 1e-10), 'r');
semilogx(f(2:numSamples/2), 20*log10(abs(Hq(2:numSamples/2)) + 1e-10), 'g');
axis([Fs/numSamples Fs/2 -50 100]);
hold off;
pause;

semilogx(f(2:numSamples/2), 20*log10(abs(Y(2:numSamples/2)) + 1e-10), 'b');
hold on;
semilogx(f(2:numSamples/2), 20*log10(abs(Hq(2:numSamples/2)) + 1e-10), 'r');
semilogx(f(2:numSamples/2), 20*log10(abs(Hx(2:numSamples/2)) + 1e-10), 'g');
axis([Fs/numSamples Fs/2 -50 110]);
hold off;


How to get output from High Pass filter on this signal? I'm trying to solve this exercise:

  Consider a High-Pass filter with a cutoff frequency of $W=2\pi$

  Calculate the output signal if the input signal is $x(t)=\displaystyle \frac{\sin(4\pi t)}{\pi t}$

The fourier transform for the input signal is:

$X(\omega)= \left\{ \begin{array}{lcc}
             1 &   \text{if}  & -4\pi \lt \omega \lt 4\pi \\[2ex]
             0 &  \text{otherwis}e \\
             \end{array}
   \right.$

Now, the output signal is :
$Y(\omega)=H(\omega)X(\omega)$

$Y(\omega)= \left\{ \begin{array}{lcc}
             1 &   \text{if}  & 2\pi \le |\omega| \le 4\pi \\[2ex]0 &  \text{otherwise} \\
             \end{array}
   \right.$

I can get here, but the solution to the exercise converts this function with inverse fourier transform, giving the result of $x(t)=2\displaystyle \frac{\sin(\pi t)}{\pi t}\cos(3\pi t)$.

I don't know how did they get to that conclusion, I have tried to search for it, but i don't know how to look for that.

I hope that the question is clear. Sorry if I'm lacking some details, it's my first time asking this type of questions.

What you get as $x(\omega)$ is two $\mathrm{sinc}(\pi t) = \displaystyle\frac{\sin(\pi t)}{\pi t}$ shifted by $+\omega_0$ and $-\omega_0$ in frequency ($\omega_0 = 3 \pi$). By the frequency shift property of the FT:

$$
\mathcal F\left\{x(t) \cdot e^{i\omega_0 t}\right\} = X\left(e^{j(\omega - \omega_0)}\right)
$$

Then, remember $e^{i\omega t} = \cos(\omega t) + i \cdot \sin(\omega t)$. Therefore:

$$
x(\omega) = \sqcap_{2\pi} (\omega - \omega_0) + \sqcap_{2\pi} (\omega + \omega_0)
$$

Doing inverse:
$$
\begin{align}
x(t) &= \mathcal F^{-1}\left\{x(\omega)\right\} \\
 &=\mathrm{sinc} (\pi t) \cdot \big(\cos(\omega_0 t) + i \cdot \sin(\omega_0 t)\big) +  \mathrm{sinc} (\pi t) \cdot \big(\cos(-\omega_0 t) + i \cdot \sin(-\omega_0 t)\big) \\
&= 2 \cdot \mathrm{sinc} (\pi t)\cdot \cos(\omega_0 t) \\
&= 2 \cdot \frac{\sin(\pi t)}{\pi t}\cdot \cos(3 \pi t)
\end{align}
$$

since $\cos(-x)=\cos(x)$ and $\sin(-x)=-\sin(x)$


Mystery signal - FSK, BPSK? I have a mystery audio signal that I need to decode, it's from a discontinued product. Inspecting the object code shows a module called decodeFSK so that suggests it is a form of FSK. I've tweaked it in MATLAB ad-nauseum but think I must be missing something.

The only documentation I can find says:

  "[Binary data is PCM encoded PPG sensor waveform. 68bits of binary data are sent in 12msec bursts at 25Hz.] Current modulation scheme used is Manchester encoder/decoder. Modulated signal is a square signal with frequency range from 4kHz to 8kHz."

Not sure if the 4-8kHz part is correct, as I don't see much there?

Inspecting the signal spectrum shows frequency peaks at 12.4kHz, 6.2kHz, and 200Hz.

Here's a closeup of the peak at the 12.4kHz carrier. Is this perhaps BPSK instead of FSK?

here is the waveform and my attempts to decode it so far:
waveform sample
| matlab file

Any help much appreciated!

After playing around in GNURadio, this waveform looks like 2 level ASK bursts (OOK bursts plus some residual carrier).

The carrier appears to be at 12.4 kHz.

At the 44.1 ksps sample rate, there appears to be about 7.6 samples/symbol or a symbol period of about 173.5 microseconds (assuming the symbols are not Manchester encoded).

See the top plot in the following image showing two filtered, rectified, and pulse filtered bursts.

Update to respond to comment

What I mean by 2-ASK with residual carrier is that, during a transmission, the original signal is probably something like:

$$(0.1 + 0.5q(t))\cos(2\pi 12400t)$$

where $q(t)$ is the train of unipolar message bit pulses.
With a square pulse matched filter in the receiver, one is going to get ISI, which you can see in the short (suspected) preamble of the bursts.

By the way, I have not checked the phase to see if the bursts are the real part of QPSK with a non-standard rotation (which would result in 4-QAM with constellation points only on the real axis), but that seems unlikely for a low parts count biometric headphone device.

The signal you provided was already low pass filtered with a corner at about 15-16kHz, indicating something was going to decimate this to a 32 ksps signal.

What I did in GNURadio is in the following flowgraph image:

I used the "Hilbert" block with a 201 tap Hilbert FIR filter to make the real signal analytic.  201 taps is obviously overkill, looking at the negative frequency side of the spectrum.

I then used a Blackman-Harris windowed high pass filter, with real taps, with a corner at 4 KHz, a transition width of 1 kHz, and an objective stop band attenuation of 20 dB.  This was to get rid of the mess near DC.  Some of the mess at DC was tracking with the signal, but some of it seemed unrelated.

I then took the envelope with the "Complex to Mag" block.

I then used a rectangular pulse matched filter with taps of $[1, 1, 1, 1, 1, 1, 1, 1]/8$ .  This introduces receiver side ISI, but reduces the noise.

I don't know how you captured this signal, but it would be better if you can do something on the analog side to improve the signal level or SNR, and also preserve the component that was likely at 18.6 kHz, since you're using a sample rate of 44.1 ksps.


Harris corner detection shape of $E(u,v)$ I am taking a computer vision class and I have just learnt about the Harris corner detection concept. A corner is detected when a small shift in a window function defined around the corner results in a large $E(u,v)$ term, which is the sum of squared difference of pixels between the previous window and the next window. There is a post on this forum that has a very nice conceptual explanation about harris corner detectors.
However, given that $E(u,v)$ should be large for small changes in $(u,v)$, why is $E(u,v)$ function usually depicted as a function that has a negative maxima ? Shouldn't $E(u,v)$ value be high when $(u,v)$ is close to $(0,0)$

$E(u,v)$ is zero for zero shift ($u=v=0$), where the two patches being compared are equal. As you increase the distance, the difference tends to increase. Depending on the local shape of the image, the $E(u,v)$ surface will have a different shape: at a corner or a small point, it will steeply increase from 0 in all directions. Along an edge, it will remain close to 0 along one direction, and increase steeply in the perpendicular direction. At flat image regions, the function will remain closer to zero in all directions.
The Harris corner detector (and other corner detectors too) look at the shape of this patch, if it increases steeply in all directions, the detector will produce a stronger response.


Plot of a DTFT is not periodic From this impulse response:

$$h[n]=a^n u[n]\quad\text{with}\quad\lvert a\rvert<1$$

Applying DTFT analysis equation I obtain this frequency response:

$$H(e^{j\omega})=\dfrac{1}{1-ae^{-j\omega}}$$

But when I plot $\mathrm{abs}[H(e^{j\omega})]$ with MATLAB I obtain something that is a-periodic:

What I am doing wrong? This is my code:

W = [-4*pi:0.1:4*pi];
a = 0.9;
H = 1.0./(1.0-e*a.^(-j*W));
plot(W,abs(H));

You switched the $a$ and $e$ in your third line of code. If you write it like this, then it will be behave as expected:

H = 1.0./(1.0-a*exp(-j*W));


Applying a mathematical function to a 16-bit Python array I am using scipy.io.wavfile to read a WAV file and playing it with pyaudio. Before the numpy array is written to a string and output through pyaudio, I would like to be able to apply some sort of function (sin, cos, etc.) to the 16-bit numpy array.

Does anyone have any experience doing something like this? Should I convert to float, apply the function to the float data, then convert back? I am also aware of fixed point arithmetic, but not sure if I can apply it here.

If you're doing ANY audio signal processing in python, I heartily recommend librosa. Your scipy route will work (read as float, process, playback), but it is less direct. Also, I recommend sticking with float, as fixed-point work is more relevant if your code is running on an embedded/mobile system sans-FPU.

You can pip install librosa and get applying math functions quickly. Check out the introduction ipython notebook. Additionally, if you're doing more serious work, there is great integration with scikit-learn machine learning pipelines.

To get you started your code might look like this:

import librosa
import numpy as np

# Load the example track
y, sr = librosa.load("myfile.wav")

# apply mathematical functions or processes to 'y' here:

# If you do this in an ipython notebook (highly recommended as a interactive python env),
# you can play back in one line!
IPython.display.Audio(data=y, rate=sr)


Calculating similarity between two signal plots I have two plots, each having frequency as x-axis and Gain as y-axis. By taking one data set as a reference I have to calculate similarity between them.

The graphs have same values on the x-axis and have same range on the x axis

Can 2D Correlation or Co-variance can do a decent job ? or
should I opt for Fréchet distance or DWT as I have read in some other posts ?

The first plot is the reference plot.

Here are the plots:

Please help!

Why not just use something like the relative "error" between the two?

For example, if your frequency magnitude responses are $G_1$ and $G_2$, then calculate:

$$
ERR = \sum \left | G_1(n) - G_2(n) \right|^2
$$

and then normalize with respect to the reference, $G_1$:

$$
NORMALIZED = \sum \left | G_1(n) - G_2(n) \right|^2 /  \sum \left | G_1(n) \right|^2
$$

Correlation is also a way to go, but it might show some cases where the same shape happens, but at very different frequencies... which might not be what you want.


Deconvolution to Remove Gaussian Blur in 1D Signal (Wiener Filtering?) I've got a set of biology data that I'm trying to denoise (effectively, a population statistic can only be measured convolved with a gaussian of known width)
My problem is this: I can measure (f*g), and I can measure g(x), but I need to know f(x). Because this is real data with slight blips, I can't just use a deconvolve() function.
I've been told by an EE friend that I should be using 1-D Wiener Deconvolution. The problem is, I have no earthly idea how to implement this. Supposedly, there's a Matlab package for it, but I don't have access to Matlab; and all the Python examples I've been able to find (e.g. scikit-learn) are 2D only.
Is there something I'm missing? Is there an easy formula I can implement?
Please, don't just point me to the Wikipedia page/a textbook for 1D Wiener filtering. I've tried to read several of these, but they look like incomprehensible gibberish (I'm a biologist).
EDIT: Also what's up with point spread functions? I assume my PSF should be the known Gaussian g(x), and yet every example I can find of Gaussian deconvolution with Wiener filtering (all in 2d) use a uniform distribution as their PSF instead of a Gaussian distribution.
Thanks!

Approaches
There are many methods for Deconvolution (Namely the degradation operator is linear and Time / Space Invariant) out there.
All of them try to deal with the fact the problem is Ill Poised in many cases.
Better methods are those which add some regularization to the model of the data to be restored.
It can be statistical models (Priors) or any knowledge.
For images, a good model is piece wise smooth or sparsity of the gradients.
But for the sake of the answer a simple parametric approach will be takes - -Minimizing the Least Squares Error between the restored data in the model to the measurements.
Model
The least squares model is simple.
The objective function as a function of the data is given by:
$$ f \left( x \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$
The optimization problem is given by:
$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$
Where $ x $ is the data to be restored, $ h $ is the Blurring Kernel (Gaussian in this case) and $ y $ is the set of given measurements.
The model assumes the measurements are given only for the valid part of the convolution. Namely if $ x \in \mathbb{R}^{n} $ and $ h \in \mathbb{R}^{k} $ then $ y \in \mathbb{R}^{m} $ where $ m = n - k + 1 $.
This is a linear operation in finite space hence can be written using a Matrix Form:
$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| H x - y \right\|_{2}^{2} $$
Where $ H \in \mathbb{R}^{m \times n} $ is the convolution matrix.
Solution
The Least Squares solution is given by:
$$ \hat{x} = \left( {H}^{T} H \right)^{-1} {H}^{T} y $$
As can be seen it requires a matrix inversion.
The ability to solve this adequately depends on the condition number of the operator $ {H}^{T} H $ which obeys $ \operatorname{cond} \left( H \right) = \sqrt{\operatorname{cond} \left( {H}^{T} H \right)} $.
Condition Number Analysis
What's behind this condition number?
One could answer it using Linear Algebra.
But a more intuitive, in my opinion, approach would be thinking of it in the Frequency Domain.
Basically the degradation operator attenuates energy of, generally, high frequency.
Now, since in frequency this is basically an element wise multiplication, one would say the easy way to invert it is element wise division by the inverse filter.
Well, it is what's done above.
The problem arises with cases the filter attenuates the energy practically into zero.  Then we have real problems...
This is basically what's the Condition Number tells us, how hard some frequencies were attenuated relative to others.

Above one could see the Condition Number (Using [dB] units) as a function of the Gaussian Filter STD parameter.
As expected, the higher the STD the worse the condition number as higher STD means stronger LPF (Values going down at the end are numerical issues).
Numerical Solution
Ensemble of Gaussian Blur Kernel was created.

The parameters are $ n = 300 $, $ k = 31 $ and $ m = 270 $.
The data is random and no noise were added.
In MATLAB the Linear System was solved using pinv() which uses SVD based Pseudo Inverse and the \ operator.

As one can see, using the SVD the solution is much less sensitive as expected.
Why is there an error?
Looking at a solution (For the highest STD):

As one could see the signal is restored very well except for the start and the end.
This is due to the use of Valid Convolution which tells us little on those samples.
Noise
If we added noise, things would look differently!
The reason results were good before is due to the fact MATLAB could handle the DR of the data and solve the equations even though they had large condition number.
But large condition number means the inverse filter amplify strongly (To reverse the strong attenuation) some frequencies.
When those contain noise it means the noise will be amplified and the restoration will be bad.

As one could see above, now the reconstruction won't work.
Summary
If one knows the Degradation Operator exactly and the SNR is very good, simple deconvolution methods will work.
The main issue of deconvolution is how hard the Degradation Operator attenuates frequencies.
The more it attenuates the more SNR is needed in order to restore (This is basically the idea behind Wiener Filter).
Frequencies which were set to zero can not be restored!
In practice, in order to have stable results on should add some priors.
The code is available at my StackExchange Signal Processing Q2969 GitHub Repository.
Resources

Deconvolution of Synthetic 1D Signals - How To?
Estimate the Filter Coefficients of 1D Filtration (Convolution).
1D Deconvolution with Gaussian Kernel (MATLAB).
How Is the Formula for the Wiener Deconvolution Derived?


Digitalization of a signal Suppose I have an analog signal I want to transmit. After sampling and quantization, I can map the samples to bits. So now I have a bit stream that represents the original analog signal.
My questions are

Why do we choose to map bits to symbols before transmitting? And
what exactly are these symbols?

One can understand a bit by saying '0' is 0 Volts and '1' is for example 5 Volts. So, 3) How exactly do the symbols modify the analog signal that is transmitted?

What you did by saying a $0$ bit is mapped to a voltage of zero Volts, and a $1$ bit is mapped to $5$ Volts is exactly what it means to map bits to symbols. Your symbols are 'zero Volts' and '$5$ Volts'. A bit is just a unit of information, and you have to decide which type of analog signal you choose to transmit the source information. You always have to map information bits to some analog signal in order to transmit it over an analog channel.
There are many ways to choose appropriate symbols. Your suggested mapping from bits to symbols requires a baseband channel that can also transmit DC, because your transmitted signal will have a DC value of about $2.5$ Volts, assuming that $0$s and $1$s are equally likely. There are baseband channels that are AC-coupled, i.e., you can use them to transmit very low frequencies, but not DC. In that case you could choose a positive voltage for transmitting a $1$, and a negative voltage of equal magnitude for transmitting a $0$ ("binary antipodal signaling"). There are many standard ways to encode bits for transmitting them over a baseband channel. Here you can start reading up on these so-called line codes.
Note that you can also map more than one bit to a symbol. E.g., you could choose voltages $-3$, $-1$, $1$, and $3$ to encode all possible two-bit sequences. This makes your symbol rate (baud rate) half the bit rate, which requires less channel bandwidth. There's a trade-off between bandwidth and noise immunity. For a fixed transmit power, the symbols are now closer to each other than before, so it is more likely to confuse one symbol for another when there's additive noise on the channel.
There are also channels that have a bandpass characteristic (such as a radio channel). In that case you need to use carrier modulation, which shifts the spectrum of the transmitted signal to the desired frequency band. You can use two orthogonal carriers at the same frequency, such as $\sin(2\pi f_c t)$ and $\cos(2\pi f_ct)$ to transmit two symbols at the same time (quadrature modulation). In that case, symbols are represented by complex numbers (because you have two symbols at a time) and you get a two-dimensional symbol constellation.
It is common to use constellations with many more than just two symbols, e.g., $256$-QAM, which uses one symbol to encode $8$ bits.
All the above gives a very basic and simplistic view of how encoding and modulation works. If you want to dig a bit deeper you should get a good textbook on digital communications.


Estimating Convolution Kernel from Input and Output Images Given an original image and a convolved version of it, I need to calculate the convolution kernel. For example, given:

and

I would like to calculate the convolution kernel that generated the second image from the first (a gaussian in this case).

The images will be similar to the example, but I can change the configuration slightly to make the process easier (change the shape, have a number of shapes etc.)

Take Fourier Transform of both original image and blurred image, then divide Fourier transform of blurred image by Fourier transform of original image. This will give you the Fourier transform of kernel (except at frequencies with zero value which lead to division by zero). Then take inverse Fourier transform and find the kernel in spatial domain.


Sampling pure tone sine waves What would happen if I am using the maximum frequency as the sample rate for sampling a pure tone sine wave?
For example, a $10\ \rm kHz$ sampling frequency for a $10\ \rm kHz$ monotone sine wave. What effect would aliasing produce at the output?

In contrast to common misunderstanding, aliasing does not affect the process of sampling, rather it affects the process of reconstruction. i.e., the samples themselves do not contain error, but their interpretation at the reconstruction (or anything related with it) will be in error.
So, if you sample a sine wave $$x(t) = A \sin( 2 \pi f_0 t + \theta) $$ at a sampling frequency of $F_s = f_0$, (which is undersampled according to Nyquist-Shannon sampling theorem for bandlimited signals), then your samples will be $$ x[n] = A \sin(\theta) $$ a constant signal! When reconstructed back into continuous waveform at the interpolator output, you will get a DC signal, $x(t) = A \sin(\theta)$, as opposed to the original sine wave.
The error is due to the false assumption at the interpolation stage that those samples $x[n]$ came from an original continuous-time signal which was band-limited to $f_0/2$, Nyquist frequency of sampling frequency $f_0$, which is not the case in reality.


am I recording LTE uplink signals? I have this app in my phone which tells me the LTE UL frequency and power:

so I went with my SDR and measured the received signal at that frequency:

Now, how can I be sure this is the transmitted signal from my phone? As I made a call to and from this phone and nothing changed.

I did a second experiment when the phone switched to LTE band 20 and the spectrum looks quite different:

I also recorded this video when the phone is just in "standby" and this other video when the phone is sending an email with an attachment.

I have the feeling I'm not seeing the correct information in the LTE band 3 because my dongle (this one) is not actually receiving anything there, although it shows an spectrum? My goal here is to be able to detect my phone's transmitted signal and then, to make sure I can distinguish it before I invest in a more expensive SDR to be able to monitor higher bands.

There's a lot of different domains of knowledge coming together here, so I'll split my answer into multiple sections, each answering an implicit question that you raise in your explicit question. Hope that helps!

Can your RTL-Dongle actually receive at 1.72 GHz?

So, first the bitter pill: There's a lot of sellers out there that offer RTL dongles and claim that they've got an R820T tuner, but actually ship an Fitipower tuner that covers far less spectrum. I find that the most likely explanation here, especially considering the product description and the "tech specs" table disagree. Said table only claims:

  Receiving Frequency:48.25-863.25 MHz Tuning Range

So, possibly you've got an RTL-Dongle with a fitipower FC0012, not a R820T (which went out of production a while ago). RTL-SDR comes with the rtl_test utility, but possibly your software can show you the tuner installed in that device, too.

However, considering your display shows the right frequencies, let's assume it actually tuned to ~1.720 GHz when you instructed the device to tune there.

Are you using an antenna that's suitable for these frequencies?

Just mentioning that, because I've supported far too many people that told me after a couple of hours that they were trying to receive FM broadcast radio at 100 MHz using a Wifi antenna, and vice versa.

Are the powers in these different applications related?

Then, first thing to notice: The powers displayed in your SDR software are not related to the power your App shows for reasons:

The power shown is the downlink power. You're observing uplink. Uplink should probably have significantly more power than downlink at your handset.
I don't know where the app gets that data. But, I'll go out on a limb and say it's not really an RX power estimate, at all. LTE specifies a whole set of quality metrics, all of which are useful to assess the error probability of a receiver. None of these is a pure power measurement (because that doesn't matter to a device, actually. What matters are SNR measures that incorporate noise, but also things like adjacent channel interference, clock errors, …).
An RTL dongle doesn't even remotely give calibrated values. What it sends to the computer are simply 8 bit numbers, relative to the ADC's full scale, whatever that means at any given gain, frequency, rate, filter setting, temperature, component tolerances, moon phase, amount of black cats in neighborhood… . I don't know why your program shows you "dBm", but unless you have run some calibration on your dongle for exactly your gain, frequency, and sampling rate settings, these numbers mean nothing. All you can say is that 3dB more means twice as much power.

Can I analyze LTE signals using an RTL Dongle, at all?

LTE uses a total bandwidth much larger than 3 MHz, so, if at all, you're looking at a small excerpt of your uplink; possibly even one where your handset doesn't even get an allocation of resources.
LTE is very flexible about assigning resources to different user equipment in the same cell, and very quick – unless you have a very aggressive means of holding spectral peaks, it might be pretty impossible to spot your own contribution in a cell.

Especially the downlink uses OFDM, and usually also power control, so if your handset is the main downlink user in a cell, all that you should be seeing if you had an SDR that could actually see the whole LTE downlink channel, would be what would look like a "flat-top" section of spectrum that is basically elevated white noise, possibly exhibiting non-flatness due to frequency-selectivity of the channel.

If you look at a section of OFDM, all you'll see is a noise floor, unless you manage, by luck or synchronization, use an FFT that exactly "hits" (integer factors) of the subcarrier spacings (in which case you'd be lucky to see staggered sinc spectra of the OFDM subcarriers). You'd need a very long FFT and very good SNR (which you don't seem to have) to see the sinc carriers unsynchronized.

Are you actually observing your phone's LTE signals in your pictures?

Doesn't look like it.

Your 800 MHz spectrum shows something 200 kHz wide with a rectangular shape; that might be a bunch of LTE subcarriers (those have 15 kHz spacing), but it's definitely not coming out of the noise floor far enough to be your phone's resource blocks (also, questionable whether you'd get a resource block allocation that's static for that long).

Considering the periodic spurs and the intereferer at ~1719.35 MHz that looks a lot like a mixed-up rectangular wave, I'll go with no, you're likely not looking at LTE here; there should be far more dynamics in the waterfall.

Is your ultimate goal achievable?

  My goal here is to be able to detect my phone's transmitted signal and then, to make sure I can distinguish it before I invest in a more expensive SDR to be able to monitor higher bands.

Detecting activity of any one phone in vicinity

As said, your TX signal should be much stronger close to your handset than it's RX (which really is just the simple physics of spatial power density going down with the square of distance). Especially when your phone displays "few bars" (i.e. it has low receive SNR), it will have to use a lot of power when transmitting (which is why I'm pretty sure you're not observing your phone in either screenshot).

A simple power detector with a high-enough threshold should do an approximately good enough job. You'd still not be sure whether it's your phone, or the phone in the next room. LTE relies on power control, i.e. your phone will only use as much power as necessary to achieve the rate it needs (that's not something that your phone is allowed to choose – it's necessary to keep the dynamic range of the receiver on the base station side manageable and to keep interference at bay).

Detecting activity of a specific phone

Now, this get's a lot more complicated, as you'll have to analyze the LTE signal. Not quite sure from the top of my head, but I'd go with "no, from the uplink signal alone, you can't say which phone is currently transmitting".

LTE is a cellular infrastructure network, not an ad-hoc burst network; hence, stations are assigned times when they should be transmitting by the base station using the downlink. Inherently, that means that the station doesn't need to indicate who it is when using the uplink; the base station knows because it assigned the uplink resources, who's using the uplink. Thus, you'd have to observe the downlink more than the uplink to know when you're phone's talking. Luckily, LTE is generally pretty well-encrypted, so you shouldn't be trivially able to extract very specific information from the downlink about your phone, but you might at least be getting unencrypted local IDs of the different stations in the cell from the downlink, if your SDR receiver software is able to demodulate the signal.


Cross-correlation with SDPSK modulation I'm trying to detect a sync word in a signal using Symmetrical Differential PSK (bit 0 = +90° phase, bit 1 = -90° phase).

Given the signal is encoded in the phase difference between two symbols, the absolute phase is irrelevant. How would I go about using cross-correlation to detect the sync word?

I have four possible orientations for the sequence considering the absolute phase, but for performance reasons, I would prefer not to run the cross-correlation more than once. Should I calculate the phase difference of each symbol with the previous (i.e. extract the bits) and correlate on that or is there a way to do it without first converting the signal to bits?

Convert your input signal to an analytic signal with a Hilbert transformer, if your input signal is real.  Then run a single complex cross correlation filter.  Assuming no frequency offset, the phase of the correlator output at the correlation peak will be the input signal's phase shift from your reference signal.


$\mathcal Z$-transform if the output is given A impulse response for a LTI system is given by:

  $$h[n]=\left(\frac{2}{3}\right)^n u[n]+2 \left(\frac{1}{5}\right)^n u [n]$$

  and if the putput for the system is given by:

  $$y[n]= \left(\frac{1}{3}\right)^n u [n] $$

  What is the input in the $z$-domain ($X(z)$)?
  And what is the output in the time-domain ?

I'm not sure what is the correct procedure for solving this type of problem. Can you just simply transform $y[n]$ into $Y(z)$ and then solve for $X(z)$?

$$\frac{Y(z)}{X(z)}=\frac{3-\frac{23}{12}z^{-1}  }{\left(1-\frac 23z^{-1}\right)\left(1-\frac 15z^{-1}\right) }$$

Then I replace $ Y(z)$ with $\frac{1}{1-1/3z^{-1}}$ and then I isolate $ X(z)$.

$$\frac{\frac{1}{1-\frac 13z^{-1}}}{X(z)}=\frac{3-\frac{23}{12}z^{-1}  }{\left(1-\frac 23z^{-1}\right)\left(1-\frac 15z^{-1}\right)}.$$

Then I get $X(z)=\frac{6z^-2-39z^{-1}+45}{23z^{-1}-11z^{-1}+135}$ and then I transform it into a sequence in the time domain.

Is it right? or is there a better method for solving this.


What function to apply a Metamorphopsia (grid warp) translation? How would I go about finding or creating a function to warp images dependent on grid lines surrounding an origin (I don't wish for a grid on my actual images):

It should work in practically the same manner as the envelope distort tool in Adobe Illustrator, though I wish for mine to be adjusted according to amplification of the warp, so small multiplier would lead to less warp, whilst higher would lead to greater warp encompassing more of the image surrounding the origin. Also of note is I wish for my images to keep their colours. I use Python. The distortion surrounding the point can be random, just not extreme.

Background info: Compounding adversarial attacks into biological relevance and then researching the biological relations of different NNs. Created images previously (through illustrator) though the amplification will allow me to go back and forth between early and late stages of diseases, and obviously just make it far simpler for me to increase my dataset tenfold, which I desperately need also!

Sorry for using medical term for the translation name, image processing is not my area of expertise in the slightest!

What you are looking for is called a "Geometric Transformation".

When applied to image processing, it comes down to a re-mapping of pixel positions. In optical systems, they can be used to represent certain distortions caused by lenses. In fact, ray transfer matrices end up being applied as geometric transformations to images, since they represent exactly how is a pixel displaced from its original position as its light ray travels within an optical system.

There are two things that need a little bit more consideration with geometric transformations:

How to express the mapping of pixels; and
How to apply the interpolation.

Expressing the mapping of the pixels might be relatively "easy" for known distortions (as mentioned earlier) but what is depicted in the body of the question here seems to be arbitrary (?). In any case, you are going to have to somehow express the function that takes a pixel from location $x,y$ to location $x',y'$ in a way that makes sense to the problem at hand. The way that image processing software (such as Illustrator that you mention) does it, is by setting up a grid that divides the image into a set of patches. These lines can be piecewise linear or bezier curves. This "control" (in terms of a user interface element) is not exactly straightforward to build.

Interpolation is absolutely necessary here because the geometric transformation might "ask" for pixel values between known pixels. There are obvious ways out of this problem (e.g. linear interpolation) but it is very likely that you will need a different method that makes sense within your application.

In Python, geometric transformations can be handled by Scipy's geometric_transformation(). Its use is straightforward provided that you know the input and output pixel locations.

Hope this helps.

EDIT:

  How would I go about remapping these pixels considering I don't want singular pixels remapped but the entire areas 'pulled' in a direction.

This is what the geometric transformation does, as long as the input/output point pairs are provided.

  Will interpolation work with colour in my image?

Yes.

  Do you have any suggestions where I could learn how to do this practically? Arbitrary can come after I've sorted this first part out!

Yes, have a look at the links i provided. From the mathematical point of view, we are talking about a series of linear transforms that map a "rectangle" to a distorted "rectangle". See the documentation of the geometric_transformation().


IIR/FIR filtering before decimation? I am trying to design a low pass filter to apply before resampling the sensor's output at a lower rate, however I am not quite sure how to choose between the IIR and FIR fitlers. Matlab by default, filters the data with an 8th order Chebyshev Type I lowpass filter with
cutoff frequency .8*(Fs/2)/R, before resampling. but how can I choose the best low pass filter type ?

It depends on your specification really.

When you're decimating, some degree of aliasing will occur. Your tolerance of this aliasing, signal delay, and changes in the phase response of your input signal will determine the best approach to take.

A very long FIR might provide a very steep cutoff and linear phase performance (depending on design), but will delay the signal (group delay), may cause pre-ringing and can be computationally demanding. FIRs however also allow for elegant 'polyphase' implementations (which reduce computational demand by effectively breaking the filter down into smaller component filters).

IIRs on the other hand have the potential to have a shorter delay, at the expense of a non linear phase response and perhaps a more shallow roll-off. You could achieve a steep rolloff with an IIR but it would likely be quite high order and hence susceptible to coefficient quantisation and/or stability problems.


Performance metrics for analysing the performance of various filters We are working with speckle noise reduction in ultrasound images by using various despeckling filtering techniques. We compare the performance of each filtering technique by using Quality metrics such as mean square error,signal to noise ratio,peak signal to noise ratio,Edge preservation index,Speckle index.

We know that mean square error metric indicates how different the images being compared are. Therefore the lower its value is,the closer the estimated image is to the original image.
The peak signal to noise ratio must be high which indicates an improvement in speckle reduction.But when we calculate these metrics, mean square error is high and Peak signal to noise ratio is small. any one please say me why this happens so.

The original image that  we have used is shown below

Then we corrupted the original image with 3% noise

The formula that we used to calculate these metrics is given below

The result which we obtain while using median filter is given below

They are getting values like that because they may be reading the image as double while calculating the metrics, you can use im2uint8 to convert the double into uint8 and estimate the metrics. you will get the desired output. Though it does not mean that their values are incorrect.


Explicit a succession using z inverse transform Is it possible to explicit $y(n)$ of this mathematical succession in recursive form using z inverse transform?:

$ y(0) = 1 \\ y(n+1) = 2y(n) + 3 $

I can't write this as a linear system in z (linear difference equation) because there is the $3$ that doesn't depends on input or output.

I know some methods to solve this but I'm trying to use z inverse transform with no success.

If you are allowed to define your system to be causal so that it is non-zero only for $n \ge 0$, then you can use the unit step function $u[n]$ so that your system becomes:

$$y[0] = 1$$
$$y[n+1] = 2y[n] + 3u[n]$$

Which does have a Z-transform representation.


Equalize Speaker using impulse response So I have this school project where my aim is to improve a loudspeaker in general, my idea is to mesure the impulse response of the room using a mic, then use this mesure to conceive a loudspeaker that can compensate the effects of the room.
My question is how could I exploit the room impulse response in conceiving that filter let's, and what kind of filters should I use? would that even be possible?

PS: sorry if the question sounds dumb, I have zero experience and knowledge with this kind of stuff.

Thanks in advance!

This depends a bit on what you want to get out of this and how much effort/work  you are willing to put in.

Doing a room EQ that actually works and makes it sound consistently better is quite complicated. There are commercial systems available but they tend to be complicated and expensive or tied to a specific product.

A really good freeware option to play around with is Room EQ Wizard https://www.roomeqwizard.com/

If you want to do this yourself, I suggest you break it down in a few separate steps that can be tackled individually

Figure out how to make a good measurement in a room (noise, reverb, microphone, ADC, SNR estimation, etc.)
Figure out what to measure (using the tools from step 1) that's representative enough of the entire room/speaker/listener situation
Extract perceptually relevant parameters from these measurements
Turn these parameters into some sort of correction that would move the parameters towards some perceptual target
Design a filter and/or algorithm based on this correction
Design the process/HW/system to actually implement and apply the filter or correction

All of these steps have been studied in detail and there is a fair bit of literature available on any of them.


How condition for existence of Fourier transform is valid? The condition for Discrete time Fourier transform to exist for function $f(n)$ is given as

$$\sum_{-\infty}^\infty |f(n)|  < \infty.$$

In case of continuous Fourier transform the difference is summation is replaced by integration in the above equation

It looks something
strange for me because on left side we are
taking the signal within 1. infinity
limits. Also,2. it is again of summation type and we are expecting right side to be less than infinity. How is it possible?

Let us take two  examples

  Suppose there is unit step signal $u(n)$
  with constant unity magnitude
  ranging from o to $\infty$ .then we
  can clearly see by putting value
  in the above formula

L.H.S=1+1+1+......upto ${\infty}$
So R.H.S. =$\infty$

So can I say that DTFT of unit step signal doesn't exist?

  Similarly, if I represent the  power signal like $cos(wn)$  using the exponential terms and compute it's value from the formula of condition of existence ,it will come out to be
  $\infty$ .

So can I say that DTFT of power signal like
 $cos(wn)$ doesn't exist?

As mentioned in Batman's answer, the condition of the sequence being absolutely summable is only sufficient but not necessary. The Fourier transform can be extended to $\ell_2$ sequences, i.e. sequences for which

$$\sum_{n=-\infty}^{\infty}|f[n]|^2<\infty$$

is satisfied. A further generalization is possible if you allow distributions and their derivatives (such as the Dirac delta function) as Fourier transforms. Then also sequences with polynomial growth have a Fourier transform in that generalized sense. Just sequences with exponential growth have no Fourier transform. However, they can be handled by the $\mathcal{Z}$-transform. Note that everything said above also holds for continuous-time functions (if your replace "$\mathcal{Z}$-transform" by "Laplace transform").

Sequences with a constant amplitude (such as the step function or a sinusoid) will have a Dirac delta in their Fourier transform. Sequences with a linear growth will have the (generalized) derivative of a Dirac delta in their Fourier transform, etc.

And as for your question how it is possible for a sum over infinitely many elements to be finite, your should as a first basic example read up on the geometric series.

EDIT:

So if we allow a Dirac delta impulse as a Fourier transform, we can ask what the corresponding time function would be. By taking the inverse Fourier transform of a shifted Dirac delta impulse $\delta(\omega-\omega_0)$ we get

$$\frac{1}{2\pi}\int_{-\infty}^{\infty}\delta(\omega-\omega_0)e^{j\omega t}d\omega=\frac{1}{2\pi}e^{j\omega_0t}$$

which follows from the sifting property of the Dirac delta impulse. So we now have the Fourier transform pair

$$e^{j\omega_0t}\Longleftrightarrow2\pi\delta(\omega-\omega_0)$$

from which the transforms of $\cos(\omega_0t)$ and $\sin(\omega_0t)$ follow directly by rewriting these functions as $\cos(\omega_0t)=\frac12(e^{j\omega_0t}+e^{-j\omega_0t})$ and $\sin(\omega_0t)=\frac{1}{2j}(e^{j\omega_0t}-e^{-j\omega_0t})$.


Transfer Function vs Frequency Response Talking about Butterworth filters my prof said that a filter is described by its frequency response $H(\omega )$.

What I think is that filter is really described by its transfer function, since there can exist many transfer functions that share the same $H(\omega )$.
Was my prof wrong?

Was my prof wrong?

No. You're right, you can find multiple systems with the same frequency response, or the same transfer function.

But if you're using something as a (linear, which is often implied) filter, then all you care about is the frequency response (otherwise, you'd probably not be calling the thing "filter"). So, when I define a filter, it's usually sufficient to specify $H(\omega)$.


Friis law and the power delay profile Assume we have communication system between transmit and receive over wireless channel.

The Friis path loss model is given by the following equation

$$P_r = P_t G_t G_r (\frac{\lambda}{4\pi D})^2$$

Where $P_{r(t)}$ is the received (transmit) power, $G_{t(r)}$ is the gain of the antennas at transmitter (receive) and $\lambda $ is the signal wavelength while $D$ is the distance between transmit receive ends.

In theoretical (academia) and simulation scenarios (for example MATLAB)I have come across examples wireless channel where a power delay profile is defined, that is time of arrival (TOA) of rays versus amplitude of those rays. For example, assume I have five rays arriving and delay spread $40 ns$ with amplitude

$$channel= [ 0.5+i0.5, 0.3+i0.3, 0.2+i0.2, 0.1+i 0.1, 0.3 ] $$

$$TOA= [0, 10*10^{-9} , 20*10^{-9}, 30*10^{-9} , 40*10^{-9}]$$

My simple question is: How are the two concepts related??

Thanks

Basically, multipath propagation causes time dispersion (hence causing a delay spread ($\sigma_{\tau}$)) and that dispersion influences the fading characteristics of the wireless channel.So delay spread is related to the fading characteristics of a channel and not really to the path loss. Particularly, if $\sigma_{\tau}$ is greater than the symbol period ($T_s$) then the wireless channel is said to be frequency selective, meaning that a signal will interfere with itself causing additional fading.

You should also note that the rms delay spread is not the maximum delay that your signal experiences as I think you may asssume, $\sigma_{\tau}$ can be calculated as follows:

$\sigma_{\tau} = \sqrt{\overline{\tau ^2} - \overline{\tau} ^2}$

where

$$\overline{\tau} =  \frac{\Sigma P(\tau _i) \tau _i}{\Sigma P(\tau _i)},$$
$$\overline{\tau ^2} =  \frac{\Sigma P(\tau _i) \tau _i ^2}{\Sigma P(\tau _i)}$$
and $P(\tau _i)$ is the power of the component with delay $\tau _i$


Can the order of filtering and downsampling be exchanged? Say for example, I have a signal with maximum frequency at 600Hz, and I originally oversampled it at 2000Hz. Afterwards I was only interested in frequency less than 500Hz. I guess I could first filter out frequency >500Hz, then downsample by taking every second sample point; or I could first downsample by taking every second sample, then filter out frequency >500Hz. Is there any problem with the two methods? Do they generate different results?

In Short: The order does matter.

Detailed:
Downsampling by a factor of 2 without filtering, will cause aliasing in frequencies of 500Hz and more.
Since your signal is bandlimited to 600Hz, the frequency range of 400Hz-500Hz will be corrupted.
So (theoretically) the frequency range 0Hz-400Hz will be the same with or without filtering.
If you were only interested in this range, it would have been possible to avoid the filtering.


What information is lost, if any, in normalizing Seismic/Earthquake time-frames in an STFT to equal variance? This might be more of a statistical question, but let us assume that we take a spectrogram, or STFT, where time is on the x-axis, and frequency is on the y-axis in this matrix. What information might be lost - or 'changed', if we decide to normalize each column by its own standard deviation, so that the variance of each column is equal to 1.

I understand that we would not be able to see sudden changes in overall power changes over time, etc, that is obvious, but what other features might be contorted and/or lost if we did this?

Thanks.

Although the specifics of your data will be the ultimate factor in how this sort of normalization will impact what you're trying to measure, I think it's useful to look at extreme cases to get a sense of what might go wrong.

For example, suppose the time samples of your STFT just happened to come out one-hot, with 1s in either the highest-frequency or lowest-frequency bins and 0s everywhere else :

STFT = [[1, 0, 0, 0, 0],
        [0, 0, 0, 0, 1],
        [1, 0, 0, 0, 0],
        ...
        ]

If you subtracted out the mean values from each time sample, then clearly you'd get some "bleeding" of the power information across frequencies. In particular, the middle-frequency bin would get some nonzero value, even though in your unnormalized STFT values, that bin never has any power.

What about just dividing each time sample by its standard deviation, like you're asking about ? In this case you'll just be changing the value of the nonzero frequency bin. This would certainly impact the mean value across time of each frequency bin, since all nonzero values in your spectrogram would now be 1+$\epsilon$ rather than just 1.

This is just one possible thought experiment that you could do to explore these effects. I'd say in general that the statistics of the frequency bins across time would definitely be affected by normalizing each time frame independently, but again that the specifics will depend heavily on the data that you're analyzing.


OFDM pilot allocation block type and comb type As known, the pilot insertion in OFDM can be either bock type or comb type, I think comb type is appropriate for double-selective channel while block type is for frequency selective channel:

In the above figure, it means that we perform the iFFt in Vertical direction and then send the data horizontally. Suppose we have the size frame is $M$ x $N$ which means we perform $M-$point iFFT for each column before transmitting the whole frame through the channel.

At receiver side, we receive the frame  as follows, in figure $(b)$ pilot location index is 1:4:end, and in figure $(a)$ pilot location index is $x$:4:end, and $x = 1 : N$. is that right?

The received frame is of length $MN$x$1$, how can we deal with it to estimate and equalize the channel? suppose we use the LS for estimation? ...  I think in figure $(b)$ it's possible to reshape the received frame into $M$x$N$ after that  perform $M$-point FFT and then get the estimated channel $h_{est}$ as $h_{est} = P_r./P_o$  such that $P_r$ and  $P_o$ denote received and original pilot, respectively. finally we interpolate to get the estimated channel. Is that Right?

If so, What's about the channel estimation in figure $(a)$ ? how can we deal with it ?

In both cases, you need to interpolate the channel between the pilots you've got.
Both cases are typically suboptimal, since they'd only work perfectly for (a) actual block-fading (which is a convenient model, but doesn't look like reality) or (b) for a channel that is perfectly interpolatable from just a few points of observation in frequency (but that would imply you have designed an OFDM system with too many subcarrier, and that has other downsides).

You'd usually need both, because over the course of a longer data transmission, the channel might change, so there's need for regularly estimating more than just a few subcarriers. And because you cannot just guess the subcarrier state from a pilot some (up to) $\frac Q2$ subcarriers away, in general.

Modern systems tend to solve this by having something like (a) at the beginning of a frame, and then distributing pilots in a pseudorandom, or in a diagonal way throughout the frame, so that you get updates on the channel behaviour over time. Due to having updated your channel state info underway, you can space your pilot blocks further apart (i.e. increase $Q$), thus not reducing your data symbol rate by a lot.

LTE does that. There's a lot of resources you can find when searching for "LTE pilots", but linking to any of these here is a bit futile, since I don't know which background you're coming from.


How to automatically identify a surge in a periodic signal Say I have a periodic signal i.e. Sine wave, which sometimes can have spikes that are 100% larger than the peak value. These spikes can be periodic as well.  I'm trying to find an algorithm that could detect these spikes without human intervention.
I was thinking for each point to average the value with the surrounding 3 or more points but it doesn't seem to be working.
Any help is appreciated.

My suggestion would be to plot the frequency response by taking the FFT of a reasonable length to see the spectrum. Apart from the two spikes in FFT output (corresponding to sinewave) you could see the DFT of the periodic spur signal with approximate amplitude of $N$ where $N$ is the period of the combined signal. The DFT of spur is a constant amplitude signal, which means, for all FFT bins, it will be constant value. If your signal is noisy, you could easily average out these values to find the amplitude of spur. You should first subtract amplitude of sinewave from the FFT before computing this. Here is the code I tried

clc
clear all
n=0:9999;
x = sin(pi/5*n);
spike = repmat([1 zeros(1,9)],1,1000);
x1 = x + spike;
plot(abs(fftshift(fft(x1,1000))))

I got the amplitude of spur signal to be $N/N_{fft}$.


Why the average signal power should be 1 When I read some books or papers related into signal modulation, for example, frequency-domain OFDM signal, it's always mentioned that the average power of signal is assumed to be 1.
Why does that assumption is taken? what's if that average signal power is more or less than one?

Assuming an average power of one for a signal, or for the noise, is a common assumption, mainly to simplify equations. It is assumed that for other powers, a scaling by the actual power can be done.
For instance, many denoising techniques set a threshold that is proportional to the noise power (the square root of it). In the theory, one often assumes (for computations) that $\sigma=1$, design an optimal $\lambda$ and set the actual threshold to $\lambda \sigma$.
The actual threshold is multiplied according to the estimated noise amplitude.


Significance of Sinusoidal signals Can anyone tell what are the significance of sinusoidal signals in the field of science and engineering?
sorry for being noob,

Significance of (pure) sinusoidals in engineering (summary):

1- Periodic signals: Pure sine-waves, represent a useful and simple tool in modeling various physical phenomena, such as standing waves, rotations, alternating currents etc.

2- Fourier analysis: Periodic signals of any engineering sort (not necessarily of pure sine waves) are conveniently decomposed into pure sine waves acording to Fourier series analysis.

3- Linear Systems: strictly speaking complex exponentials $e^{j\omega t}$, are eigenfunctions of LTI systems but sine waves are deeply related. More specifically, the frequency of the sinusoid is preseverd but the amplitude and phase are changed by LTI systems, according to the system's frequency response $H(\omega)$. For instance, for a zero phase system (a symmetric noncausal impulse response) where $H(\omega)$ is real, then the pure sinusoids will also display the eigenfunction property.

4- Nonperiodic Signals: By the extension of Fourier series analysis to integral analysis of continuum frequencies, almost any practical signal is analysed by periodic sine waves.


Matlab version rpmordermap and rpmfreqmap Are those functions only available in Matlab 2015b version or they are included in 2015a version also? Thanks in advance.

At the very bottom of MATLAB manual pages there is always an info about the time when given function was first introduced. For both of rpm*map functions:


Does baseband signal has negative frequency? Does baseband signal has negative frequency? In wiki baseband signal means that can include frequencies that are very near zero. I wonder whether baseband signal has negative frequency.
Thank your answers..

It might help to watch my YouTube video explaining negative frequency: https://youtu.be/gz6AKW-R69s
I'm a professor and I've been building what I think of as an on-line video textbook for Signals, Systems, and Digital Communications.


IDTFT of $\sum_{k=-\infty}^{+\infty}(u(\Omega+\pi)+u(\Omega+\frac{\pi}{4})-u(\Omega-\frac{\pi}{4})-u(\Omega-\pi))\star \delta(\Omega-2k\pi)$ Compute the IDTFT of the following signal:

  $$X(\Omega)=\sum_{k=-\infty}^{+\infty}\left(u(\Omega+\pi)+u\left(\Omega+\frac{\pi}{4}\right)-u\left(\Omega-\frac{\pi}{4}\right)-u(\Omega-\pi)\right)\star \delta(\Omega-2k\pi)$$

Using the IDTFT definition, I obtain:

$$x[n]=\frac{1}{2\pi}\int_{-\pi}^\pi X(\Omega) e^{j\Omega n}d\Omega = \frac{1}{2\pi}\left(\int_{-\pi}^{-\frac{\pi}{4}} e^{j\Omega n}d\Omega + 2\int_{-\frac{\pi}{4}}^{\frac{\pi}{4}} e^{j\Omega n}d\Omega + \int_{\frac{\pi}{4}}^{\pi} e^{j\Omega n}d\Omega\right) = \frac{1}{2\pi jn}\left( e^{-j\pi n/4} -e^{-j\pi n} +2e^{j\pi n/4} -2e^{-j\pi n/4} + e^{j\pi n} -e^{j\pi n/4} \right) = \frac{\sin(\pi n/4)}{\pi n}$$

However, from the transformation table we know that this $x[n]$ has a different transform:

$$X_d(\Omega)=\sum_{k=-\infty}^{+\infty}\left(u\left(\Omega+\frac{\pi}{4}\right)-u\left(\Omega-\frac{\pi}{4}\right)\right)\star \delta(\Omega-2k\pi)$$

which is clearly different from $X(\Omega)$. So, I am asking: what have I done wrong here in this IDTFT definition application?

You made a mistake in this last equality

$$x[n] = \frac{1}{2\pi jn}\left( e^{-j\pi n/4} -e^{-j\pi n} +2e^{j\pi n/4} -2e^{-j\pi n/4} + e^{j\pi n} -e^{j\pi n/4} \right) \neq \frac{\sin(\pi n/4)}{\pi n}$$

The right way to do this would be:

$$\begin{align}
x(n)&=\frac{1}{2\pi jn}\left( e^{-j\pi n/4} -e^{-j\pi n} +2e^{j\pi n/4} -2e^{-j\pi n/4} + e^{j\pi n} -e^{j\pi n/4} \right) \\
&=\frac{1}{\pi n}\left( \frac{e^{-j\pi n/4}-e^{j\pi n/4}}{2j}+ \frac{e^{j\pi n} -e^{-j\pi n}}{2j} +\frac{2e^{j\pi n/4} -2e^{-j\pi n/4}}{2j}  \right) \\
&=\frac{1}{\pi n}\left( -\sin\left(\pi n/4\right)+\sin\left(\pi n\right)+2\sin\left(\pi n/4\right) \right)\\
&=\frac{1}{\pi n}\left( \sin\left(\pi n\right)+\sin\left(\pi n/4\right) \right)
\end{align}$$

Which makes sense as one can see that the DTFT is the sum of two rectangular windows, one of width $2\pi$ and one of width $\pi/2$, each corresponding to each $\mathrm{sinc}()$.

EDIT:

I've just noticed judging by your comment that the mistake you made was to think that

$$\frac{\sin(\pi n)}{\pi n} =0$$

As you already know, the numerator is $0$ for all $n$... except for $n=0$. In that case, the denominator also is zero, so you have an indetermination. The same happens in the continuous case. Remember that we assume

$$\mathrm{sinc}(0)=1$$

by taking the limit when $t\to0$. To respect the fact that the discrete $\mathrm{sinc}$ is a sampled version of the continuous one, then it equals $1$ at the origin too. So we can state that:

$$\frac{\sin(\pi n)}{\pi n} =\delta(n)$$

Notice that the DTFT you wrote can also be expressed (using the fact that convolution is distributive) as:

$$X(\Omega)=\sum_{k=-\infty}^{+\infty}\left(u(\Omega+\pi)-u(\Omega-\pi)\right)\star \delta(\Omega-2k\pi) + \sum_{k=-\infty}^{+\infty}u\left(\left(\Omega+\frac{\pi}{4}\right)-u\left(\Omega-\frac{\pi}{4}\right)\right)\star \delta(\Omega-2k\pi)$$

The first term is a window of width $2\pi$ that is $2\pi$-periodic... So it's basically $1 \ \forall \Omega$.

$$X(\Omega)=1 + \sum_{k=-\infty}^{+\infty}u\left(\left(\Omega+\frac{\pi}{4}\right)-u\left(\Omega-\frac{\pi}{4}\right)\right)\star \delta(\Omega-2k\pi)$$

Now it's easier to see that the IDTFT I got at the beginning of the question corresponds indeed to the given DTFT.


Phase response of a FIR Hilbert transformer The graph below shows a linear phase response of FIR Hilbert transformer with a constant diff of -0.06135923. I am not sure if I understand the concept of the phase response correctly, but as the quadrature is the input signal shifted by PI/2, why is the phase response not a straight horizontal line at PI/2?

The code to generate the figure was copied directly from Scipy's documentation (the example at the bottom of the page) except that I did not unwrap the phase angle and I used the following filter coefficients:

N = 10
b = [2*np.sin(np.pi*n/2)**2/(np.pi*n) for n in range(1, N+1)]
b = np.append(-1*np.flip(b, 0), np.append([0], b))*np.hamming(2*N+1)

I tried the same with different values of N and each time the behavior of the phase response line was the same, just the steepness was varying.

The slope represents a delay of $N$ samples. You can center the impulse response to time zero by multiplying the transfer function by $e^{i\omega N}$ (transfer function of a negative delay of $N$ samples) or np.exp(1j*w*N) in Python:

from scipy import signal
import numpy as np
N = 10
b = [2*np.sin(np.pi*n/2)**2/(np.pi*n) for n in range(1, N+1)]
b = np.append(-1*np.flip(b, 0), np.append([0], b))*np.hamming(2*N+1)
w, h = signal.freqz(b)
h = h*np.exp(1j*w*N)
import matplotlib.pyplot as plt
fig = plt.figure()
plt.title('Digital filter frequency response')
ax1 = fig.add_subplot(111)
plt.plot(w, 20 * np.log10(abs(h)), 'b')
plt.ylabel('Amplitude [dB]', color='b')
plt.xlabel('Frequency [rad/sample]')
ax2 = ax1.twinx()
angles = np.angle(h)
plt.plot(w, angles, 'g')
plt.ylabel('Angle (radians)', color='g')
plt.grid()
plt.axis('tight')
plt.show()

Most of the above is from Scipy documentation which you link to. This gets plotted:

Figure 1. Frequency response of the Hilbert transformer centered at time zero.

The phase goes to zero at frequency zero, because the phase frequency response transitions from $-\pi/2$ for positive frequencies to $\pi/2$ for negative frequencies.


How Does the RMS of White Noise Change with Sampling Frequency? There is an analog system which includes the continuous-time linear equalizer (CTLE). With some .noise analysis the power-spectral density (PSD) of the noise in that system is provided. So let's not care about the details how PSD is provided, the only thing that matters is that I have access to the noise PSD of the circuit. The noise PSD is not necessarily white and it can have an arbitrary shape, but it dies out after a certain frequency (let's say most energy is contained within 100 GHz).
By having the noise PSD, I integrate it over a large enough frequency range (in this example from 0 to 100 GHz) to get the (almost) total power of the noise, let's name it $\sigma^2$.

Now, in my MATLAB simulation I have a clean sampled signal $s(n \,\,\Delta t)$ with sampling frequency $f_s$, where $f_s = 1 / \Delta t$ (To add more details, $s(t)$ is a continuous-time signal that has a bandwidth of 30GHz and I am oversampling $s(t)$ by taking 32 samples per symbol time. So the sampling frequency in this example would be $ f_s $ = 32 x 30GHz = 960 GHz . Maybe these details don't matter).

I want to add discrete white Gaussian noise to $s(n \,\,\Delta t)$ such that this discrete Gaussian noise has the same power (i.e. variance) as the continuous (non-white) noise which was calculated to be $\sigma^2$. I know I can do something like $N = \sigma$ x randn( size(s) ) in MATLAB to generate noise with power $\sigma^2$, but my main question is that what is the dependency of the discrete noise in this example to the sampling rate $f_s$. If I change $f_s$ how the randn() function or the $\sigma$ should change.

Cause I see in some pages it is mentioned that:

"Given a continuous White Noise (Wide Sense) with variance $\sigma^2$ and you want sample it at rate of $f_s$ you should generate discrete noise samples with variance of $f_s \, \sigma^2$",

but I don't know what is the mathematical reasoning behind it and if the same applies in my case or not.

If you have any good read on this matter, that would be really appreciated.

It changes according to the bandwidth of the LPF Filter in the sampling chain.
Namely if you use the same LPF yet with two different samplings rates the RMS of the noise is the same.

Now, the theoretical RMS of White Noise is given by its PSD which is a constant.
The RMS of the sampled noise is the constant value multiplied by the BW of the LPF in the sampling chain.


Non-coherent receiver for DBPSK: Implementation details I am implementing a non-coherent receiver for DBPSK.
I have a computer science background with little knowledge of DSP so my questions here may seem obvious.
The design of the receiver I am trying to implement is the following:
Design of DPSK receiver

And what I am doing right now is:

On the tx side

Differentially encoding my bit sequence
BPSK modulating it
Applying a high-pass filter

On the rx side

Applying a band-pass filter around the carrier frequency
Delaying (by one bit worth of data) and multiplying the data by itself
Applying a low-pass filter
Integrating over one bit worth of data (i.e. sum the sample values)
Making the decision depending on the sign of my integral

By 'data' I mean samples at 44.1kHz.
When I simulate this on Octave, it works (i.e. I can detect the bit sequence even with high AWGN).
The problem is I am trying to move to an infinitely long data sequence (i.e. I do not know the exact moment where the tx began).
So here are my questions:

How can I achieve time and frame synchronization ? (I tried with a barker code but I am not sure I am getting it right)
This is the thing that's most confusing to me. My guess is you'd have to get time synced before the integration block, and frame synced after the decision block.
How can I express the bandwidth in function of carrier frequency f and bit rate ?
Does the relation between bit rate and carrier frequency have any influence on the receiver ?
What are the optimal values for the filters ? (number of taps)
How do I calculate Eb/N0 for in my case ? (I want to compare my BER curve with others to see how good is my receiver)

1) A typical approach would be to have a frame alignment word (and by "word" I mean bit sequence) at the beginning of every frame.  You could time synchronize by looking for that FAW through cross-correlations.  Barker codes are a good choice because they have good autocorrelation properties.

The problem with the approach is that the data, being random, can look like the FAW sometimes.  You can use the periodicity of the FAW to screen out those false positives.  For instance, if your frame is 1000 bits long, you would look for multiple hits spaced 1000 bits apart.  Once you have found them you have both time and frame synchronization.

2) There is not a direct relationship between carrier frequency and bandwidth.  There is a direct relationship between bit rate and bandwidth, but there are other factors involved.  The easiest way to determine your bandwidth is to take an FFT of your signal and measure it.

3) Your question is unclear.

4) You haven't provided enough information, and "optimal" can mean a lot of different things.  If you want a good answer to this you should start a new thread with more details.

5) It sounds like you are adding noise via AWGN, so you know what your SNR is.  With BPSK Eb/No is equal to the SNR.


Deriving the CRLB for an active system with non-square shaped bandwidth I'm attempting to calculate the CRLB for a bandlimited time-delay system which has a triangular shaped signal spectrum, instead of the usual square one.

Currently I'm working on an active system.

Following the procedures done by Quazi et al. in  "An overview on the time delay estimate in active and passive systems for target localization",

I can reach two different results. The first one, using his calculations for a High SNR, active system, don't really match my data, while the second one, computing for a Passive system and considering a $2$ factor (as he explains after equation 23) due to only one of the sources being contaminated by noise (which makes intuitive sense to me, as long as the SNR is high), actually fits my data very well.

As such, I am assuming I'm doing the first set of calculations wrong, and I was hoping some of you could help me spot my mistake.

So here goes, for both the calculations

PASSIVE

$$\frac{1}{\sigma^2} = 2T \int_{0}^{B}(2\pi f)^2 \frac{S(f)}{2N(f)}df$$

where my noise is square shaped, and my signal is triangle shaped, in baseband, so

$$\frac{1}{\sigma^2} = 2T \int_{0}^{B}(2\pi f)^2 \frac{S_0(1-f/B)}{2N_0}df$$

$$\frac{1}{\sigma^2} = 2T(2\pi)^2 \frac{S_0}{2N_0}\int_{0}^{B} f^2(1-f/B) df$$
$$\frac{1}{\sigma^2} = T \ 4\pi^2 \  SNR \  (B^3/3-B^3/4) $$
$$ \sigma^2= \frac{12}{4\pi^2T} \frac{1}{SNR} \frac{1}{B^3}$$

which, for an active system, dividing the variance by two since the source isn't noisy, I end up with
$$ \sigma^2= \frac{12}{8\pi^2T} \frac{1}{SNR} \frac{1}{B^3}$$

This fits very well to my data, and the last argument is logical to me.

ACTIVE

So, he starts by defining the CRLB as $\sigma^2 = \frac{1}{d^2\beta^2}$ where

$$d^2 = \frac{2E}{N0}$$
where $E$ is the energy of the signal, so $E = ST$, where $T$ is the signal observation time

$$\beta^2 = \frac{\int_{-\infty}^{+\infty}\omega^2|F(\omega)|^2 d\omega}{\int_{-\infty}^{+\infty}|F(\omega)|^2 d\omega}$$

Now, the calculation of $\beta$ is similar to what I've done before, setting my $S(f) = S_0(1-\frac{f}{B})$

$$\beta^2 = \frac{\int_{0}^{B}(2\pi f)^2S_0(1-\frac{f}{B}) \ \ 2\pi df}{\int_{0}^{B} S_0(1-\frac{f}{B}) \ \  2\pi df}$$

$$\beta^2 = \frac{(2\pi)^2\int_{0}^{B} f^2(1-\frac{f}{B}) \ \ df}{\int_{0}^{B}
(1-\frac{f}{B}) \ \ df}$$

$$\beta^2 = \frac{(2\pi)^2\frac{B^3}{12}}{ \frac{B}{2}} = (2\pi)^2\frac{B^2}{6}$$

solving now for the cramer-rao

$$\frac{1}{\sigma^2}=\frac{2E}{N_0}(2\pi)^2\frac{B^2}{6}$$

multiply top and bottom by $B$ so $N = N_0B$, and $S/N = SNR$

$$\frac{1}{\sigma^2}=2T(2\pi)^2 \ SNR\  \frac{B^3}{6} = 8\pi^2T \ SNR \frac{B^3}{6}$$

$$ \sigma^2= \frac{6}{8\pi^2T} \frac{1}{SNR} \frac{1}{B^3}$$

Which does not fit my data as tightly.

I'd like to ask you where my intuition is going wrong with both solutions, or if I am making any dumb mistakes.

Thanks

Found the solution, I believe.

It's a dumb mistake. I'm calculating the SNR as $S_0 / N_0$ in one case and as $S/N$ in the other.

Since the noise is square shaped, and the signal is triangle shaped, the first is obviously wrong.
Correcting to $\frac{S_0 B/2}{N_0 B}$ gives the factor that was missing from one to the other.


FFT - mixed radix - bit reversal I wonder what is bit reversal for mixed radix FFT. Is there any algorithm that compute the bit reversal for various mixed radix FFT?
Here should be mention why I want that. For example if I have signal with N=48 points. There is some method (which I still understand very little) to make DFT by prepare matrix 2x2x2x2x3, and make 5 dimensional DFT calculations. And because I can’t understand that idea but I understand radix-2 and radix-4 I wonder if I can make the matrix like 3x16 and on each 16 point row make the radix-2 or radix for. But not sure how to prepare bit reversal for input data (or for output if I make decimation in frequency).

3 by 16 or 16 by 3 or 4 by 12 or 12 by 4

Hope this helps, not the mixed case , but I think this is the right direction

I don't have a copy of
Brigham, E. Oran, and E. Oran Brigham. The fast Fourier transform and its applications. Vol. 1. Englewood Cliffs, NJ: prentice Hall, 1988.

handy. I have a vague recollection that he had some mixed fft examples but I think was radix 4.  If you go to the fftw website and look at the papers section, there might be more help.

x=reshape(1:25',5,5)'  % read in row order fortran index
x =
     1     2     3     4     5
     6     7     8     9    10
    11    12    13    14    15
    16    17    18    19    20
    21    22    23    24    25
reshape(x,25,1)        % read out column order
ans =
     1
     6
    11
    16
    21
     2
     7
    12
    17
    22
     3
     8
    13
    18
    23
     4
     9
    14
    19
    24
     5
    10
    15
    20
    25

[y,i]=digitrevorder(0:24,5)  % c indexing
y =
  Columns 1 through 20
     0     5    10    15    20     1     6    11    16    21     2     7    12    17    22     3     8    13    18    23
  Columns 21 through 25
     4     9    14    19    24
i =
     1
     6
    11
    16
    21
     2
     7
    12
    17
    22
     3
     8
    13
    18
    23
     4
     9
    14
    19
    24
     5
    10
    15
    20
    25
echo off