How to use matlab toolbox for 2qam dfe equalizer? I made Channel equalization with matlab,
You can find the code below.

clc; clear all; N = 1e3; MQAM = 2;

d = randi([0 1],N,1);
x = qammod(d,MQAM,'gray','InputType','bit');
h = [0.1 0.5, 0.1 0.2];
r = awgn(complex(filter(h, 1, x), zeros(size(x))), 10, 'measured');

eq = comm.DecisionFeedbackEqualizer;
eq.Constellation = qammod(0:1, 2);
eq.StepSize = maxstep(eq,r);
eq.ReferenceTap = 2;

[y,e,w] = eq(r, x);

dec = qamdemod(y,MQAM,'gray','OutputType','bit');

errors = sum( dec~= d)

For some reasons it is not working and I don't know why the error is still the half of the transmitted symbols number?

I think it is the delay of the equalizer that I am not understanding ... somehow the outputs are delayed and I need to find the correct delay

From the MATLAB documentation (https://uk.mathworks.com/help/comm/ref/comm.decisionfeedbackequalizer-system-object.html): "The equalizer uses the reference tap location to track the main energy of the channel". The equalizer object has a property called the reference tap which has default value of one. Looking at your channel coefficients, it looks like the tap with the main energy is the second tap so you should set eq.ReferenceTap = 2.

You can also work on choosing an appropriate step size. The equalizer returns the error signal back, you call it e, you can plot that and use it as a guide to pick a step size. You can also use the in-built function maxstep which gives you the maximum step size needed for a LMS-like algorithm to converge.

Last suggestion I have is unrelated but when calling awgn and the input is real, then the function does not add complex noise so your noise will only appear in the in-phase/real axis. You can get around it just by creating your own complex noise, or by specifying the input as complex with imaginary part equal to zero: r = awgn(complex(filter(h, 1, x), zeros(size(x))), SNR_dB, 'measured').

Edit

The OP is still having issues so I will add the settings that I use to get results.

eq = comm.DecisionFeedbackEqualizer;
eq.Constellation = qammod(0:MQAM-1, MQAM);
[~, eq.ReferenceTap] = max(h);
eq.StepSize = maxstep(eq, r)/2; % maxstep returns largest so choose something smaller to be safe

And here are the resulting error signal and before/after constellation plot.


Can discrete wavelet transform for denoising purposes be implemented in real time? I have been researching effective algorithms for denoising biomedical signals (non-stationary) that can be implemented in real time either using FPGA or DSP. I can across many suggestions for algorithms and effectiveness but found DWT to be the one that can denoise most of the noises I am interested in (power line interference, surge, power dips..)
I used the wavelet toolbox in matlab and performed denoising and wanted to make sure that this method had been implemented effective before as real-time.
I read few papers too but want more opinions.

Even though I'm not much of a wavelet expert, I can answer with a definitive YES!!
Now, do you have any messy details to complain about, like the size of the box, the power that it consumes, or the number of frames of delay between the taking of the image and its display?  Because that's your real problem.
Take the particular wavelet transform you want to use.  Sit down and with pencil and paper figure out how many calculations (floating point operations, basically multiplies and adds) that need to be performed per frame.  Look to see if the transform you want to use works with the fast wavelet transform.  Multiply the FLOP per frame by your frame rate -- you should have an answer in MFLOPS or GFLOPS.  Now hold that up against available processors or FPGAs, and figure out what needs to be in your box.  Then multiply by three, because it always takes more real estate to do a calculation than you think.


Coherence bandwidth/delay spread in OFDM Say we want to design an OFDM system and we have some parameters given to us, i.e. the coherence time $T_c$ and the max delay spread $T_m$.

It's clear to me why we need to ensure $T_s(N+N_{cp}) \ll T_c$ where $N_{cp}$ is the cyclic prefix in samples - we want the channel to stay coherent during the entire symbol. What is not clear is what happens to the constraint $T_s(N+N_{cp}) \gg T_m \Leftrightarrow B_N \ll B_c$. Apparently, in OFDM it suffices with $T_s N_{cp} \geq T_m$.

It makes intuitive sense in a way, but I'm still not sure exactly how the coherence bandwidth interacts with the subchannel bandwidths. Why don't we need to ensure $B_N \ll B_c$?

Indeed, the constraint $T_s(N+N_{cp}) \gg T_m$ should be derived as $T_s(N+N_{cp}) \gg T_s N_{cp} \geq T_m$.

The first inequality to make sure that the CP overhead is as small as possible. The second inquality is what you have said, to avoid ISI between OFDM symbols and to convert linear convolution to circular convolution so that we can use single-tap equalizer.

Thus, $B_N << B_c$ is a consequence of the prior design criterion. However, if subcarrier spacing $B_N$ is much smaller than coherence bandwidth $B_c$, the channel estimation is facilitated by adding several pilots in $B_c$. For info, in LTE, two pilot are separated by 45 kHz while the rms delay spread is about 1us gives $B_{c,50\%} \approx 200kHz$ and $B_{c,90\%} \approx 20kHz$; thus 45kHz pilot spacing is good compromise between pilot overhead and estimation


Time-shift confusion Say input-output of a system is defined as:

$$
x[n] \longrightarrow x[nM]
$$

then what will be the output of $x[n-1]$?

will it be:

\begin{align}
x[n-1] \longrightarrow  &x\left[(n-1)M\right] = x\left[nM - M\right]\\
&\textrm{or}\\
x[n-1] \longrightarrow  &x\left[nM - 1\right]
\end{align}

And if someone says shift a signal by 1, if the signal is $x[2n-1]$ then the result will be

$$x[2n-1-1]\quad \textrm{or}\quad x[2(n-1)-1]\quad ?$$

I'm really confused between these. Can anyone make me understand it, please?

Perhaps the easiest way to understand it is by seeing what happens as we plot graphs of an arbitrary sequence $x[n]$:

As it goes through the system and take one every second sample (i.e. $x[0]$, $x[2]$, $x[4]$, $x[6]$, $x[8]$ and $x[10]$ in the graph above), you get $y[n] = x[nM]$ which can be depicted as:

Now when you shift the input you get $x'[n] = x[n-1]$:

If we again take one out of every second sample we would then get $y'[n] = x'[nM]$:

We then note that this corresponds to the sequence $x[1]$, $x[3]$, $x[5]$, $x[7]$ and $x[9]$ from the original graph. For this simple case where $M=2$ you can see that the result is
$$
\begin{align*}
  y'[n] &= x'[nM] \\
        &= x[2n-1] \\
        & =x[nM-1] &\mbox{given $M=2$}
\end{align*}
$$

A similar process could be followed to arrive at the general result whereby the output of that system given an input of $x[n-1]$ would be $x[nM-1]$.


Is there a counterpart for zero forcing in communications in DSP? If yes, what is forced to zero?

If you understand zero forcing from the perspective of vector spaces, you would be able to identify it different systems including DSP. In essence you are projecting your equivalent (received vector in communications) in the null space of a specific subspace which makes signals lying in that subspace to be zeroed. In communications, this would be the null space of interfering users (in a multi user scenario).

What's more, you can extend this idea from vector space to Hilbert spaces. For ex: consider the vector/Hilbert space of N degree polynomial, now you might want to project your polynomial on the null space  of the subspace spanned by some polynomials which are not of your interest.

For ex: you might want to approximate a N degree polynomial by a M degree polynomial, $M<N$, then you can find the null space of the polynomials from degree $M+1$ uptil N (using gram Schmidt procedure if possible) and then project the N degree polynomial along that to zero force those terms. So you see the idea is generic. This last example of polynomial approximation I quoted is also an example of zero forcing in DSP, less trivial one though.

Mbaz has quoted a reference for a more trivial explanation.


Strange phase spectrum from stft I've tried to implement STFT in functional way in Python. I've got code such as below:

def genSpectrums_dft(x, w, N, dft_func = fp.fft):
    hN = N / 2
    M = len(x)
    hM1 = int(math.floor((M + 1) / 2))
    hM2 = int(math.floor(M / 2))
    wM = w(M)
    w_sum = sum(wM)
    xw = x * wM / w_sum
    dft_buffer = np.concatenate(
        (xw[hM2:], np.zeros(N - hM1 - hM2), xw[:hM2])
    )
    X = dft_func(dft_buffer)
    mX = 20 * np.log10(abs(X[:hN]))
    pX = np.unwrap(np.angle(X[:hN]))
    return mX, pX

def genSpectrums_stft(x, w, M, N, H, dft_func = fp.fft):
    hM1 = int(math.floor((M + 1) / 2))
    hM2 = int(math.floor(M / 2))
    x_ext = np.concatenate((np.zeros(hM2), x, np.zeros(hM2)))
    lv = np.arange(len(x) / H + 1)
    wM = w(M)
    w_sum = sum(wM)
    xlv = np.array([x_ext[l * H:l * H + M] / w_sum for l in lv])
    return tuple(
        np.transpose(np.array(X))
        for X in zip(
            *[genSpectrums_dft(xl, w, N, dft_func) for xl in xlv]
        )
    )

Full working code is available here (Python, numpy, scipy and matplotlib required to run it).

When I run my code for any wave, STFT magnitude spectrum seems to be quite good, but STFT phase spectrum always looks like this:

But as far as I know it should rather look like this:

Both these spectrums have been generated from the same wave and the second parameters. The second one has less frequencies. Why my phase spectrum has these strange vertical lines?

I've got a few remarks rewarding my code:

I am almost positive that finding dft spectrums works correctly
Though the spectrums distinct, when I inverse the process, I get something really similar to the input

If you are not familiar with functional programming, I don't mind if you correct my mistakes in structural way, but the most important thing for me is to understand the problem with my code.

You are unwrapping the phase along the wrong axis. You are obviously interested in the phase change of the frequency bins along the time axis in order to say something about the local phase coherence.


Trying to write a signal "unsmoothing" algorithm Up front: I'm just a hobbyist, and I know zero about signal processing; I'm trying to write a simple(?) javascript algorithm and hoping I don't need a PhD in DSP to make progress on it.

I have an input signal that often shows only the tiniest difference between the max value and the min value. I'm trying to figure out a way to convert the values such that the max is proportionally much larger than the min, with a steep curve in between. I've tried everything I can think of with exponents and logarithms, even just subtracting the min from the rest of the curve, but nothing seems to work. Part of the problem, as can be seen in my awesome diagrams below, is that the signal strength varies quite a bit; I don't know how to account for that. Also, sometimes the signal is already pretty spikey, and I don't know how to account for that either.

I've poked around with google, but I don't really know what questions to ask. I'm not sure whether there is a class of algorithms that does this sort of thing. I tried "signal unsmoothing", with no luck, and I don't know what else it could be called.

Some hand-drawn diagrams, placeholders while I try to learn how to draw curves on stackexchange. The two on the left show the kinds of signals I often get: varying strength, but the max is barely bigger than the min. The one on the right is what I'm trying to get: a big spike in the middle with a steep curve down to the min.

Is there a name for a class of algorithms that does this sort of thing? Is there a better way to formulate my questions so I can google more effectively? Is there a ridiculously trivial way to do this?

In case it helps, my ultimate goal: I'm modeling the movement of mosquitos. I'm working on getting them to move toward the smell of their prey while accounting for negative factors such as toxins in the air, or high temperature, as well as their jittery motions to avoid being squashed. It's going pretty well, except I have this chronic problem of the input signals being too flat for meaningful comparisons.

What you are trying to do is Contrast Enhancement. Contrast is the ratio of maximum to the minimum.

Lets consider as example that your signal's intensity is ranging from 30 to 40 but you need a range of 0 to 255 with good contrast.

Lets consider your 1D signal as I(x). First find the minimum of the signal and subtract min from the whole curve. In our example 30 is the minimum.So subtract 30 from the whole curve.The signal will be ranging from 0 to 10 now.

  J(x) = I(x) - Min(I(x))

But we want the intensity of the signal from 0 to 255. So find the ratio between the maximum value allowed(255 in our case) to maximum of your signal(our signal max is 10 now) and multiply the whole curve with that factor. Now our signal's intensity will be ranging from 0 to 255.

  K(x) = J(x)*(Max_value_allowed/Max(J(x))

At this point you should have got the spikey curve you need. But if you feel the curve is not spikey enough, contrast can be increased by applying negative slope to intensities in first half of histogram(In our case 0 to 127) and positive slope to intensities in second half  of histogram(128 to 255 in our case) which will make lower intensities much lower and the high intensities much higher.Good luck.Hope this helps you.

Added:Instead of above step3 you could also do,

  L(x) = power(K(x),n)/power(Max(K(x)),n-1)

n should be greater than 1.As n increases the spikiness also increases. In your case choosing around 2 or 3 should be enough I guess.


AM-AM and AM-PM plots I have complex baseband signal and its output when passed through power amplifier. How to plot the AM-AM and AM-PM plots?

The AM-AM and AM-PM plots are based on a model that assumes the input to the power amplifier is a single unmodulated tone.  According to the model, the output of the power amplifier is a tone of the same frequency with a gain and phase shift that are functions of the input amplitude.  If the input is $A \exp{(j(\omega t+\theta))}$, the output is $g_A(|A|)A \exp{(j(\omega t+\theta+g_\theta(|A|)))}$.  Here $g_A()$ and $g_\theta()$ are functions that compute the gain and phase shift as a function of the input magnitude.

Therefore, in order to build an AM-AM and AM-PM model, you need to input a tone into the system and sweep its power over a wide range, measuring the gain and phase shift at each input power.  When the drive level is low enough that the amplifier operation is linear, you should find that the gain and phase shifts stop changing as the input power changes.  At higher drive levels, the gain and phase shifts will change substantially as a function of the input drive level.  You need to make measurements that span the drive levels at which you might operate the power amplifier, down to the level at which the amplifier becomes linear.

Although power amplifier AM-AM and AM-PM curves are typically measured using tones, you could use another waveform instead.  However, developing a model using a waveform with a time-varying envelope is a more open-ended problem, because the input amplitude varies both with time and with drive level.

Keep in mind that the AM-AM and AM-PM curves represent a simplified model.  That simplicity makes it easy to use, but also limits its accuracy.


Relation between probability of a signal and its magnitude spectrum I am a newbie. In these days I am searching if there are reference about properties that can relate the probability of the signal and the magnitude spectrum. I know what are both. But I was curious wheter there are connection between them. Can you help me
Given a signal y(t), the magnitude spectrum  is
$$Y(f)=\left\lvert\int y(t)e^{-2\pi ft}\,\mathrm dt\right\rvert,$$
in other words the absolute value of the FT of the signal.

Given a signal y(t), the magnitude spectrum  is
$$Y(f)=\left\lvert\int y(t)e^{-2\pi ft}\,\mathrm dt\right\rvert,$$
in other words the absolute value of the FT of the signal.

But that's only applicable for deterministic signals¹. If your signal itself is random, then you can't know $y(t)$ – you only know some stochastic properties of it. There's no $y(t)$ to transform if you don't actually know $y(t)$. For random signals, $y(t)$ can't be given, otherwise the signal wouldn't be random. You might observe a realization, but that realization doesn't describe the full stochastic process!
The definition of the power spectral density for random signals like your $y$ is the Fourier transform of the autocorrelation function. For that to be sensibly defined, you need your autocorrelation function to depend only on one variable, the time shift, which means your random signal needs to be wide-sense stationary. But, if that's the case, then the autocorrelation function (ACF) is defined as
$$\phi_{yy}(\tau) := \mathbb E \left\{ y(t) y^*(t+\tau) \right\},$$
and then the power spectral density, which gives your magnitude spectrum is its Fourier transform, $\Phi_{yy}(f) :=  \int_{-\infty}^{\infty} \phi_{yy}(\tau) e^{-i2\pi f\tau}\,\mathrm d\tau$.
Now, in the expectation $\mathbb E \left\{ y(t) y^*(t+\tau) \right\}$, there's your random signals probability density function:
\begin{align}
\phi_{yy}(\tau) &:=
\mathbb E \left\{ y(t) y^*(t+\tau) \right\} \\
&=
\int_{-\infty}^\infty
  \int_{-\infty}^\infty
    y'_1\,y'_2
    \,\cdot\,
    f_{
      (y({t_1}),y({t_1+\tau}))
    }
   (y'_1,y'_2)
    \,
  \mathrm d y'_2
  \;
\mathrm d y'_1
\end{align}
Notice that here, the density is a function of two values, and parameterizes over the time shift.

¹ you're also missing imaginary unit in the exponent and the integration bounds, but that's not the problem here. But, really, make sure to be able to write down the Fourier transform correctly.


DIY DSP projects ideas I have a Masters Degree in DSP but all I learnt there is theory. I would like to complete some projects at home which I can add in my resume to showcase my skills and to get jobs in industry. But I have not found anything after googling. So can someone here please guide me what should I do?
P.S.:

I am not interested in image processing.

I would like to start projects on Matlab/Octave or Python and then writing the same code in Embedded C for a DSP hardware.

There's an unlimited amount of projects you can do... The extend depends on how much you want to (or you can) turn your home into a lab ?
Audio projects can range from software synthesizers, multimedia players, to aduio compression algorithms, micophone array based advanced imaging, and to all sorts of acoustic magic. You would need alot of electronics and acoustic skills, in addition to DSP, for these projects though.
Hardware projects will also cost extra bucks & mess, compared to pure software projects...
Also, eventhough you've got your masters degree, just still need to polish your DSP skills, and make sure that you don't miss any serious topics in the fundamentals. The best way to do is to code the concepts in a book such as Discrete-Time Signal Processing. Implement (or verify) various algorithms and theorems there. Pay attention to convolutions, LCCDE recursions, filter design, sampling, DFT, and FFT topics. Write your own functions, simulations and libraries. Make theoretical predictions, and try verifying them with simulations... Use Matlab or Octave at this stage.
Have you done some Adaptive Filtering? Why not implement them in realistic applications? Try channel equalisation, echo cancellations, noise reduction, system identification, line enhancement, prediction, compression etc. Try LMS, RLS or their variants.
In the low-level software side, try improving your C / C++ and Assembler programming-debugging skills. Chose platforms such as Intel/AMD or ARM-based CPUs. I bet you will benefit from learning how to program MMX, SSE, AVX SIMD-vector extensions of these CPUs. And get a general glimpse on compilers, switches, libraries, open source development, various tricks and hardware properties and limitations of PCs.
Learn some operating system programming too. Windows,Linux,OS, and lately but more attractively Mobile Android. Learn how to use basic compilers for these platforms and their working mechanisms.
Start with the simplest and cleanest hardware: Arduino. Learn how to code it in its IDE. That's the most user friendly mCU ever. Though, it's a limited one, but you can still do a lot with it. It has an ADC-DAC; you may wish to try sampling some sources and processing them. Best is to go with a simple microphone shield (an attachement board to an Arduino). Record audio, process it, and output some signals based on it. You will learn a lot about USB connection, Serial (RS-232) communication, SPI, I2C communication protocols, motor drive examples, various sensors involving 9DOF acceleromoter-gyro-magnetic sensors, distance measurement sensors, light measuring sensors, temperature, pressure, gas sensors... You will have a lot to do with them... You may also learn practical Kalman filtering for perfecting your measurement results too.
You should support your Arduino with some basic external hardware, such as buffers, filters and actuators. It's good to learn connecting/interfacing them to the mcu, for practical embedded DSP.
Expand your hardware with a Raspberry PI, the most friendly electrically interactable PC platform with ARM based multi-core CPU and a Linux based OS in it.
Add some simple programmable SDR-RTL (software defined radio) kits to your Raspberry Pi to turn it into an RF receiver and data processing machine. You can receive signals from radio stations, Satellites, AirPlanes, wireless phones or any RF source in its range! You can write C-code to implement your analog/digital decoding algorithm. That makes a set of plain advanced DSP projects.
You can also progress with some more industrial MCUs from ST, NXP, Texas, ADI or any similar vendor. They are producing ARM complient cores. HEnce once you learn ARM architecture you will be confident in all of these platforms.
One good thing is, Arduino IDE also supports an STM32 based development card just as it supports ative Arduino UNO cards. So You can transit into ST programming from Arduino IDE and go more professional with ST native development later.
You can further add ADI BlackFin, SHArC, or Texas TMS based DSP processors.
Furthermore, last decade have withnessed great progress in Big-Data processing using GPU computing; so you will benefit form learning some parallel programming using Nvidia CUDA or openCL...
All these will take some time and patience.


What are some methods for detecting a sinusoid in noise? I'm looking for a broad classification on the kinds of techniques available. Something I can use to begin a literature survey.

Some details:

should be feasible for real-time implementation.
Binary decision, I should be able to input the signal continuously, and the algo has to detect the start and stop times of the sinusoid in the waveform.
There is no other signal, i.e. if the sinusoid is not present, then there will only be noise
the input is band-limited, and the sinusoid, if present, is assured to be within that band.
Tradeoffs are between speed (how soon after the appearance of the sinusoid can the algo detect its presence) and false positives (should be minimal)
I can't provide any exact numbers about acceptable performance measures, because I'm not sure myself. I intend to implement all your suggestions for my application and find out myself. I'm just looking for the "standard" techniques for tackling this problem.

Further details:

The input is the output of a band pass filter, so the noise is
also significant only in the pass band.
There is no surety when
the sinusoid may appear. The duration of the sinusoid is in the
range of 50-100 ms. The amplitudes of the sinusoid will fluctuate.

Depending on the frequency of sampling, a FFT (Fast Fourier Transform) will work. For example, if your sampling rate is only once every 20 ms, you'll only get a few samples of the sinusoid, but if you're sampling every 0.5 ms then you'd get a lot more samples. FFTs usually work best with a large number of samples to average over. In that case even if your signal is the sum of a few sinusoids, they can be accurately determined.

Alternatively, you could have a look at the MUSIC algorithm. I'm not too sure about the details of how it is implemented, but it has been implemented in several real-time detection scenarios. An alternative to MUSIC is the Esprit algorithm.

Still, if your sample size is large enough and the noise doesn't swamp the signal entirely, then an FFT will be (on average) the fastest option.


Determine if system is linear time variant The system equation is given as:
$$y(n)=(n-1)x(n-1)+(n+1)x(n+1)$$
I solved that the system is time variant:
\begin{align}
y(n-k)&=(n-k-1)x(n-k-1)+(n-k+1)x(n-k+1)\\
H[x(n-k)]&=(n-1)x(n-k-1)+(n+1)x(n-k+1)
\end{align}
And for the linear/non-linear part, here is what I have so far:
\begin{align}
H\big[a_1x_1(n)+a_2x_2(n)\big]&=(n-1)\big[a_1x_1(n-1)+a_2x_2(n-1)\big]\\&+(n+1)\big[a_1x_1(n+1)+a_2x_2(n+2)\big]
\end{align}
It is ok?

To be time invariant you would have to show that:
$$y(n-k) = (n-1)x(n-k-1) + (n+1)x(n-k+1)$$
(You shift the $x$ and you get a shift in $y$).
But you can't get rid of the $(n-k-1)$ and $(n-k+1)$.
In general any time you have the time term $n$ by itself it will be time variant unless it cancels out somehow.


Can't seem to get welch function in scipy work correctly I'm trying to learn DSP through using Python's scipy package.
I have measure some 200Hz signal from a machine. Now I want to inspect the spectral density of the signal. Here's the signal plotted out:

As you can see, the signal is obviously not a low frequency signal. However, I still get the PSD graph like below.
Here's the code I have been using:

def plot_and_save(sample,i=""):
filepath = sample + "/" + "merged_data.csv"
df = pd.read_csv(filepath)

sns.set_style("darkgrid")
sns.set_colodef plot_and_save(sample,i=""):r_codes("dark")
fig1, axes = plt.subplots(3,2)
fig1.figsize = [1920,1080]

axes[0,0].plot(df["Time"], df["Voltage"], marker="",color="blue")
axes[0,0].set_title("Plot of Voltage against time")
axes[0,0].set_xlabel("Time (us)", ha= "right", x=1.0)
axes[0,0].set_ylabel("Voltage (V)")

axes[0,1].plot(df["Time"], df["Current"], marker="",color="red")
axes[0,1].set_title("Plot of Current against time")
axes[0,1].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[0,1].set_ylabel("Current (A)")

sampling_freq = 1e6
vol_f, vol_psd = welch(df["Voltage"].to_numpy(), fs=sampling_freq)
axes[1,0].plot(vol_f, vol_psd, marker="",color="green")
axes[1,0].set_title("Plot of Voltage's Power Spectral Density")
axes[1,0].set_xlabel("Frequency", ha= "right",x=1.0)
axfilepath = sample + "/" + "merged_data.csv"
df = pd.read_csv(filepath)

sns.set_style("darkgrid")
sns.set_color_codes("dark")
fig1, axes = plt.subplots(3,2)
fig1.figsize = [1920,1080]

axes[0,0].plot(df["Time"], df["Voltage"], marker="",color="blue")
axes[0,0].set_title("Plot of Voltage against time")
axes[0,0].set_xlabel("Time (us)", ha= "right", x=1.0)
axes[0,0].set_ylabel("Voltage (V)")

axes[0,1].plot(df["Time"], df["Current"], marker="",color="red")
axes[0,1].set_title("Plot of Current against time")
axes[0,1].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[0,1].set_ylabel("Current (A)")

sampling_freq = 1e6
vol_f, vol_psd = welch(df["Voltage"].to_numpy(), fs=sampling_freq)
axes[1,0].plot(vol_f, vol_psd, marker="",color="green")
axes[1,0].set_title("Plot of Voltage's Power Spectral Density")
axes[1,0].set_xlabel("Frequency", ha= "right",x=1.0)
axes[1,0].set_ylabel("PSD")

cur_f, cur_psd = welch(df["Current"], fs=sampling_freq)
axes[1,1].plot(cur_f, cur_psd, marker="",color="yellow")
axes[1,1].set_title("Plot of Current's Power Spectral Density")
axes[1,1].set_xlabel("Frequency", ha= "right",x=1.0)
axes[1,1].set_ylabel("PSD")

axes[2,0].plot(df["Time"], df["Current"]*df["Voltage"], marker="",color="red")def plot_and_save(sample,i=""):
axes[2,0].set_title("Plot of power against time")
axes[2,0].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[2,0].set_ylabel("Watt")

sns.scatterplot(x="Current", y="Voltage", data=df, ax=axes[2,1], color = 'r', edgecolor="None", linewidth=0)
axes[2,1].set_title("Plot of Voltage against Current")
axes[2,1].set_xlabel("Current (A)", ha="right", x=1.0)
axes[2,1].set_ylabel("Voltage (V)")

fig1.suptitle("Sample number: " + str(i+5),fontsize=16)
# plt.savefig("formatted/"+date+"/figure/"+str(i+5)def plot_and_save(sample,i=""):".svg",format="svg")
plt.show()es[1,0].set_ylabel("PSD")

cur_f, cur_psd = welch(df["Current"], fs=sampling_freq)
axes[1,1].plot(cur_f, cur_psd, marker="",color="yellow")
axes[1,1].set_title("Plot of Current's Power Spectral Density")
axes[1,1].set_xlabel("Frequency", ha= "right",x=1.0)
axes[1,1].set_ylabel("PSD")

axes[2,0].plot(df["Time"], df["Current"]*df["Voltage"], marker="",color="red")
axes[2,0].set_title("Plot of power against time")
axes[2,0].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[2,0].set_ylabel("Watt")

sns.scatterplot(x="Current", y="Voltage", data=df, ax=axes[2,1], color = 'r', edgecolor="None", linewidth=0)
axes[2,1].set_title("Plot of Voltage against Current")
axes[2,1].set_xlabel("Current (A)", ha="right", x=1.0)
axes[2,1].set_ylabel("Voltage (V)")

fig1.suptitle("Sample number: " + str(i+5),fontsize=16)
# plt.savefig("formatted/"+date+"/figure/"+str(i+5)".svg",format="svg")
plt.show()

Thank you so much for the help!

I revised my code and realised the result was correct the whole time. The reason is that not plotting on a log scale, I thought the frequency with the most power density is at 0Hz (it should be 200Hz)!


Interpreting microphone output signal I have a question about interpretation of the signal produced by a microphone.

The output of a microphone is a voltage level that corresponds to the pressure on the membrane. The part that confuses me is why are the output signal values located above and below the average (picture attached). Do the values "above and below" correspond to the membrane deflecting, and then going back to its original position? If that's the case, should I just concider the values above the average for my further computation of the Fft?

Thank you!

Sound behaves (more or less) like a wave. When the positive peak hits the membrane it is moved in one direction and afterwards "sucked" in the other. Just as waves on the beach do.
Each positive peak generates a equal negative peak (in a perfect world...). Superimposition of different waves makes it hard to see this directly.

The average is not necessarily the atmospheric pressure, but is a combination of it in combination with the property of the electrical circuit of the microphone.

For your fft you have to consider positive as well as negative values. However, you may want to get rid of the average value if you do not care for very low frequencies.


What is a memory less system? As far as I have grasped the concept is that if we are given the expression

$$y[n]=(2x[n]-x^2[n])^2\,,$$

it is a memory less system because even if we give negative values of $n$, we still get the overall result in the positive sign, right?

And this expression is memory system because we get the negative answer for negative values of $n$, right?

$$y[n]=(2x[n]-x^2[n])\,.$$

A confusion may arise, for causal systems, from mistaking "having negative signal values (amplitudes)" and "depending on  negative time indices".

A strict memory-less system does depend neither on past (for the causal case) nor on future values, only the current ones, to determine the current value of the output. The output at $n$ only depends on inputs at $n$. For a strict memory-less and causal system, you could not depend on former (negative, in a way relative to the current index) time indices.

By extension, some allow memory-less systems to allow a constant delay, i.e. to hold only one value at a given lag: the output at $n$ only depends on inputs at $n-k$  of a fixed integer $k$: So, globally, memory-less system can afford the following shape, with $k$ any integer:

$$y[n] = f(x[n-k])\,.$$

For realizable causal systems, that cannot look ahead in the future, on:y $k\ge 0 $ are allowed.


How to segment phone call audio into silence/non silence? My problem is that I don't know the energy of the background noise, so I can't just threshold the energy. The processing is done in real time, and I have about 500msec to decide.
Ideally, I'd want quiet consonants considered non-silence.

There are a bunch of parameters that you can look at:

Overall energy
Short term spectrum: Speech has a fairly distinctive "pink-like" spectrum and noise (which is happening during the non-speech parts) tends to be white if it's electrically dominated or "red" (i.e. low frequency heavy) if it's acoustic background noise or microphone noise
Amplitude statistics. Most noise signals have a Gaussian distribution, speech is closer to a Laplace distribution

I think a combination of these three should give a fairly robust detection scheme.


Instantaneous velocity and displacement from acceleration signal using a proper filtering method first I need to mention I'm new to signal processing. here is the situation:
I have an acceleration time-series derived from an accelerometer

I wanted to imply a filtering method like high pass filter to denoise the signal. Using Fast Fourier Transform, I need first to decite on Cutoff frequency value.
After looking at fft(acc) I thought the cutoff frequency should be 5hz
I must note for x-axis (frequency) I just invert time (Fr=1/t) hope it's correct.

the estimated FFT's values are complex (not real numbers) which I force to work on the real part not imaginary. here is the code I used to filter acceleration signal noise by fft-highpass filter:

        % Matlab_code
   % Filtering signals by fft-HighPass
%%
acc=xlsread('accel_signals.xlsx',1,'B2:B4036');
t=xlsread('accel_signals.xlsx',1,'A2:A4036');
figure(1)
plot(t,acc,'b')
xlabel('Time(s)');ylabel('Acceleration (m/s^2)');
hold on
%%
Ts=mean(diff(t)); % Sampling rate
Fs=1/Ts;          % Sampling Frequency
Fc=5;             % Cutoff frequency = 5 hertz
fft_aac=real(fft(acc));
signal_temp=[acc,1./t];
signal=signal_temp;
for i=1:length(signal_temp)
    if signal_temp(i,2)<5
        signal(i,1)=0;
    else
        signal(i,1)=signal_temp(i,1);
    end
end
filtered_acc=real(ifft(signal(:,1)));
plot(t,filtered_acc,'r')

    %

Now the graph below is the result. the blue line is the noisy one and the red line is filtered signal. the filtered acceleration signals aren't even in the range of acceleration data but it's simply a noisy straight line.

here I listed my questions:
1- why is that happening? filtering didn't work well!
2- Do I need to use Hi pass or low pass filter by the way
3- What is the right way to choose cutoff frequency

please also help me with commenting on code and the way I approached filtering signal noise.

And after having a good filtering do I need to use a numerical integration (like trapezoidal integration method) to measure the instantaneous velocity and position?
Thanks in advance

There a a number of conceptual problems with you code.  The first has to do with filtering by zeroing out FFT bins.  This is covered in:

Why is it a bad idea to filter by zeroing out FFT bins?

The second has to do with just taking the real part of the fft.  I don't understand why you did this but this is the wrong way to do the wrong thing.

If your near term priority is to process your data.  I suggest you discard most of your code and take an entirely different approach with using the built in Matlab timeseries class and use one of the built in filter methods on the object you create. The plot methods are also easy to use and you can set units on the object you create.  This is the most direct way to a "correct" result with your data.

If you want to understand how to write your own correct code, this will take more time.  There is more than one way to filter data and choosing which takes some time to learn.  You can look at:

What resources are recommended for an introduction to DSP?

Essentially, I'm pointing you at a number of books that are several hundreds of pages of study, which is the long term answer to your problem. You will know what is correct.

If you do some Googling, you can probably find something in the middle that suits your needs but I can't recommend anything because I don't know what you know and don't know.

If you may permit the analogy, If one sees someone drowning, do you throw them a floatation device or try to teach them to swim.

Someone else here may provide a suitable solution.

It is to your credit that you have shown genuine effort.


What is difference between frequency and pitch?in nutshell and simple words? Although both are some how directly related
Such that when one increases,the other one also increases

I have read different Google results but couldn't understand .

I want to know what is exact difference between both these terms?
Please kindly explain your answer with example preferably having graphs/wave forms for demonstration of difference

Frequency is mathematically defined as the number of cycles per second. So it is a more strict word mathematically. It is represented numerically by the unit called Hertz.

$f=1/T$, where $T$ represents the one-period length of a waveform. This makes frequency quantifiable.

Pitch on the other hand, is a perceptual characteristic of a sound frequency, so it's a word used to describe how humans perceive a tone as different or not. That's why it is mostly used in the field of music more than mathematics. It is closely associated to the frequency of a sound, but not all frequencies can by described as different pitches (e.g. 440 Hz and 441 Hz would practically sound the same pitch, but a 440 Hz and 466 Hz would sound as two different pitches).

That's how musical scales were invented---to create a scale that classifies the fundamental frequencies and label them as distinct pitches (or notes). A good musical scale represents the foundation of music theory for different styles of music.


is y[n] = x[n] + n time invariant? My steps were as follows:

$\ x_2[n] = x[n-k] $

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

and

$\ y_2[n] = x_2[n] + n = x[n-k]+(n-k)$

Does this mean that it is indeed time invariant?

No; the system given by

$$ y[n] = x[n] + n $$ is time-varying, due to the added term $n$.

Your mistake is in the line :
$$\ y_2[n] = x_2[n] + n = x[n-k]+(n-k) $$ which should be instead

$$\ y_2[n] = x_2[n] + n = x[n-k]+n$$ and therefore implies that

$$y_2[n] \neq y[n-k] $$.


Largest possible downsampling factor With the intention of converting a single rate LPF FIR filter to multi rate equivalent: what is the largest possible downsampling factor? I found one source that said it was the sampling frequency divided by the frequency at the centre of the transition band, but am seeking verification.

There is not one correct answer to this question, other than the theoretical one which would be that the maximum downsampling factor $M$ is given by

$$M=\frac{f_s}{f_c}$$

where $f_c$ is the highest frequency of the (ideally band-limited) signal. Since there's no ideally band-limited signal, it depends on the application how much aliasing can be tolerated. The frequency $f_c$ must be chosen as the frequency above which the signal (after low pass filtering) has no significant energy. What 'significant' means is determined by the application.


Analyzing a signal that contains frequency content at Fs/2 doesn't seem to work unless there is a phase shift I am trying to write a basic program that samples a 4 kHz sinewave at a sampling rate of 8 kHz and takes the FFT of the signal and plots it.

From everything I have read, as long as the signal you are sampling has frequency content that is less than or equal to Fs/2 no aliasing will occur and the results will be accurate. However writing a simple example seems to be more complicated than I thought.

Using Python I wrote up a basic example:

import matplotlib.pyplot as plt
import numpy as np

Fs = 8000.0
Ts = 1/Fs
N = 8
t = np.arange(0, N*Ts, Ts)

# 4 kHz sinewave
y = np.sin(2*np.pi*4000*t)

# Bandwidth of the signal (Hz)
BW = Fs/2

# Spectral Lines (number of frequency samples)
SL = round(N/2 + 1)

# Frequency scale
f = np.linspace(0, BW, SL)

Y = np.fft.fft(y)*2/N
Y = Y[:SL]

fig, ax = plt.subplots(2, 1)

ax[0].stem(t, y, use_line_collection=True)
ax[0].set_xlabel('Time')
ax[0].set_ylabel('Amplitude')

ax[1].stem(f, abs(Y), use_line_collection=True)
ax[1].set_xlabel('Freq (Hz)')
ax[1].set_ylabel('|Y(freq)|')

plt.show()

Looking at the output you can see that this does not properly capture the frequency content at 4 kHz in the FFT plot.

But if I simply add a phase shift to the signal by changing this line

y = np.sin(2*np.pi*4000*t)

to

y = np.sin(2*np.pi*4000*t + np.pi/2)

I end up with much better results...

My question is, how come this doesn't work without a phase shift? It makes me wonder how anyone can be confident in the result of an FFT when the signal being sampled has some frequency content that is at the frequency Fs/2. How do you ensure that the signal you are analyzing is at a phase alignment that will not cause issues?

the discrete function $$x_q[n]=\sin(\pi n)$$ is always zero for all of the integers $n$.

the discrete function $$x_i[n]=\cos(\pi n)$$ is always $(-1)^n$ for integer $n$.

so this general sinusoid at Nyquist that has a phase term:

$$\begin{align}
x[n] &= A \cos(\pi n + \theta) \\
     &= A \big( \cos(\pi n) \cos(\theta) - \sin(\pi n) \sin(\theta) \big) \\
     &= \big(A\cos(\theta)\big) \cos(\pi n)  +  \big(-A\sin(\theta)\big) \sin(\pi n) \\
     &= \big(A\cos(\theta)\big) (-1)^n  +  \big(-A\sin(\theta)\big) \cdot 0 \\
     &= B (-1)^n \\
\end{align}$$

So now, how do you tell the difference between a sampled sinusoid that had amplitude $A$ and a phase angle $\theta$ and another sampled sinusoid that has amplitude $B \triangleq\big(A\cos(\theta)\big)$ and no phase shift from the cosine?


Localized Gamma Correction Question:

Is there an algorithm that does (what I would call) "localized" gamma correction?

That is, applying different gamma corrections on different parts of an image, using a masking filter or something along those lines.

I have two images that seem to be correlated by a mix of exposure and gamma corrections, but not equally on the whole image:

First image

Second image

For further details on how these images were created, please refer to this (unanswered) earlier post.

Gamma Correction is Pixel Wise operation.
Hence what you can do is estimate it per pixel and then average it per local area.


How do I find the ROC of a system if it has no poles The output of a system of discrete time $y[n]$ is corellated with the input $x[n]$ through the equation $y[n]$.

$$y[n] = \frac 13\big(x[n-1]+x[n]+x[n+1]\big)$$

It then asks me to find the system function, the ROC, the zeroes and the poles.

I found out that $$Y(z) = \frac 13\left(z^{-1}X(z)+X(z)+zX(z)\right)$$

I found out that the system function is : $$H(z) = \frac {Y(z)}{X(z)} = \frac 1{3z}(z^{2}+z+1)$$

The pole here is 0. How do I tell if this sequence is left/right/two sided to figure out the ROC?

This is an unsolved exercise given by my professor.

Since your time domain equation is of finite length, its Z-transform sum converges for (almost) any given $z$. Hence, the ROC of your system function spans the entire complex plane except for

$z=0$, which is included only if your function is anti-causal and contains only negative potencies $z^{-n}$ which evaluate as $\frac{1}{z^n} \to \infty$ for $z=0$, or
$|z| = \infty$, which will only be in the ROC for purely causal functions. Causal components $z^n$ will approach infinite values when $z$ is $\pm\infty$ (or any complex number of infinite radius), hence the absolute value.

Since your function has both causal and anti-causal parts, the sequence converges for $0 < |z| < \infty$. The region is limited by the two poles at $z=0$ and $|z|=\infty$.


How we can encode/decode sparse signals? I have question and looking for help. Suppose we have a vector of real values (lat's say 64 length resulting from factorization 8*8 block image). We got a sparse representation of that vector (let's suppose that the sparse vector is 128 length with very few active elements, about 4 coefficients). I want to send this sparse vector through a channel, so I need bit stream of 1,0. Are there any efficient encoding/decoding technique for the sparse signals ?

Since you have tagged your question compressive-sensing, I would guess that you are interested in a compressive sensing solution to this. In principle, this could work by taking a number (less than 64) of random linear combinations (measurements) of your length 64 vector and transmitting these. At the receiver, the fact that this vector has a sparse representation allows you to reconstruct it again from these measurements. However, compressive sensing does not work that well on low-dimensional vectors (64 is relatively small) and so I do not expect this scheme to work very well.

A more natural solution, since you mention that you know the sparse length 128 representation of the vector, is to use run-length encoding (Wikipedia explanation) of the sparse representation and transmit that. Since the sparse vector is almost all zeros, this should be very efficient.


Trying to implement MATLAB $\tt pwelch$ in Python I'm implementing an existing MATLAB script in Python and noticing some differences in the behaviour of MATLAB's pwelch function compared to scipy.signal.welch.

From the docs, scipy.signal.welch uses a default window size of 256, whereas MATLAB's pwelch doc says:

  By default, x is divided into the longest possible segments to obtain as close to but not exceed 8 segments with 50% overlap".

Both use Hamming windows and overlap by 50%, so I assume the differences I'm seeing are due to window sizes.

However, if I set the same window size for both, I still see different frequency resolution in the results. e.g. for data series of length 6144, sampled every 0.292969 seconds, and set window size as follows:

n = 1000
freqs, power = signal.welch(timeseries, fs=1.0/0.292969, scaling='density', nperseg=n)
print(freqs[2]) # output 0.0034

n = 1500
freqs, power = signal.welch(timeseries, fs=1.0/0.292969, scaling='density', nperseg=n)
print(freqs[2]) # output 0.00227555370371

for MATLAB:

n = 1000;
dt = 0.292969;
[psd,freq] = pwelch(timeseries,n,[],[],1/dt);
freq(2) % output 0.0033

n = 1500;
[psd,freq] = pwelch(timeseries,n,[],[],1/dt);
freq(2) % output 0.0017

For higher window size the difference in output resolution is quite stark, and results in quite different amplitudes for peaks in the spectrum. Am I missing some other parameter which needs adjusting?

This was an oversight on my part - I needed to adjust nfft to be the next power of 2 above window size.


How do you do signal averaging on a realtime data? Smoothing algorithms implemented by me

Exponential smoothing - https://www.openprocessing.org/sketch/771368
Moving Mean - https://www.openprocessing.org/sketch/540696

For example - a1+a2+a3+a4+... aN = A/N ..... is an average of all the a's which isn't same as ---> (((a/1+a2)/2 + a3)/3 +a4)/2 .... so on...

can anyone explain me how wrong I am or explain me signal averaging correctly?

@Greyfrog. Here are the descriptions of four different kinds of averaging operations:


How do I find variance from the PSD of a stochastic process? I have a time series that consists of noise and a signal, shown here windowed and Wiener filtered:

and the PSD of just the noise (used in filtering):

I want to find the variance of the noise using the PSD in order to estimate the significance of a measurement. I'm not super well-versed in signal processing, so I apologize in advance for any conceptual blunders! From digging around Wikipedia and StackExchange, I found that I can estimate the variance of a deterministic process like so:

$\sigma_{x} ^{2} = \int_{-\infty}^{\infty} S_{x}(f)df - \mu_{x}^{2}$

where $S_{x}(f)$ is the PSD and $\mu_{x}$ is the mean (side note: should there be a factor of $\frac{1}{\pi}$ or $\frac{1}{2\pi}$ in front of the integral?).

I then found this answer that confirmed that.

My question is: what modifications, if any, do I have to make to that formula to apply it to a stochastic process (the noise)? Also, the answer linked above specifies a stationary process- do I have to detrend the noise or can I assume that it is stationary?

The factor of $1/\pi$ or $1/(2\pi)$ in front depends on the convention chosen for $S_x$, i.e., one-sided vs two-sided PSD, and also expression of frequency in Hz vs rad/s. Therefore, it's a matter of using consistent conventions between the PSD and the integral. You can check that by seeing if the Wiener-Khinchin theorem is satisfied, https://en.wikipedia.org/wiki/Wiener%E2%80%93Khinchin_theorem.

The bigger issue is the evaluation of the PSD integral near $f=0$, which has to do also with your point about nonstationarity. You need to numerically evaluate the integral over the PSD, which is derived from the DFT of the original signal. The $\mu_x$ term removes the exact $f=0$ term in the DFT. However, there is a subtlety about the low frequency components: in practice, you rarely want to estimate the variance of the entire time-domain signal. I suspect that in applications, you would be more interested in the variance of the signal near the center of your window. Including the window edges will give a lower variance. By limiting the integral as

$\sigma_x(\tau) = \left(\int_{1/\tau}^{\infty} + \int_{-\infty}^{-1/\tau} \right) S_x(f) df$

you effectively detrend the data on timescales below $\tau$.

None of this depends on whether the signal is deterministic vs. stochastic - there's no way for the mathematics to tell the difference.


Why is beamforming needed in 5G? With the advent of 5G Era, base station density
increases, the cell radius decreases, and LOS propagation
gets more dominant. So, why do we need for beamforming in 5G NR, although we have smaller cells and we have a small number of devices in a cell?

One reason is that higher frequencies are envisioned. With higher frequencies, the path loss grows (cf. Friis equation). Also, the wavelength is reduced and thus, $\lambda/2$ radiators start becoming quite small. The power you can radiate from a small aperture is limited as well so overall our power budget suffers. The only way out is directivity (cf. Friis equation again): we need to make sure our energy is directed towards the users of interest. This is possible by employing large antenna arrays (i.e., arrays comprised of a large number of elements). If each element radiates a certain power, $N$ elements can add up coherently to the $N$-fold power, provided we send it in the right direction. This is what beamforming does.
Another reason is that beamforming allows us to spatially separate users so that communication links can share the same time and frequency resources and be separated only spatially. In theory, an $N$-element antenna array can serve $N$ data streams at the same time and frequency, either to one user (that receives $N$ spatially multiplexed data streams, as it is common in WiFi) or to $N$ users (that each receive one). This increases network capacity even more.
Both factors play well together: higher frequencies and smaller wavelengths lead to a more rapid decorrelation of the channel in space so that even relatively closely spaced users can still be spatially multiplexed. Think of very narrow "pencil" beams, each pointed to a different receiver (if we are serving their lines of sight) or path (if we are using multipath).
While this sounds good, it also poses significant challenges. If we only slightly depoint the beam, the signal may be lost entirely. For moving users, this can be quite a challenge. Initial acquisition becomes harder too because signals transmitted/received omnidirectionally might be too faint to detect.


OFDM time vs. frequency domain channel estimation/equalization In OFDM, the majority of Equalizers are used in frequency domain.  I mean the signal is transmitted in time domain (after performing iFFT), then to estimate the channel we should perform FFt to use any type of equalization such that LS or MMSE.

My question, what's about if we perform the equalization in time domain? it means at receiver side, we perform any type of estimation/equalization before performing fft (estimate and equalize in time domain)?   What's the advantages of using frequency domain OR  using time domain estimation and equalization?

Thank you

Traditionally, OFDM became popular in WiFi and LTE because the channel model consisted of multi-path. That is, the radio signal transmitted in 1-6GHz frequencies bounced from various obstacles (walls, trees, cars, humans) at the receiver. Of course this is time varying because obstacle position or transmitter/receiver position also changes. But to simplify the computation, if we consider a multi-path channel at time $t$, its baseband model for $L$ tap multipath is

$$
h[n] = \sum_{k=0}^{L-1}a_k\delta[n-k]
$$
where $a_k = r_ke^{j\phi_k}$ is the complex tap for path $k$. If you plot the FFT of this channel, you will find it is not a flat but rather having different gains at different frequencies. To equalize this in time domain, you would use adaptive filtering techniques. The problem is complexity and time taken for adaptive filtering to estimate the channel. By the time adaptive filter error tapers down, channel would have changed (due to movement of transmitter/receiver/obstacle). So for each OFDM burst, there is a need to immediately equalize to demodulate the packet and ultimately show the data on the device! This is where advantage of OFDM comes into picture. Even though the complexity of OFDM is higher (overhead of preamble, pilot symbols, need to maintain orthogonality), the equalization simplifies to a single tap equalizer in frequency domain as the whole frequency selective channel is split into small flat-fading channel due to IFFT -> Cyclic Prefix -> FFT method. You make sure the circular convolution happens by inserting cyclic prefix so in frequency domain it turns into a point-by-point equalization $Y[k]/H_{eq}[k]$ where $Y[k]$ is the received symbol at subcarrier $k$, and $H_{eq}[k]$ is the equalizer tap you computed using pilots.

EDIT: I understand the channel model given above is very simplified view of multi-path model. The simplification is explained in the famous book by Prof. David Tse and Prof. Pramod V (https://web.stanford.edu/~dntse/Chapters_PDF/Fundamentals_Wireless_Communication_chapter2.pdf)


Phase correlation vs. normalized cross-correlation I asked this over at Mathematics Stack Exchange, but since this sort of lies on the border of the questions normally asked over there and the questions you see over here I'll ask it here as well. (As of now, there's been no activity on my question over there.)

In 2-dimensional discrete signal analysis (specifically image processing), a definition I found for the normalized cross-correlation between two images, both of size $M\times N$ $g_1(x, y)$ and $g_2(x, y)$ is:

$$r_1 = (g_1 \star g_2)(x, y)_{\rm Normalized} = \frac{\sum\limits_{m=0}^{M-1} \sum\limits_{n=0}^{N-1}\left[g_1(m, n) - \overline{g_1}\right]\left[g_2(x + m, y + n) - \overline{g_2}\right]}{\sqrt{\sum\limits_{m=0}^{M-1} \sum\limits_{n=0}^{N-1} \left[g_1(m, n) - \overline{g_1}\right]^2 \left[g_2(x + m, y + n) - \overline{g_2}\right]^2}}$$

This is supposedly performed by the function normxcorr2 in the Signal Processing Toolbox in MATLAB, although using the Fast Normalized Cross-correlation algorithm by J. P. Lewis. Compared to the Phase Correlation method (with normalized cross-power spectrum) suggested by Kuglin and Hines:

\begin{align}
G_1(u, v) &= \mathcal{F}\left\{g_1(x, y)\right\}\\
G_2(u, v) &= \mathcal{F}\left\{g_2(x, y)\right\}\\
r_2 &= \mathcal{F}^{-1}\left\{\frac{G_1^\ast (u, v) G_2(u, v)}{\lvert G_1^\ast (u, v) G_2(u, v) \rvert}\right\}
\end{align}

Without the element-wise normalization before the inverse Fourier transform, $r_2$ is the same as non-normalized cross-correlation, with the exception that the Fourier transform assumes that the signal repeats in the spatial domain. It is clear that $r_1 \neq r_2$ by looking at the resulting correlation "images" in MATLAB, which I expected, but $r_2$ almost seems discontinuous from the images I am testing the method on, while $r_1$ always turns out much smoother (the resulting $r_2$ images are always "spotty", $r_1$ are not). Why is this? I assume it has something to do with the element-wise normalization in the Fourier domain, which I suppose is basically like cropping the signal in a way, but I don't know how to conclude this from any known properties of the (inverse) Fourier transform, if this is even a correct assumption.

Here are examples of the differences ($r_2$ has been amplified to clearer show the difference in characteristics of the images). They are performed on the same 2 images.

$r_1$: Normalized cross-correlation

$r_2$: Phase correlation

Summarizing my question: Why is $r_2$ so "spotty", while $r_1$ is not?

I came across this conundrum in a 1 dimensional case, which is how I will present it.  Consider two signals that you want to correlate.  Signal 1 (figure panel a) is a damped sine wave and signal 2 (figure panel b) is two instances of signal 1 but at slightly different amplitudes.

Now consider using a normalized cross correlation as you defined in your question (figure panel c).  The result is a smoothly oscilating function that has a peak amplitude when the peak from signal 1 and the peak from signal 2 are aligned.  Conversely the normalized cross correlation function has troughs when the peak from signal 1 lines up with the troughs from signal 2.

Then consider using a phase correlation as you defined in your question (figure panel d).  A phase correlation involves a division of the absolute value (or amplitude) of the two signals.  In a 1D case that is a vector and in a 2d case, as your question is, that is a matrix.  In either case the values in that divisor can be very close to zero and when you divide by something close to zero, you get spikes or discontinuities.

There are several strategies around this and include low pass filtering the input signals or the phase correlated function.  This page on stack overflow might be helpful: https://stackoverflow.com/questions/30630632/performing-a-phase-correlation-with-fft-in-r


DC Coefficient of JPEG Block Suppose I take a grayscale image and break it into 8x8 pixel blocks, then take the DCT of those blocks.  What is the DC coefficient of the block in relation to the whole image?

To simplify I'm thinking of this in 1D.  I have a row of 64 pixels and I break it into length 8 vectors.  I pick one of these 8-vectors and take its DCT. The average of the 8-vector is not the DC term of the entire 64-row. Now suppose I took the DCT of the entire 64-row. What coefficient in the 64-pixel-DCT would the DC term of the 8-pixel-DCT be projected onto?

EDIT:
Thank you for the replies so far.  The context of this question is I'm trying to assign a cyles/pixel value to the DC coefficient of the 8-pixel sub-block, but in the context of the whole 64-block.  Its clear to me now that the DC coefficient of the 64-block is the average of the DC coefficients of each of the 8-blocks. But while the 8-block is a sample of the 64-block it doesn't make sense to me to say that the 8-block DC coefficient is actually 0 cycles/pixel. For example, one could represent nonzero frequencies using only these 8-block DC coefficients.

The DC coefficient, of the 2D-DCT (discrete cosine transform) of an 8 x 8 image block, represents the average value of the samples within the 8 x 8 block.

For your question, if you compute 8 different 1D-DCTs of 1 x 8 rows or 8 x 1 columns of this 8 x 8 block, then the DC coefficient of 2D-DCT of this 8 x 8 block will be the average of individual DC coefficients of 8 x 1 (or 1 x 8) 1D-DCTs.


How to convolve an arbitrary signal with a causal decaying exponential? I need to simplify the following convolution

$$x(t)\star [e^{-2t} u(t)]$$

where $u(t)$ is the unit step function. I'm very confused with this. Using the definition of convolution of continuous-time signals, I obtained

$$\int\limits_{-\infty}^t {x(\tau)e^{-2(t-\tau)}u(t-\tau)}{d\tau}
$$

However, I have no idea what to do at this step. I'm stuck. Is there any property / theorem / trick that can help me simplify this?

Indeed you have reached what can be reached, may be the following additional line can be obtained by moving the $t$ function $e^{-2t}$ out of the integral and replacing the $u(t-\tau)$ by $1$ as:

$$
y(t) = e^{-2t} \int\limits_{-\infty}^t {x(\tau)e^{2\tau}{}d\tau}
$$

In the general case you cannot simplify it any further...


Why is privacy scrambling a standalone block in a communications transmitter? In some telecommunication systems, scrambling of transmitted data is applied to provide some level of security. For example, the communication system might look like:

Data -> Scrambling -> Compression -> Error correction encoding

Could one possibly match the scrambling with the error-correcting code to obtain some benefit? Why is scrambling typically separated from the other blocks of the transmitter?

There are two possible scramblers, one meant for security and
privacy of the information, and the other for making sure that
the sequence of symbols transmitted over the channel has
properties that are desirable for digital transmission purposes.
Either or both might be present in a system, or both could be
absent too. Furthermore, the scramblers might well be
non-interchangeable since the scrambling needed for privacy
need not be the same as that required for transmission purposes.

The question then is, where should these scramblers be placed
in the transmitter.  The channel scrambler obviously should be
between the error-correction encoder (whose output is to be
transmitted) and the modulator that creates the actual signals
that are broadcast over an antenna or transmitted down a cable
or twisted pair of wires (or recorded on a magnetic medium or
CD or DVD or what have you).  The question of where the
privacy scrambler should be placed is more complicated.
@JasonR has contended in comments that putting it before
the data compression unit makes the input to the compressor
less compressible.  This is certainly something to think
about, but I am not completely sure about that it is always
a major concern
because the scrambler for security and privacy might be doing
other things than breaking up long sequences of repeated
symbols.  In any case, the location of the privacy scrambler
before the data compression unit may also be dictated by
privacy and security concerns that may have nothing to do with
efficient communication.  If I were in charge of privacy and
security, I would make darn sure that the chain looked like

data --> ENCRYPTION -- > scrambling -->

in my domain before the dsp and comm people even get to see
the stream of symbols.  Let them worry about
how to transmit the stream I give them as efficiently as
possible; my job is security and so I make sure that things are
nailed down tight when they leave my shop!


Estimating noise level without 'quiet' segment I have a signal which consists of an onset and a release period.

The onset, sounding like a thud, excites frequencies across the entire frequency range (0-22 kHz). It lasts for about 0.1 s.

The release, sounding like a ring, excites 2-5 specific narrow frequency bands. Their amplitudes decay exponentially, albeit at slightly different rates.

I want to find the end of these narrow peak tracks. The method I am using now tracks the max amplitude in the narrow frequency bands and takes a windowed derivative. When this derivative passes through a threshold value, the end of the peak track is marked. This doesn't seem to be working consistently for particularly noisy signals.

My second idea was to track the amplitude of the peak tracks and compare them to the noise level. As soon as they reach the noise level, mark it as the end.

How do I estimate the noise level if there is no non-signal segment? The segment  before the strike, is often very short or has some ramp in amplitude making it unsuitable as an average.

Edit: Here's the spectrogram:

If you want to determine if the 2-5 narrowband frequencies are present use a matched filter with a threshold to determine their presence. In matlab the function is xcorr.
If you want to determine what the noise level is when only signal is present, try filtering out the signal. Create a filter for each frequency then convolve the filters together to get a final filter.


Does spectral leakage include negative frequencies? When choosing a window function, window duration, and/or transmission frequency (assuming sampling rate satisfies Nyquist), one may want to understand what sort of spectral leakage would occur at a frequency of interest.

It is known that a finite-duration window corresponds to a non-finite bandwidth frequency response (e.g. rect <-> sinc), and it is also known that a multiplication in the time domain corresponds to a convolution in the frequency domain.

Consider the simple case of a non-windowed, constant-frequency sinusoid in the time domain, which corresponds to a frequency response of two delta spikes centered around 0 Hz.

Applying a window function would convolve the frequency response of the window function (e.g. a sinc function) over the delta spikes.

1. Does this convolution of the window function extend to the negative frequencies when calculating the spectral leakage in the positive frequency components (and vice versa)?

I would say yes based on the above (time domain multiplication <-> frequency domain convolution), and the following two images (source) which I annotated in red.  But it leads me to question 2, which I find a bit concerning.

2. If spectral leakage does extend to opposite-sign frequencies, doesn't that imply window functions without a zero crossing at 2x of a frequency of interest would result in constructive and/or destructive interference ("spectral leakage") in the frequency response of the windowed duration at that frequency?  I.e., frequencies could interfere with themselves?

Here is an image (source) showing that even for a given window duration, the occurrence of such "self-interference" would depend on the choice of window function:

Yes spectral leakage includes both negative and positive frequencies; i.e, the convolution is from $-\pi$ to $\pi$.

However its effect depends on what exists at the negative frequency region, when you want to compute the resulting leakage on a given positive frequency.

On a typical application where you have a real valued sine wave of frequency $f_0$ as the windowed signal, then there is a dirac impulse at the negative frequency $-f_0$ as well as the positive frequency $f_0$. Then both impulses shift the spectrum of the window Fourier transform to their own locations and create leakage effects on any other frequency (unless they have zero value at that particular frequency). However the leakage is expected to be smaller from a negative frequency as it's side lobes are already decaying more than the positive frequency contributions. (This can be the opposite when both impulses are close to $\pi$ or close to $0$ radians.)

But if you have a complex signal, such as an analytic signal which has zero negative spectral content, then effectively the negative frequencies will not create leakage.


Unsure if my unit step function is correct in relation to a diagram Here is my function:

$$3(t)(u(t)-1)) + (u(t-1)+3) - (u(t-3)+3) - 2(t)(u(t-3)+3) + u(t-4) -2(t)(u(t-6))$$

Is my function correct? If not, could you please tell me where I am wrong.

Thank you in advance.

If the diagram was given to you than your function is incorrect. First of all you don't need the small brackets for the time. If you look at the diagram you can see that you can disassemble it into smaller parts:
Part I: from 0 to 1 you have a curve which starts at (0,0) and ends at (1,3). If you plug in these $(x,y)$ values into the equation $y=kx+n$ you will get $n=0, k=3$ which means that the curve function is equal to $y=3x \Rightarrow f(t)=3t$. Since the interval is from 0 to 1 we get $u(t)-u(t-1)$ and multiplied by the curve function $3t[u(t)-u(t-1)]$.
Part II: from 1 to 3 you have a rectangular pulse with an amplitude of 3 which is $3[u(t-1)-u(t-3)]$.
Part III: from 3 to 4 you have again a curve whose function you retrieve the same way as in part I. You will get $(9-2t)[u(t-3)-u(t-4)]$.
Part IV: from 4 to 6 you have again a rectangular pulse but with an amplitude of 1 so you don't have to multiply the difference of the step functions because it's 1 by default. You will get this $u(t-4)-u(t-6)$.
Part V: Assuming that your signal equals to 0 for $t>7$ you again have a curve whose function you retrieve as in Part I and you get $(2t-11)[u(t-6)-u(t-7)]$.  Summing up all obtained parts your function would be:
$f(t)=3t[u(t)-u(t-1)]+3[u(t-1)-u(t-3)]+(9-2t)[u(t-3)-u(t-4)]+u(t-4)-u(t-6)+(2t-11)[u(t-6)-u(t-7)]$  If you use MATLAB here's the code and the plot:

t=0:0.01:8;
x = 3*t.*(stepfun(t,0)-stepfun(t,1))+3*(stepfun(t,1)-stepfun(t,3))+(9-2*t).*(stepfun(t,3)-stepfun(t,4))+stepfun(t,4)-stepfun(t,6)+(2*t-11).*(stepfun(t,6)-stepfun(t,7));
plot(t,x,'r');
axis([0 8 0 4])
grid on;


Finding period of a square wave with varying sampling frequency I have a square wave (0-1.8V) with a varying sampling frequency (from a circuit simulator). It is also not a perfect square wave (the high and low signal could be very close to but not precisely zero or 1 Volt) and the fall time is slightly slower than the rise time.

Another factor is that I will retake the data and the period (and possibly duty cycle as well) will change when I alter another parameter, so I am not working with a fixed number of cycles per window.

I tried autocorrelation but it seems that it gave me double or so the actual period. And fft seems to need a fixed sampling rate.  What would be the best solution to finding the period of this square wave?

Pick a nice center value, say 1.0, scan across looking for crossings from down to up (or the reverse).  Mark the spot of your first crossing, count up some fixed number of crossings, call the count N.  Take your system timer value at your first crossing and your last crossing.  Divide the time interval by N and you'll have a very good estimate of the period.  How uneven the crossings are in the interval is immaterial.  If you can assure the two crossing that you are timing happen in "the same state", you will lose one more level of uncertainty.

You can't do it much simpler than that.  Increased accuracy will come from interpolating the crossing point rather than taking the closest sample's time.

----------------------

Follow up

You'll have to figure the Interpolate routine out yourself.  It's not that hard.

The logic is very straightforward and significantly less computationly heavy than an auto-correlation technique.  More accurate in this case too.

Pseudo-code

    FirstIndex = -1
    N = 0

    For i = 0 to PointCount
      If signal[i] <= 1.0 and signal[i+1] > 1.0 Then
         If FirstIndex < 0 Then
            T1A = System.Time(i)
            T1B = System.Time(i+1)
            FirstIndex = i
         Else
            TNA = System.Time(i)
            TNB = System.Time(i+1)
            LastIndex = i
            N += 1
         End If
      End If
    Next

    T1 = Interpolate( T1A, T1B, Signal[FirstIndex], Signal[FirstIndex+1] )
    TN = Interpolate( TNA, TNB, Signal[LastIndex], Signal[LastIndex+1] )

    Period = ( TN - T1 ) / N

You may not need to interpolate, (TNA - T1A)/N may be accurate enough.


What is the difference between Remez exchange and Parks-McClellan filter design? A question recently came up regarding Parks-McClellan and some of the comments point out that the wikipedia article on Parks-McClellan states,

  ...the Parks-McClellan algorithm is a variation of the Remez algorithm or Remez exchange algorithm.

Is there a difference in performance?  Is the remez function in octave equivalent to the Matlab function firpm?

The Remez exchange algorithm is a generic iterative procedure to approximate any function optimally in the $L^\infty$ sense (i.e., give the best worst-case approximation or in other words, minimize the maximum error or minmax).  The Parks-McClellan algorithm (PM) is a variation of the Remez exchange algorithm, applied specifically for FIR filters. From the wiki article you cited:

  Thomas [Parks] drove from Houston to Princeton in order to attend a conference. At the conference, he heard Ed Hofstetter's presentation about a new FIR filter design algorithm (Maximal Ripple algorithm). He brought the paper by Hofstetter, Oppenheim, and Siegel, back to Houston, thinking about the possibility of using the Chebyshev approximation theory to design FIR filters. He heard that the method implemented in Hofstetter's algorithm was similar to the Remez exchange algorithm and decided to pursue the path of using the Remez exchange algorithm.

Without going into too much detail, the primary difference between the two algorithms is that the Remez exchange algorithm (RE) gives you conditions to design the optimal filter (specifically, see #3 here: the errors must be of equal weighted magnitude and alternating in sign).  RE implements an iterative procedure to calculate polynomial coefficients (which can be mapped to FIR filter coefficients) that satisfy the above criterion, which is the "Alternation Theorem".  The "E" in RE is the part of the procedure where the maxima in the error, that are used in the iterative procedure, are replaced by the new maxima, which are closer to the optimal values.  PM uses Tchebyshev polynomials to convert polynomial coefficients to coefficients governing a series of cosine functions which are directly translated to symmetrical FIR coefficients.

Coming to your question about the remez function in Octave and the firpm function in MATLAB, I believe they're the same. MATLAB used to have a remez, which was phased out in favour of firpm. Octave probably still sticks to the former. Typing help remez in MATLAB R2011b gives the following:

  REMEZ Parks-McClellan optimal equiripple FIR filter design.

REMEZ is obsolete.  REMEZ still works but may be removed in the future.
Use FIRPM instead.

See also FIRPM.


how do I derive the system equation for a simple delay with feedback? I am a software engineer, and just learning digital signal processing formally, though I've hacked around before a fair amount.

I'm implementing a delay audio VST and I'm trying to wrap my head around how one creates equations of the output given a system diagram. I'm not used to circuit diagrams (they're weirdly parallel to me), so I'm trying to figure out how one goes from this system diagram:

to an equation. As best I can tell, the equation for this is:

y[t] = x[t] + y[t - T] * gain

where x[t] is the input on the far left, y[t] is the output on the far right, and the delay time in samples is T.

Is this correct?

I don't fully see how the system diagram implies the second term is y[t - T] rather than x[t - T], but intuitively, I know it must be, because an audio delay unit can't really incorporate feedback if the original signal delay only happens once (feedback must enable the ever-decaying ringing with time). Otherwise we'd get just a shifted version of the original, added to the original. This would be boring.

In such cases it helps to define an additional signal $w[n]$ at the input of the delay. Now you can write down two equations describing the system:

$$\begin{align}w[n]&=x[n]+g\cdot w[n-N]\\y[n]&=x[n]+w[n-N]\end{align}\tag{1}$$

where $N$ denotes the delay (in samples) and $g$ is the feedback gain. From the second equation you can express $w[n-N]$ in terms of $x[n]$ and $y[n]$:

$$w[n-N]=y[n]-x[n]\tag{2}$$

Plugging this into the first equation gives

$$y[n+N]-x[n+N]=x[n]+g\cdot \left(y[n]-x[n]\right)\tag{3}$$

Rearranging and subtracting $N$ from all indices results in

$$y[n]=x[n]+(1-g)x[n-N]+g\cdot y[n-N]\tag{4}$$

Eq. $(4)$ is the (single) difference equation describing the system. However, it's more efficient to implement the system using the two equations given in $(1)$.


Ideal sampling using sinc funcion Let $ x(t) $ be a bandwidth limited signal such as $ \forall |\omega|>\frac{\pi}{T} : X^F(\omega)=0 $ while $ X^F(\omega)$ is the signal's Fourier Transform.

Let us denote $y[n]=\int_{-\infty}^{\infty}x(\tau)sinc(\frac{\tau-nT}{T}) d\tau$

I noticed that my professor said that because the signal is bandwidth limited and the width of the $sinc$ function is $T$ and it's orthogonal property, and with regard to ideal sampling and reconstruction we can say that $y[n]=T\cdot x[nT]$.

I understand that we can treat $sinc$ as a delta function somehow but I'm can't find the right mathematical explenation which had lead to this result.

I'd happy for some insights.

Your professor is using sloppy notation. When one writes the sample value $x[\cdot]$, it is generally accepted that the argument is an integer (which $nT$ might not be) and we relate this to the continuous-time signal $x(t)$ via $x[n] = x(nT)$. So, what your professor should have said is that $y[n]$ has value $T\cdot x[n] = T\cdot x(nT)$, and not that $y[n] = T\cdot x[nT]$ (even if $T$ happens to be an integer).  An outline of a proof that $y[n] = T\cdot x[n] = T\cdot x(nT)$ is given below.

If $x(t)$ is a continuous-time signal bandlimited to $\left(-\frac{1}{2T}, \frac{1}{2T}\right)$, that is, its Fourier transform $X(f)$ has the property that $X(f) = 0$ for all $|f| \geq \frac{1}{2T}$,then it can be reconstructed from its samples spaced $T$ seconds apart. So, let's consider the signal $x(t)\operatorname{sinc}\left(\frac{t}{T}\right)$ where
$$\operatorname{sinc}(t) = \begin{cases}\frac{\sin(\pi t)}{\pi t},& t \neq 0,\\1, & t = 0.\end{cases}$$
The Fourier transform of $\operatorname{sinc}(t)$ is $\operatorname{rect}(f)$ and so the Fourier transform of $\operatorname{sinc}\left(\frac{t}{T}\right)$ is $T\operatorname{rect}(fT)$, and the Fourier transform of $x(t)\operatorname{sinc}\left(\frac{t}{T}\right)$ is the convolution  $X(f)\circledast T\operatorname{rect}(fT)$ of the Fourier transforms $X(f)$ and $T\operatorname{rect}(fT)$ of $x(t)$ and $\operatorname{sinc}\left(\frac{t}{T}\right)$ respectively.  Now, in principle, the convolution of a function and a rectangular pulse is easy to compute but I refuse to write out the answer because it is not needed: what we need to find is the value of $\int_{-\infty}^\infty x(t)\operatorname{sinc}\left(\frac{t}{T}\right) \mathrm dt$ which, if you think about it a bit, is actually the value of the Fourier transform of $x(t)\operatorname{sinc}\left(\frac{t}{T}\right)$ evaluated at $f=0$. That is,
\begin{align}
\int_{-\infty}^\infty x(t)\operatorname{sinc}\left(\frac{t}{T}\right) \mathrm dt &= X(f)\circledast T\operatorname{rect}(fT)\bigg |_{f=0}\\
&= T\int_{-\infty}^\infty X(\lambda)\operatorname{rect}((0-\lambda) T) \, \mathrm d\lambda\\
&= T\int_{-\infty}^\infty X(\lambda)\operatorname{rect}(\lambda T) \, \mathrm d\lambda & \scriptstyle\text{rect is an even function}\\
&= T\int_{-\frac{1}{2T}}^{\frac{1}{2T}}X(\lambda) \, \mathrm d\lambda &\scriptstyle\text{rect has finite support}\\
&= T\int_{-\infty}^\infty X(\lambda) \, \mathrm d\lambda &\scriptstyle{X(\lambda)=0 ~\text{for}~|\lambda| > \frac{1}{2T}}\\
&= T
\cdot x(0).
\end{align}
Thus, $y[0] \overset{\text{def}}{=} \int_{-\infty}^\infty x(t)\operatorname{sinc}\left(\frac{t}{T}\right) \mathrm dt$ equals $T\cdot x(0)$. I leave it as an exercise for the reader to work out the Fourier transform of $\operatorname{sinc}\left(\frac{t-nT}{T}\right)$ and show that $y[n] \overset{\text{def}}{=} \int_{-\infty}^\infty x(t)\operatorname{sinc}\left(\frac{t-nT}{T}\right) \mathrm dt$ works out to be $T\cdot x(nT)$ and so if we define the samples of $x(t)$ as $x[n] = x(nT)$, then we have that $y[n] = T\cdot x[n]$ for all $n$.


Why is my laser beam intensity profile not Gaussian? I am kind of new to laser physics and signal processing, I hope I am in the right place to ask this question. I am also sorry for my not-great English.

I have been asked to study the transverse intensity profile of a laser beam produced by a $\rm He$-$\rm Ne$ source using a commercial webcam. I am placing the webcam sensor (without any lens) in front of the source and taking a snapshot. Then I am assigning each pixel an intensity value that is simply the recorded value in the red channel ($0-255$). If I choose a row of pixels and plot such value as a function of the pixel, I should expect a Gaussian profile. However it is not what I am observing.

In red you can see a Gaussian fit of the data. The actual data is shown in blue. I cannot really call this a good Gaussian shape, although it is some kind of a bell.

What could cause such a shape?
Could it be due to the laser source which could not be perfectly single-mode?
Or could it be due to the sensor?

It is interesting to note that I have also produced an image after making the beam pass through a single mode optical fiber, obtaining similar results; this makes me think that the problem lies in the sensor. However, the only modification the sensor causes and that I am aware of is Gamma correction, and as far as I know, that means that the red (or any channel) value is elevated to a power which is usually $<1$, but that should not change the actual shape. A Gaussian function elevated to a power is still a Gaussian function, although its parameters are different. Can anyone help?

What could cause such a shape?
  Could it be due to the laser source which could not be perfectly single-mode?
  Or could it be due to the sensor?

There is nothing particularly wrong with the "shape", but there are a few things you can do on the sensor and data processing side, to improve the extraction of an accurate profile.

Your "biggest" problem is going to be the saturation of your web cam's sensor because of the intensity of the beam. This is evident as a flattening of the Gaussian curve's peak. There is a bit of that in this example, but if the $Y$ axis is not scaled, then values in the vicinity of 50 don't look like being saturated....However...

Please make sure that "Automatic Exposure Control" on the web cam is off and then manually set the exposure for minimum distortion (it will probably have to be set to its lowest setting). This will also affect the slopes of the curve as exposure control will strive to maintain contrast across the image.

In addition, you might still need some "analogue" attenuation too. To achieve this, use one or more neutral density filters in front of the camera.

Finally, as it is evident from this profile, it is composed of integer quantities and you are trying to fit a continuous distribution on it. There are two things that can be done here towards improvement: The first, is to shoot a video and then create an average image of the profile over a few hundred frames. This will give you a much smoother real ($\mathbb{R}$) number profile which will vastly improve the fit. The second, is to still shoot a video but produce the beam profile as a boxplot. This will give you the opportunity to observe exactly what the profile "does" over a few hundred (or more) frames. Very briefly, a boxplot would present in one view, information about the distribution of the brightness values at each pixel of the profile over time. You can then fit a Gaussian over the means (or medians) of the profile.

Hope this helps.


Why does speeding up audio files with VLC Player not match with doing the same thing with RPG Maker or Audacity? When I speed up 150% an audio file with VLC Player, why doesn't it match the higher-pitched version I get when I speed the same file 150% with RPG Maker VX Ace's Sound Test?

When I speed up with VLC Player, I don't get a pitch increase

When I speed up with RPG Maker or Audacity, I do...

How would their methods of speed increase differ, such that that happens?

Also, I'm on a Windows, if that matters at all...

i don't have any of those products, but evidently the VLC Player does what we call "time-compression" when you speed up.  and likely "time-expansion" when you slow down.  together we call those operations "time-scaling".  time-scaling and sample-rate-conversion, together, make a pitch-shifter.  so a time-scaler can change tempo without changing pitch, a sample-rate converter can change tempe and pitch together, and a pitch shifter can can pitch without changing tempo.

the other products, RPG and Audacity, likely just use the same interpolation we use for sample rate conversion to do the speeding up.   essentially the same as speeding up playback of analog tape.


how to calculate the dc component of an energy signal? Does the energy signal or finite energy signal have 0 dc component always?How to confirm this from the frequency spectrum?

An energy signal $f(t)$ satisfies

$$\int_{-\infty}^{\infty}|f(t)|^2dt<\infty\tag{1}$$

i.e. it has finite energy. If $f(t)$ had a non-zero DC component $f_0$ defined by

$$f_0=\lim_{T\rightarrow\infty}\frac{1}{T}\int_{-T/2}^{T/2}f(t)dt\tag{2}$$

the integral in (1) wouldn't converge. So an energy signal must have a zero DC component.

Note that there is some confusion about the term "DC component". Some people would call $F(0)$, i.e. the value of the Fourier transform of $f(t)$ at frequency $\omega=0$, the "DC component" of $f(t)$. However, note that $F(0)$ given by

$$F(0)=\int_{-\infty}^{\infty}f(t)dt\tag{3}$$

can be non-zero (but finite), even if (1) is satisfied. One example of such a signal satisfying (1) but with $F(0)\neq 0$ is a sinc pulse.

Also have a look a this answer to a related question.


Can multi-rate operations of decimation and interpolation can be used to implement a rate change by a factor an irrational number? I know that common rates are integer values like 2 or 3, but can multi-rate operations of decimation and interpolation can be used to implement a rate change by a factor of an irrational number? Common sense is telling me no, because irrational numbers can't be made from any set of integers by definition, so no combination of interpolation and decimation filters can generate the irrational number. There is also the fact that the MATLAB command resample.m uses the poly phase structure for implementing a sampling rate change by a rational factor. All signs are pointing to irrational numbers can't be done. Am I wrong in this thought process?

You're exactly right, with rational operations, you can't achieve an irrational ratio.
The only way out there is using some interpolative number, i.e. an arbitrary resampler.
However, in signal processing, you basically always deal with system that only run for a finite time and within that have a boundable tolerance against jitter or frequency error - so one common solution to that dilemma is, e.g. in receivers, to do the resampling with a slightly wrong, but rational ratio, and then let the symbol synchronizer, which often does some form of arbitrary resampling internally, anyways, take care of the rest.


Processing Wifi Packets Im trying to use GNUradio 3.8 with an sdr to pick up some wifi packets 2.4 and 5 GHz but im seeing a ton of dropped samples from my sdr. Im wondering if my dropped samples have to do with the sample rate. Ive seen online sdrs have a mixer that downconvert the frequency band to baseband and then its sampled. But i believe one cannot configure arbitrary sample rates, as the receiver expects a given number of samples per symbol. Is this correct. If im trying to sample 2.4Ghz and 5Ghz a/b/g/n how do i figure out my sample rate/dropped samples issue. Thanks

At the very least, you'll need a sampling rate that's high enough to capture the full bandwidth of your Wifi transmission. In a typical WiFi receiver, that puts the sample rate at exactly the OFDM width: 5, 10, 20 or 40 MHz, depending on your Wifi.
Your dropped samples simply happen because your computer isn't fast enough to keep up with that many samples per second. There's only two solutions to that:

Get a faster computer
Make your computer do less per sample.

Also note that the sustained write rates of permanent storage devices can be lower than your sample data rate, so even just writing to disk can be a fatal bottleneck.
Sadly, we don't know the processing you're currently doing, but maybe have a look at gr-ieee802-11, which is optimal enough that a modern midrange laptop can receive 20 MHz wide Wifi with it, continuously.


How to take into account phase differences when doing Fourier Transform? I apologize if the title seems very vague, I don't really know how to put it a different way.
From what I understand, the Fourier transform is essentially just doing a dot product between multiple sine waves of different frequencies with the original signal. When the "chosen" sine wave has the same frequency as a component in the signal, it will be in phase and produce peaks. If the "chosen" wave isn't in the original signal, then it will either be cancelled by the original signal, or produce smaller peaks.
Let's say that there's a 100Hz sine wave component in a signal with a phase of $\frac{\pi}{2}$. If I were to take the dot product with a pure sine wave of 100Hz with no phase, it would appear as if two signals are present, one at 95Hz and one at 105Hz:

Is there a way to negate the effects of this? I was thinking of taking the maximum of pure cosine and sine wave in order to take into account phase differences, but I get much more noise around the frequency component:

Note: I'm not using any windowing function (Hanning, etc.) I can post my Python code if that will help.
EDIT: I tried using the formula as per Marcus' comment, and I'm getting weird results:

def do_fft(y, rate, f):
    N = len(y)
    summ = 0

    for n in range(N):
        summ += y[n] * np.exp((0 - 1j) * 2 * np.pi * n * (f/rate))
    return float(summ.real)
    #return float(abs(summ))

def selective_fft(y, rate, frequencies):
    ar = []
    for f in frequencies:
        ar.append(do_fft(y, rate, f))
    return ar

"frequencies" is an array of frequencies that I want to "test" for.
EDIT: Here is how I'm doing the pseudo-fourier transform (originally):
t = np.linspace(0, 0.1, 4410)

ar = np.arange(0, 500, 1).reshape((-1,1)) #The frequencies to test

m = np.tile(ar, (1, t.shape[0])) #create matrix for frequencies
time_matrix = np.tile(t, (ar.shape[0], 1)) #create matrix for t vector

a = 100 #frequency components
b = 200
c = 300

y = np.sin(2 * np.pi * a * t + np.pi/2) + np.sin(2 * np.pi * b * t) + np.sin(2 * np.pi * c * t)

yc = y.copy() #copy signal

y = y.reshape((-1,1)).T #Add column and transpose to get one row and 4410 columns

y_stacked = np.tile(y, (ar.shape[0], 1)) #stack the function based on how many frequencies to test (500 from before)

sub = np.sin(2 * np.pi * np.multiply(time_matrix, m))
sub2 = np.cos(2 * np.pi * np.multiply(time_matrix, m)) #use cosine

dotted= np.abs(np.sum(np.multiply(y_stacked, sub), axis = 1).reshape((-1,1))) #basically the dot product
dotted2 = np.abs(np.sum(np.multiply(y_stacked, sub2), axis = 1).reshape((-1,1)))

This is what it looks to me what you're doing: you're preparing a signal with mixed frequencies, then you're making a sin/cos matrix of multiple sines which are to be multiplied with the test signal, and then you're using the absolute value of each sine and cosine multiplications and then average them. You then consider the variation of these averages the "FFT". If I didn't get this wrong, then you're not doing it the correct way.
What you're actually performing is you're calculating the $a_n$ and $b_n$ terms for the Fourier transform. Those represent the real and imaginary parts of the frequencies. What you're plotting is the variation of those terms as you sweep a frequency, but the FFT plot is the magnitude.
And your problem is an XY one: you're saying that you get phase errors due to the FFT, but you're not using an FFT, just something that resembles it. At this point I still don't what exactly you're after: to find a particular frequency in a test signal, or to perform an FFT? Because, as shown in the comments, there are no two peaks for multiplying $\sin(t)\cdot\sin(t+\phi)$.
So I'll consider that this is what you want: to find the spectrum of a single frequency from within a test signal of multiple sines, as its phase varies. This is just the basic of the lock-in amplifier/filter linked in the comments, and this is the math of it ($s(t)$ is a sum of sines, using unity amplitudes):
$$\begin{align}
x(t)&=s(t)+\sin(2\pi ft+\phi)&=s(t)+\sin(2\pi\dfrac{t}{T}+\phi) \tag{1} \\
a&=\dfrac{1}{T}\int_0^T{x(t)\sin(2\pi\dfrac{t}{T})\mathrm{d}t}&=\dfrac{\cos(\phi)}{2} \tag{2}\\
b&=\dfrac{1}{T}\int_0^T{x(t)\cos(2\pi\dfrac{t}{T})\mathrm{d}t}&=\dfrac{\sin(\phi)}{2} \tag{3}\\
\end{align}$$
Since you want the peaks for the amplitude, both terms are multiplied by 2. If $x(t)$ has a different amplitude, it will be a constant that comes out of the integration and simply multiply the final result. The magnitude will be $\sqrt{a^2+b^2}$. As for the frequency of that particular peak, it will be given by the second term, $\cos(2t)$:
$$\sin(t)\sin(t+\phi)=\dfrac{\cos(\phi)}{2}-\dfrac{\cos(2t+\phi)}{2} \tag{4}$$
Whether you are perfoming the above steps, or an FFT, the plotted result should be the same. For example, in Octave, for a simple sine with displacement (*n just normalizes the magnitude for the plot):
n = 128;
N = 8192;
wt = 2*pi*[0:n-1]/32;
x = sin(wt + pi/6);
s = sin(wt);
c = cos(wt);
xs = x.*s;
xc = x.*c;
a = mean(xs).*ones(1, N/4);
b = mean(xc).*ones(1, N/4);
X = fft(x, N);
Xs = fft(xs, N);
comp = hypot(a, b)*n;
cosphi = cos(pi/6)/2*ones(1, N/4)*n;
plot(abs(X), "", abs(Xs), "", comp, "", comp/2, "", cosphi)

And the result shows:

X, blue: the spectrum of the original signal, matching...
comp, yellow: the magnitude calculated through the Fourier coefficients
Xc, orange: formula (4), or the spectrum of the signal multiplied with a reference sine (s, resulting xs), matching...
comp/2, purple: half the magnitude of X
cosphi, green: the DC value of xs -- note that it's DC, not a frequency peak

Just in case it's not clear, since you have the real and imaginary parts of the Fourier coefficients, then you have the magnitude (as seen in the code) and the phase as $\mathrm{atan2}(b,a)$. To prove this, too, here is a different perspective:

In the right pane, upper left, you see V(x) which is the frequency that's added to to the other ones. It's not the fundamental, 1 kHz, it's the first harmonic (3rd for only odd ones, 3 kHz). Its amplitude is ⅓ of the fundamental, and its angle is an odd multiple of the base angle which is seen at the top, in .func phi(x) {x*23}. It's the 3rd harmonic, thus it's three times that value. That FM/AM block multiplies the input signal, V(y) (plotted in green on the left), with unity sine and cosine of the desired frequency (3 kHz), averages them over one fundamental period (the whole bottom-left qurter of the schematic), giving the time quantities V(a) and V(b), representing the $a$ and $b$ coefficients. These, in turn, are used to generate the magnitude, V(mag) (black trace), matching the amplitude of the sine, V(x) (blue trace), and the angle, V(phi) (red trace). The magnitude and the angle are measured and their results are displayed in the small blue rectangle at the top (m for magnitude, a for angle). The results are not perfect because of roundings, waveform compression, no imposed timestep, etc, but the results are very representative of what I have been talking about so far: you don't need to have matching amplitudes in order to find the magnitude and the angle, since that is the whole purpose of the $a$ and $b$ coefficients. If this doesn't prove anything, I don't know what else to say, and you're better off waiting for someone else's answer.

Your comments give the impression that you haven't read anything that was said to you. Maybe you don't understand the math behind it (but then how come you said you read the article about the lock-in detection), or the code in Octave (though I kept it very simple, and quite close to what formulas you might see in books).
This is in Python. You used it in your OP, so you should be able to see what's going on. It follows the exact mathematical proof as above:
import numpy as np
t=np.linspace(0,2-1/500,1000)
phi=np.pi/18               # 10 deg
x=np.sin(2*np.pi*t+phi)    # fundamental
for i in range(2,10):      # add 10 harmonics of decreasing amplitude
                           # and increasing phase
  x=x+np.sin(2*np.pi*i*t+i*phi)/i
  i=i+1
n=3                        # check the 3rd harmonic, it should have
                           # 1/3 amplitude, 3*10 deg
s=np.sin(2*np.pi*n*t)      # reference sine, note the n*t
c=np.cos(2*np.pi*n*t)      # reference cosine
xs,xc=x*s,x*c
a,b=2*np.mean(xs),2*np.mean(xc)

At this point you have all you need. If the magnitude is not zero, the signal is there, and you can see what phase it has:
mag=np.hypot(b,a)
rad=np.atan2(b,a)
deg=rad*180/np.pi
mag,rad
(0.33333333333333326, 29.999999999999876)

As you can see, the amplitudes of the signals are not the same (reference sine is 1, 3rd harmonic is ⅓), yet the results are spot on. This means that the method works (not surprisingly).
At this point, the most basic of logic will tell you that, if this method gives the magnitude and the phase of a certain signal then a magnitude of zero will most certainly mean that there is no signal. Sure enough, since there are only 10 harmonics, searching for a 13th harmonic gives these:
s13=np.sin(2*np.pi*13*t)
c13=np.cos(2*np.pi*13*t)
xs13,xc13=x*s13,x*c13
a13,b13=2*np.mean(xs13),2*np.mean(xc13)
a13,b13
(4.369837824924616e-16, 7.993605777301127e-17)

Given your comments, I find it surprising that I have to insist in explaining, since you, yourself, said in the comments:

if that specific frequency doesn't exist, the signal will average out to zero over time

so then why do you ask how to tell if there is a signal when the magnitude is different than zero? Are you serious, or trolling? I hope it's not the latter, because it would mean you're mocking the time people dedicate, out of their own free time, to help others.